basic DMRG demo applied to the Fermi-Hubbard model

CMakeLists.txt

@@ -13,9 +13,12 @@ set(CHEMTENSOR_DIRS "src" "src/tensor" "src/mps" "src/operator" "src/algorithm"
 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")
 
-add_executable(chemtensor_test ${CHEMTENSOR_SOURCES} ${TEST_SOURCES})
-
+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})
 target_link_libraries(     chemtensor_test PRIVATE ${BLAS_LIBRARIES} ${LAPACKE_LIBRARIES} lapacke ${HDF5_LIBRARIES})
 
+add_executable(            basic_dmrg_fermi_hubbard ${CHEMTENSOR_SOURCES} "examples/dmrg/basic_dmrg_fermi_hubbard.c")
+target_include_directories(basic_dmrg_fermi_hubbard PRIVATE ${CHEMTENSOR_DIRS} ${BLAS_INCLUDE_DIRS} ${LAPACKE_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS})
+target_link_libraries(     basic_dmrg_fermi_hubbard PRIVATE ${BLAS_LIBRARIES} ${LAPACKE_LIBRARIES} lapacke ${HDF5_LIBRARIES})
+
 add_test(NAME chemtensor_test COMMAND chemtensor_test)

examples/dmrg/basic_dmrg_fermi_hubbard.c

