renaming subfolder "mps" -> "state", and additional minor revisions

CMakeLists.txt

@@ -10,9 +10,9 @@ find_package(LAPACK REQUIRED)
 find_package(HDF5 REQUIRED COMPONENTS C)
 find_package(Python3 REQUIRED COMPONENTS Development NumPy)
 
-set(CHEMTENSOR_DIRS "src" "src/tensor" "src/mps" "src/operator" "src/algorithm" "src/util")
-set(CHEMTENSOR_SOURCES "src/tensor/dense_tensor.c" "src/tensor/block_sparse_tensor.c" "src/tensor/qnumber.c" "src/tensor/clebsch_gordan.c" "src/tensor/su2_recoupling.c" "src/tensor/su2_tree.c" "src/tensor/su2_tensor.c" "src/mps/mps.c" "src/operator/op_chain.c" "src/operator/local_op.c" "src/operator/mpo_graph.c" "src/operator/mpo.c" "src/operator/ttno_graph.c" "src/operator/ttno.c" "src/operator/hamiltonian.c" "src/algorithm/bond_ops.c" "src/algorithm/operation.c" "src/algorithm/dmrg.c" "src/algorithm/gradient.c" "src/util/util.c" "src/util/queue.c" "src/util/linked_list.c" "src/util/hash_table.c" "src/util/abstract_graph.c" "src/util/bipartite_graph.c" "src/util/integer_linear_algebra.c" "src/util/krylov.c" "src/util/pcg_basic.c" "src/util/rng.c")
-set(TEST_SOURCES "test/tensor/test_dense_tensor.c" "test/tensor/test_block_sparse_tensor.c" "test/tensor/test_clebsch_gordan.c" "test/tensor/test_su2_tree.c" "test/tensor/test_su2_tensor.c" "test/mps/test_mps.c" "test/operator/test_mpo_graph.c" "test/operator/test_mpo.c" "test/operator/test_ttno_graph.c" "test/operator/test_ttno.c" "test/operator/test_hamiltonian.c" "test/algorithm/test_bond_ops.c" "test/algorithm/test_operation.c" "test/algorithm/test_dmrg.c" "test/algorithm/numerical_gradient.c" "test/algorithm/test_gradient.c" "test/util/test_queue.c" "test/util/test_linked_list.c" "test/util/test_hash_table.c" "test/util/test_bipartite_graph.c" "test/util/test_integer_linear_algebra.c" "test/util/test_krylov.c" "test/run_tests.c")
+set(CHEMTENSOR_DIRS "src" "src/tensor" "src/state" "src/operator" "src/algorithm" "src/util")
+set(CHEMTENSOR_SOURCES "src/tensor/dense_tensor.c" "src/tensor/block_sparse_tensor.c" "src/tensor/qnumber.c" "src/tensor/clebsch_gordan.c" "src/tensor/su2_recoupling.c" "src/tensor/su2_tree.c" "src/tensor/su2_tensor.c" "src/state/mps.c" "src/operator/op_chain.c" "src/operator/local_op.c" "src/operator/mpo_graph.c" "src/operator/mpo.c" "src/operator/ttno_graph.c" "src/operator/ttno.c" "src/operator/hamiltonian.c" "src/algorithm/bond_ops.c" "src/algorithm/operation.c" "src/algorithm/dmrg.c" "src/algorithm/gradient.c" "src/util/util.c" "src/util/queue.c" "src/util/linked_list.c" "src/util/hash_table.c" "src/util/abstract_graph.c" "src/util/bipartite_graph.c" "src/util/integer_linear_algebra.c" "src/util/krylov.c" "src/util/pcg_basic.c" "src/util/rng.c")
+set(TEST_SOURCES "test/tensor/test_dense_tensor.c" "test/tensor/test_block_sparse_tensor.c" "test/tensor/test_clebsch_gordan.c" "test/tensor/test_su2_tree.c" "test/tensor/test_su2_tensor.c" "test/state/test_mps.c" "test/operator/test_mpo_graph.c" "test/operator/test_mpo.c" "test/operator/test_ttno_graph.c" "test/operator/test_ttno.c" "test/operator/test_hamiltonian.c" "test/algorithm/test_bond_ops.c" "test/algorithm/test_operation.c" "test/algorithm/test_dmrg.c" "test/algorithm/numerical_gradient.c" "test/algorithm/test_gradient.c" "test/util/test_queue.c" "test/util/test_linked_list.c" "test/util/test_hash_table.c" "test/util/test_bipartite_graph.c" "test/util/test_integer_linear_algebra.c" "test/util/test_krylov.c" "test/run_tests.c")
 
 add_executable(            chemtensor_test ${CHEMTENSOR_SOURCES} ${TEST_SOURCES})
 target_include_directories(chemtensor_test PRIVATE ${CHEMTENSOR_DIRS} ${BLAS_INCLUDE_DIRS} ${LAPACKE_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS})