@@ -0,0 +1,124 @@
+#include <cblas.h>
+#include <lapacke.h>
+#include <stdio.h>
+#include "hamiltonian.h"
+#include "dmrg.h"
+#include "aligned_memory.h"
+#include "rng.h"
+#include "util.h"
+
+
+int main()
+{
+	// number of spin-endowed lattice sites (local dimension is 4)
+	const int nsites = 6;
+	printf("number of spin-endowed lattice sites: %i -> Hilbert space dimension: %li\n", nsites, ipow(4, nsites));
+
+	// Hamiltonian parameters
+	const double t  = 1.0;
+	const double u  = 4.0;
+	const double mu = 1.5;
+	printf("Fermi-Hubbard Hamiltonian parameters: t = %g, u = %g, mu = %g\n", t, u, mu);
+
+	// overall quantum number sector of quantum state (particle number and spin)
+	const qnumber q_pnum = 7;
+	const qnumber q_spin = 1;
+	printf("psi restricted to overall quantum number sector with particle number %i and spin %i (integer convention)\n", q_pnum, q_spin);
+	const qnumber qnum_sector = fermi_hubbard_encode_quantum_numbers(q_pnum, q_spin);
+
+	// maximum allowed MPS bond dimension
+	const long max_vdim = 32;
+
+	struct mpo hamiltonian;
+	{
+		struct mpo_assembly assembly;
+		construct_fermi_hubbard_1d_mpo_assembly(nsites, t, u, mu, &assembly);
+		mpo_from_assembly(&assembly, &hamiltonian);
+		delete_mpo_assembly(&assembly);
+	}
+	if (!mpo_is_consistent(&hamiltonian)) {
+		fprintf(stderr, "internal consistency check for Fermi-Hubbard Hamiltonian MPO failed");
+		return -1;
+	}
+
+	// initial state vector as MPS
+	struct mps psi;
+	{
+		struct rng_state rng_state;
+		seed_rng_state(42, &rng_state);
+
+		construct_random_mps(hamiltonian.a[0].dtype, hamiltonian.nsites, hamiltonian.d, hamiltonian.qsite, qnum_sector, max_vdim, &rng_state, &psi);
+		if (!mps_is_consistent(&psi)) {
+			fprintf(stderr, "internal MPS consistency check failed");
+			return -1;
+		}
+	}
+
+	// run two-site DMRG
+	const int num_sweeps      = 4;
+	const int maxiter_lanczos = 25;
+	double tol_split          = 1e-8;
+	double* en_sweeps = aligned_alloc(MEM_DATA_ALIGN, num_sweeps * sizeof(double));
+	double* entropy   = aligned_alloc(MEM_DATA_ALIGN, (nsites - 1) * sizeof(double));
+	printf("Running two-site DMRG (using num_sweeps = %i, maxiter_lanczos = %i, tol_split = %g)... ", num_sweeps, maxiter_lanczos, tol_split);
+	if (dmrg_twosite(&hamiltonian, num_sweeps, maxiter_lanczos, tol_split, max_vdim, &psi, en_sweeps, entropy) < 0) {
+		fprintf(stderr, "'dmrg_twosite' failed internally");
+		return -2;
+	}
+	printf("Done.\n");
+
+	printf("||psi||: %g\n", mps_norm(&psi));
+
+	printf("psi bond dimensions: ");
+	for (int l = 0; l < nsites + 1; l++) {
+		printf("%li ", mps_bond_dim(&psi, l));
+	}
+	printf("\n");
+	printf("splitting entropies for each bond: ");
+	for (int l = 0; l < nsites - 1; l++) {
+		printf("%g ", entropy[l]);
+	}
+	printf("\n");
+
+	// matrix representation of Hamiltonian
+	struct dense_tensor hmat;
+	{
+		struct block_sparse_tensor hsparse;
+		mpo_to_matrix(&hamiltonian, &hsparse);
+
+		// convert to dense tensor
+		block_sparse_to_dense_tensor(&hsparse, &hmat);
+		delete_block_sparse_tensor(&hsparse);
+
+		// dimension and data type consistency checks
+		assert(hmat.ndim == 4);
+		// dummy virtual bond dimensions are retained
+		assert(hmat.dim[0] == 1 && hmat.dim[3] == 1);
+		assert(hmat.dim[1] == ipow(hamiltonian.d, hamiltonian.nsites));
+		assert(hmat.dim[1] == hmat.dim[2]);
+		assert(hmat.dtype == DOUBLE_REAL);
+	}
+
+	// reference eigenvalues (based on exact diagonalization of matrix representation)
+	double* w_ref = aligned_alloc(MEM_DATA_ALIGN, hmat.dim[1] * sizeof(double));
+	int info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', hmat.dim[1], hmat.data, hmat.dim[2], w_ref);
+	if (info != 0) {
+		fprintf(stderr, "LAPACK function 'dsyev' failed, return value: %i\n", info);
+		return info;
+	}
+	printf("reference ground state energy: %g\n", w_ref[0]);
+
+	printf("energy after each DMRG sweep:\n");
+	for (int i = 0; i < num_sweeps; i++) {
+		printf("%i: %g, diff to reference: %g\n", i + 1, en_sweeps[i], en_sweeps[i] - w_ref[0]);
+	}
+
+	aligned_free(w_ref);
+	delete_dense_tensor(&hmat);
+	aligned_free(entropy);
+	aligned_free(en_sweeps);
+	delete_mps(&psi);
+	delete_mpo(&hamiltonian);
+
+	return 0;
+}

src/algorithm/dmrg.c

@@ -246,7 +246,7 @@ int dmrg_singlesite(const struct mpo* hamiltonian, const int num_sweeps, const i
 /// The input 'psi' is used as starting state and is updated in-place during the optimization.
 ///
 int dmrg_twosite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, const double tol_split, const long max_vdim,
-	struct mps* psi, double* en_sweeps, double* restrict entropy)
+	struct mps* psi, double* restrict en_sweeps, double* restrict entropy)
 {
 	// number of lattice sites
 	const int nsites = hamiltonian->nsites;

src/algorithm/dmrg.h

@@ -10,4 +10,4 @@
 int dmrg_singlesite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, struct mps* psi, double* en_sweeps);
 
 int dmrg_twosite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, const double tol_split, const long max_vdim,
-	struct mps* psi, double* en_sweeps, double* restrict entropy);
+	struct mps* psi, double* restrict en_sweeps, double* restrict entropy);

src/mps/mps.c

@@ -55,6 +55,97 @@ void delete_mps(struct mps* mps)
 }
 
 
+//________________________________________________________________________________________________________________________
+///
+/// \brief Construct a matrix product state with random normal tensor entries, given an overall quantum number sector and maximum virtual bond dimension.
+///
+void construct_random_mps(const enum numeric_type dtype, const int nsites, const long d, const qnumber* qsite, const qnumber qnum_sector, const long max_vdim, struct rng_state* rng_state, struct mps* mps)
+{
+	assert(nsites >= 1);
+	assert(d >= 1);
+
+	// virtual bond quantum numbers
+	long* dim_bonds  = aligned_alloc(MEM_DATA_ALIGN, (nsites + 1) * sizeof(long));
+	qnumber** qbonds = aligned_alloc(MEM_DATA_ALIGN, (nsites + 1) * sizeof(qnumber*));
+	// dummy left virtual bond; set quantum number to zero
+	dim_bonds[0] = 1;
+	qbonds[0] = aligned_alloc(MEM_DATA_ALIGN, sizeof(qnumber));
+	qbonds[0][0] = 0;
+	// dummy right virtual bond; set quantum number to overall quantum number sector
+	dim_bonds[nsites] = 1;
+	qbonds[nsites] = aligned_alloc(MEM_DATA_ALIGN, sizeof(qnumber));
+	qbonds[nsites][0] = qnum_sector;
+	// virtual bond quantum numbers on left half
+	for (int l = 1; l < (nsites + 1) / 2; l++)
+	{
+		// enumerate all combinations of left bond quantum numbers and local physical quantum numbers
+		const long dim_full = dim_bonds[l - 1] * d;
+		qnumber* qnums_full = aligned_alloc(MEM_DATA_ALIGN, dim_full * sizeof(qnumber));
+		for (long i = 0; i < dim_bonds[l - 1]; i++) {
+			for (long j = 0; j < d; j++) {
+				qnums_full[i*d + j] = qbonds[l - 1][i] + qsite[j];
+			}
+		}
+		dim_bonds[l] = lmin(dim_full, max_vdim);
+		qbonds[l] = aligned_alloc(MEM_DATA_ALIGN, dim_bonds[l] * sizeof(qnumber));
+		if (dim_full <= max_vdim) {
+			memcpy(qbonds[l], qnums_full, dim_bonds[l] * sizeof(qnumber));
+		}
+		else {
+			// randomly select quantum numbers
+			for (long i = 0; i < dim_bonds[l]; i++) {
+				qbonds[l][i] = qnums_full[rand_interval(dim_full, rng_state)];
+			}
+		}
+		aligned_free(qnums_full);
+	}
+	// virtual bond quantum numbers on right half
+	for (int l = nsites - 1; l >= (nsites + 1) / 2; l--)
+	{
+		// enumerate all combinations of right bond quantum numbers and local physical quantum numbers
+		const long dim_full = dim_bonds[l + 1] * d;
+		qnumber* qnums_full = aligned_alloc(MEM_DATA_ALIGN, dim_full * sizeof(qnumber));
+		for (long i = 0; i < dim_bonds[l + 1]; i++) {
+			for (long j = 0; j < d; j++) {
+				qnums_full[i*d + j] = qbonds[l + 1][i] - qsite[j];
+			}
+		}
+		dim_bonds[l] = lmin(dim_full, max_vdim);
+		qbonds[l] = aligned_alloc(MEM_DATA_ALIGN, dim_bonds[l] * sizeof(qnumber));
+		if (dim_full <= max_vdim) {
+			memcpy(qbonds[l], qnums_full, dim_bonds[l] * sizeof(qnumber));
+		}
+		else {
+			// randomly select quantum numbers
+			for (long i = 0; i < dim_bonds[l]; i++) {
+				qbonds[l][i] = qnums_full[rand_interval(dim_full, rng_state)];
+			}
+		}
+		aligned_free(qnums_full);
+	}
+
+	allocate_mps(dtype, nsites, d, qsite, dim_bonds, (const qnumber**)qbonds, mps);
+
+	for (int l = 0; l < nsites + 1; l++) {
+		aligned_free(qbonds[l]);
+	}
+	aligned_free(qbonds);
+	aligned_free(dim_bonds);
+
+	// fill MPS tensor entries with pseudo-random numbers, scaled by 1 / sqrt("number of entries")
+	for (int l = 0; l < nsites; l++)
+	{
+		// logical number of entries in MPS tensor
+		const long nelem = integer_product(mps->a[l].dim_logical, mps->a[l].ndim);
+		// ensure that 'alpha' is large enough to store any numeric type
+		dcomplex alpha;
+		assert(mps->a[l].dtype == dtype);
+		numeric_from_double(1.0 / sqrt(nelem), dtype, &alpha);
+		block_sparse_tensor_fill_random_normal(&alpha, numeric_zero(dtype), rng_state, &mps->a[l]);
+	}
+}
+
+
 //________________________________________________________________________________________________________________________
 ///
 /// \brief Internal consistency check of the MPS data structure.