src/operator/hamiltonian.c

@@ -2693,7 +2693,7 @@ static void spin_molecular_mpo_graph_add_term(const struct spin_molecular_mpo_gr
 /// \brief Construct a molecular Hamiltonian as MPO assembly, assuming a spin orbital basis and
 /// using physicists' convention for the interaction term (note ordering of k and l):
 /// \f[
-/// H = \sum_{i,j,\sigma} t_{i,j} a^{\dagger}_{i,\sigma} a_{j,\sigma} + \\frac{1}{2} \sum_{i,j,k,\ell,\sigma,\tau} v_{i,j,k,\ell} a^{\dagger}_{i,\sigma} a^{\dagger}_{j,\tau} a_{\ell,\tau} a_{k,\sigma}
+/// H = \sum_{i,j,\sigma} t_{i,j} a^{\dagger}_{i,\sigma} a_{j,\sigma} + \frac{1}{2} \sum_{i,j,k,\ell,\sigma,\tau} v_{i,j,k,\ell} a^{\dagger}_{i,\sigma} a^{\dagger}_{j,\tau} a_{\ell,\tau} a_{k,\sigma}
 /// \f]
 ///
 /// If 'optimize == true', optimize the virtual bond dimensions via the automatic construction starting from operator chains.

src/mps/mps.c → src/state/mps.c

@@ -318,95 +318,85 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 
 //________________________________________________________________________________________________________________________
 ///
-/// \brief Compute the addition of two MPS chi and psi.
+/// \brief Compute the logical addition of two MPS `chi` and `psi` (summing their virtual bond dimensions).
 ///
 void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret)
 {
 	// number of lattice sites must agree
 	assert(chi->nsites == psi->nsites);
-
 	// number of lattice sites must be larger than 0
 	assert(chi->nsites > 0);
 
-	// physical quantum numbers must agree
-	assert(chi->d == psi->d); 
-	assert(qnumber_all_equal(chi->d, chi->qsite, psi->qsite)); 
-
-	const long d = chi->d;
-	const int L = chi->nsites;
-	const enum numeric_type dtype = chi->a[0].dtype;
-
-	// initialize return mps
-	allocate_empty_mps(L, d, chi->qsite, ret);
-
-	if (L == 1) {		
-		const struct block_sparse_tensor chi_a = chi->a[0];
-		const struct block_sparse_tensor psi_a = psi->a[0];
-		const int ndim = chi_a.ndim;
-
-		assert(chi_a.ndim == psi_a.ndim);
-
-		// dummy bond quantum numbers must agree
-		assert(qnumber_all_equal(
-				chi_a.dim_logical[0], 
-				chi_a.qnums_blocks[0], 
-				psi_a.qnums_blocks[0])); 
-		assert(qnumber_all_equal(
-				chi_a.dim_logical[2],
-				chi_a.qnums_blocks[2],
-				psi_a.qnums_blocks[2]));
+	const int nsites = chi->nsites;
 
+	// physical quantum numbers must agree
+	assert(chi->d == psi->d);
+	assert(qnumber_all_equal(chi->d, chi->qsite, psi->qsite));
+
+	// leading and trailing (dummy) bond quantum numbers must agree
+	assert(chi->a[0].dim_logical[0] ==
+	       psi->a[0].dim_logical[0]);
+	assert(qnumber_all_equal(
+			chi->a[0].dim_logical[0],
+			chi->a[0].qnums_blocks[0],
+			psi->a[0].qnums_blocks[0]));
+	assert(chi->a[nsites - 1].dim_logical[2] ==
+	       psi->a[nsites - 1].dim_logical[2]);
+	assert(qnumber_all_equal(
+			chi->a[nsites - 1].dim_logical[2],
+			chi->a[nsites - 1].qnums_blocks[2],
+			psi->a[nsites - 1].qnums_blocks[2]));
+
+	// initialize return MPS
+	allocate_empty_mps(nsites, chi->d, chi->qsite, ret);
+
+	if (nsites == 1)
+	{
 		// copy sparse tensor into resulting tensor
-		copy_block_sparse_tensor(&chi_a, &ret->a[0]);
+		copy_block_sparse_tensor(&chi->a[0], &ret->a[0]);
 
-		// add individual dense_tensors
-		const long nblocks = integer_product(ret->a->dim_blocks, ndim);
+		// add individual dense tensors
+		const long nblocks = integer_product(ret->a[0].dim_blocks, ret->a[0].ndim);
 		for (long k = 0; k < nblocks; k++)
 		{
-			struct dense_tensor* a = psi_a.blocks[k];
-			struct dense_tensor* b = ret->a->blocks[k];
-			if (a != NULL && b != NULL) {
-				dense_tensor_scalar_multiply_add(numeric_one(dtype), a, b);
+			struct dense_tensor* a = psi->a[0].blocks[k];
+			struct dense_tensor* b = ret->a[0].blocks[k];
+			if (a != NULL) {
+				assert(b != NULL);
+				dense_tensor_scalar_multiply_add(numeric_one(a->dtype), a, b);
 			}
 		}
-	} else {
-		// leading and trailing (dummy) bond quantum numbers must agree
-		assert(qnumber_all_equal(
-				chi->a[0].dim_logical[0], 
-				chi->a[0].qnums_blocks[0], 
-				psi->a[0].qnums_blocks[0]));
-		assert(qnumber_all_equal(
-				chi->a[chi->nsites - 1].dim_logical[2], 
-				chi->a[chi->nsites - 1].qnums_blocks[2], 
-				psi->a[chi->nsites - 1].qnums_blocks[2])); 
-
-		struct block_sparse_tensor tlist[2];
-
+	}
+	else  // nsites > 1
+	{
 		// left-most tensor
 		{
 			const int i_ax[1] = { 2 };
-			tlist[0] = chi->a[0];
-			tlist[1] = psi->a[0];
-
+			struct block_sparse_tensor tlist[2] = {
+				chi->a[0],
+				psi->a[0],
+			};
 			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret->a[0]);
 		}
 
 		// intermediate tensors
-		for (int i = 1; i < L - 1; i++) {
+		for (int i = 1; i < nsites - 1; i++) {
 			const int i_ax[2] = { 0, 2 };
-			tlist[0] = chi->a[i];
-			tlist[1] = psi->a[i];
-
+			struct block_sparse_tensor tlist[2] = {
+				chi->a[i],
+				psi->a[i],
+			};
 			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &ret->a[i]);
 		}
 
 		// right-most tensor
 		{
 			const int i_ax[1] = { 0 };
-			tlist[0] = chi->a[L - 1];
-			tlist[1] = psi->a[L - 1];
-
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret->a[L - 1]);
+			struct block_sparse_tensor tlist[2] = {
+				chi->a[nsites - 1],
+				psi->a[nsites - 1],
+			};
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret->a[nsites - 1]);
 		}
 	}
 }