@@ -318,7 +409,7 @@ void mps_local_orthonormalize_qr(struct block_sparse_tensor* restrict a, struct
 
 //________________________________________________________________________________________________________________________
 ///
-/// \brief Right-orthonormalize a local MPS site tensor by RQ decomposition, and update tensor at next site.
+/// \brief Right-orthonormalize a local MPS site tensor by RQ decomposition, and update tensor at previous site.
 ///
 void mps_local_orthonormalize_rq(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_prev)
 {

src/mps/mps.h

@@ -29,6 +29,8 @@ void allocate_mps(const enum numeric_type dtype, const int nsites, const long d,
 
 void delete_mps(struct mps* mps);
 
+void construct_random_mps(const enum numeric_type dtype, const int nsites, const long d, const qnumber* qsite, const qnumber qnum_sector, const long max_vdim, struct rng_state* rng_state, struct mps* mps);
+
 bool mps_is_consistent(const struct mps* mps);
 
 

src/operator/hamiltonian.c

@@ -14,7 +14,7 @@
 
 //________________________________________________________________________________________________________________________
 ///
-/// \brief Construct a Hamiltonian as MPO based on local operator chains, which are shifted along a 1D lattice.
+/// \brief Construct a Hamiltonian as MPO operator graph based on local operator chains, which are shifted along a 1D lattice.
 ///
 static void local_opchains_to_mpo_graph(const int nsites, const struct op_chain* lopchains, const int nlopchains, struct mpo_graph* graph)
 {
@@ -244,10 +244,10 @@ void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, c
 	const qnumber qs[4] = { 0, -1,  1,  0 };
 	assembly->d = 4;
 	assembly->qsite = aligned_alloc(MEM_DATA_ALIGN, assembly->d * sizeof(qnumber));
-	assembly->qsite[0] = (qn[0] << 16) + qs[0];
-	assembly->qsite[1] = (qn[1] << 16) + qs[1];
-	assembly->qsite[2] = (qn[2] << 16) + qs[2];
-	assembly->qsite[3] = (qn[3] << 16) + qs[3];
+	assembly->qsite[0] = fermi_hubbard_encode_quantum_numbers(qn[0], qs[0]);
+	assembly->qsite[1] = fermi_hubbard_encode_quantum_numbers(qn[1], qs[1]);
+	assembly->qsite[2] = fermi_hubbard_encode_quantum_numbers(qn[2], qs[2]);
+	assembly->qsite[3] = fermi_hubbard_encode_quantum_numbers(qn[3], qs[3]);
 
 	assembly->dtype = DOUBLE_REAL;
 
@@ -343,16 +343,13 @@ void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, c
 	assembly->coeffmap = aligned_alloc(MEM_DATA_ALIGN, sizeof(coeffmap));
 	memcpy(assembly->coeffmap, coeffmap, sizeof(coeffmap));
 
-	// cast to unsigned integer to avoid compiler warning when bit-shifting
-	const unsigned int n1 = -1;
-
 	// local two-site and single-site terms
 	// spin-up kinetic hopping
-	int oids_c0[] = { 3, 2 };  qnumber qnums_c0[] = { 0, ( 1 << 16) + 1, 0 };
-	int oids_c1[] = { 4, 1 };  qnumber qnums_c1[] = { 0, (n1 << 16) - 1, 0 };
+	int oids_c0[] = { 3, 2 };  qnumber qnums_c0[] = { 0, fermi_hubbard_encode_quantum_numbers( 1,  1), 0 };
+	int oids_c1[] = { 4, 1 };  qnumber qnums_c1[] = { 0, fermi_hubbard_encode_quantum_numbers(-1, -1), 0 };
 	// spin-down kinetic hopping
-	int oids_c2[] = { 5, 8 };  qnumber qnums_c2[] = { 0, ( 1 << 16) - 1, 0 };
-	int oids_c3[] = { 6, 7 };  qnumber qnums_c3[] = { 0, (n1 << 16) + 1, 0 };
+	int oids_c2[] = { 5, 8 };  qnumber qnums_c2[] = { 0, fermi_hubbard_encode_quantum_numbers( 1, -1), 0 };
+	int oids_c3[] = { 6, 7 };  qnumber qnums_c3[] = { 0, fermi_hubbard_encode_quantum_numbers(-1,  1), 0 };
 	// number operator - mu (n_up + n_dn) and interaction u (n_up-1/2) (n_dn-1/2)
 	int oids_c4[] = {  9 };    qnumber qnums_c4[] = { 0, 0 };
 	int oids_c5[] = { 10 };    qnumber qnums_c5[] = { 0, 0 };

src/operator/hamiltonian.h

@@ -12,6 +12,15 @@ void construct_heisenberg_xxz_1d_mpo_assembly(const int nsites, const double J,
 
 void construct_bose_hubbard_1d_mpo_assembly(const int nsites, const long d, const double t, const double u, const double mu, struct mpo_assembly* assembly);
 
+//________________________________________________________________________________________________________________________
+///
+/// \brief Encode a particle number and spin quantum number for the Fermi-Hubbard Hamiltonian into a single quantum number.
+///
+static inline qnumber fermi_hubbard_encode_quantum_numbers(const qnumber q_pnum, const qnumber q_spin)
+{
+	return (q_pnum << 16) + q_spin;
+}
+
 void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, const double u, const double mu, struct mpo_assembly* assembly);
 
 void construct_molecular_hamiltonian_mpo_assembly(const struct dense_tensor* restrict tkin, const struct dense_tensor* restrict vint, struct mpo_assembly* assembly);

src/tensor/block_sparse_tensor.c

@@ -1744,16 +1744,6 @@ void block_sparse_tensor_dot(const struct block_sparse_tensor* restrict s, const
 }
 
 
-//________________________________________________________________________________________________________________________
-///
-/// \brief Minimum of two integers.
-///
-static inline long minl(const long a, const long b)
-{
-	return (a <= b) ? a : b;
-}
-
-
 //________________________________________________________________________________________________________________________
 ///
 /// \brief Compute the logical QR decomposition of a block-sparse matrix.
@@ -1790,7 +1780,7 @@ int block_sparse_tensor_qr(const struct block_sparse_tensor* restrict a, struct
 			assert(b->ndim == 2);
 			assert(b->dtype == a->dtype);
 
-			const long k = minl(b->dim[0], b->dim[1]);
+			const long k = lmin(b->dim[0], b->dim[1]);
 
 			// append a sequence of logical quantum numbers of length k
 			for (long l = 0; l < k; l++) {
@@ -1929,7 +1919,7 @@ int block_sparse_tensor_rq(const struct block_sparse_tensor* restrict a, struct
 			assert(b->ndim == 2);
 			assert(b->dtype == a->dtype);
 
-			const long k = minl(b->dim[0], b->dim[1]);
+			const long k = lmin(b->dim[0], b->dim[1]);
 
 			// append a sequence of logical quantum numbers of length k
 			for (long l = 0; l < k; l++) {
@@ -2068,7 +2058,7 @@ int block_sparse_tensor_svd(const struct block_sparse_tensor* restrict a, struct
 			assert(b->ndim == 2);
 			assert(b->dtype == a->dtype);
 
-			const long k = minl(b->dim[0], b->dim[1]);
+			const long k = lmin(b->dim[0], b->dim[1]);
 
 			// append a sequence of logical quantum numbers of length k
 			for (long l = 0; l < k; l++) {