src/mps/mps.h → src/state/mps.h

@@ -14,8 +14,8 @@
 struct mps
 {
 	struct block_sparse_tensor* a;  //!< tensors associated with sites, with dimensions D_i x d x D_{i+1}; array of length 'nsites'
-	long d;                         //!< local physical dimension of each site
 	qnumber* qsite;                 //!< physical quantum numbers at each site
+	long d;                         //!< local physical dimension of each site
 	int nsites;                     //!< number of sites
 };
 

test/mps/test_mps.c → test/state/test_mps.c

@@ -4,9 +4,12 @@
 #include "aligned_memory.h"
 
 
+#define ARRLEN(a) (sizeof(a) / sizeof(a[0]))
+
+
 char* test_mps_vdot()
 {
-	hid_t file = H5Fopen("../test/mps/data/test_mps_vdot.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
+	hid_t file = H5Fopen("../test/state/data/test_mps_vdot.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
 	if (file < 0) {
 		return "'H5Fopen' in test_mps_vdot failed";
 	}
@@ -129,33 +132,34 @@ char* test_mps_vdot()
 
 char* test_mps_add()
 {
-	const long d = 2, max_vdim = 16;
-	const int site_array[6] = {1, 3, 5, 8, 12, 16};
-
-	const qnumber qnum_sector = 8;
-	const qnumber qsite_chi[2] = {1, 2};
-	const qnumber qsite_psi[2] = {1, 2};
-
-	struct mps chi, psi, res;
-	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
-	struct dense_tensor chi_vec_dns, psi_vec_dns, res_vec_dns;
-
+	const long d = 3;
+	const qnumber qsite[3] = { 1, 0, -2 };
+
+	const int nsites_list[] = { 1, 3, 5, 8, 12 };
+
 	struct rng_state rng_state;
 	seed_rng_state(42, &rng_state);
-
-	int nsites;
-	for (int i = 0; i < 6; i++) {
-		nsites = site_array[i];
 
-		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_chi, qnum_sector, max_vdim, &rng_state, &chi);
-		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_psi, qnum_sector, max_vdim, &rng_state, &psi);
+	for (int i = 0; i < (int)ARRLEN(nsites_list); i++)
+	{
+		const int nsites = nsites_list[i];
+
+		struct mps chi, psi;
+		const qnumber qnum_sector = (i % 4);
+		const long max_vdim = 16;
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &chi);
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &psi);
 
+		// perform logical addition
+		struct mps res;
 		mps_add(&chi, &psi, &res);
 
+		struct block_sparse_tensor chi_vec, psi_vec, res_vec;
 		mps_to_statevector(&chi, &chi_vec);
 		mps_to_statevector(&psi, &psi_vec);
 		mps_to_statevector(&res, &res_vec);
 
+		struct dense_tensor chi_vec_dns, psi_vec_dns, res_vec_dns;
 		block_sparse_to_dense_tensor(&chi_vec, &chi_vec_dns);
 		block_sparse_to_dense_tensor(&psi_vec, &psi_vec_dns);
 		block_sparse_to_dense_tensor(&res_vec, &res_vec_dns);
@@ -164,7 +168,7 @@ char* test_mps_add()
 
 		// compare dense tensors
 		if (!dense_tensor_allclose(&res_vec_dns, &chi_vec_dns, 1e-13)) {
-			return "addition of mps does not match reference";
+			return "logical addition of two MPS does not match reference";
 		}
 
 		// free memory
@@ -187,7 +191,7 @@ char* test_mps_add()
 
 char* test_mps_orthonormalize_qr()
 {
-	hid_t file = H5Fopen("../test/mps/data/test_mps_orthonormalize_qr.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
+	hid_t file = H5Fopen("../test/state/data/test_mps_orthonormalize_qr.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
 	if (file < 0) {
 		return "'H5Fopen' in test_mps_orthonormalize_qr failed";
 	}
@@ -321,7 +325,7 @@ char* test_mps_orthonormalize_qr()
 
 char* test_mps_compress()
 {
-	hid_t file = H5Fopen("../test/mps/data/test_mps_compress.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
+	hid_t file = H5Fopen("../test/state/data/test_mps_compress.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
 	if (file < 0) {
 		return "'H5Fopen' in test_mps_compress failed";
 	}
@@ -472,7 +476,7 @@ char* test_mps_compress()
 
 char* test_mps_split_tensor_svd()
 {
-	hid_t file = H5Fopen("../test/mps/data/test_mps_split_tensor_svd.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
+	hid_t file = H5Fopen("../test/state/data/test_mps_split_tensor_svd.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
 	if (file < 0) {
 		return "'H5Fopen' in test_mps_split_tensor_svd failed";
 	}
@@ -605,7 +609,7 @@ char* test_mps_split_tensor_svd()
 
 char* test_mps_to_statevector()
 {
-	hid_t file = H5Fopen("../test/mps/data/test_mps_to_statevector.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
+	hid_t file = H5Fopen("../test/state/data/test_mps_to_statevector.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
 	if (file < 0) {
 		return "'H5Fopen' in test_mps_to_statevector failed";
 	}