Fermi-Hubbard model Hamiltonian and DMRG

The following code uses the Python interface of chemtensor to construct the Fermi-Hubbard Hamiltonian as matrix product operator and to run DMRG.

import numpy as np
import chemtensor

# maximum number of OpenMP threads (0 indicates that OpenMP is not available)
chemtensor.get_max_openmp_threads()

16

Construct Fermi-Hubbard Hamiltonian

# number of spin-endowed lattice sites (local dimension is 4)
nsites = 6

# Hamiltonian parameters
t  = 1.0
u  = 4.0
mu = 1.5

# construct the Fermi-Hubbard Hamiltonian as matrix product operator (MPO)
help(chemtensor.construct_fermi_hubbard_1d_mpo)
hamiltonian = chemtensor.construct_fermi_hubbard_1d_mpo(nsites, t, u, mu)

Help on built-in function construct_fermi_hubbard_1d_mpo in module chemtensor:

construct_fermi_hubbard_1d_mpo(...)
    Construct an MPO representation of the Fermi-Hubbard Hamiltonian with nearest-neighbor hopping on a one-dimensional lattice.
    Syntax: construct_fermi_hubbard_1d_mpo(nsites: int, t: float, u: float, mu: float)

# local dimension
hamiltonian.d

4

# local physical quantum numbers (particle number and spin)
[chemtensor.decode_quantum_number_pair(qnum) for qnum in hamiltonian.qsite]

[(0, 0), (1, -1), (1, 1), (2, 0)]

# virtual bond dimensions
hamiltonian.bond_dims

[1, 6, 6, 6, 6, 6, 1]

Run two-site DMRG

# overall quantum number sector of quantum state (particle number and spin)
q_pnum = 7
q_spin = 1
qnum_sector = chemtensor.encode_quantum_number_pair(q_pnum, q_spin)

# run two-site DMRG
help(chemtensor.dmrg)
psi, en_sweeps, entropy = chemtensor.dmrg(hamiltonian, num_sweeps=4, maxiter_lanczos=25, tol_split=1e-8, max_vdim=32, qnum_sector=qnum_sector)

Help on built-in function dmrg in module chemtensor:

dmrg(...)
    Run the two-site DMRG algorithm for the Hamiltonian provided as MPO.
    Syntax: dmrg(mpo, num_sweeps=5, maxiter_lanczos=20, tol_split=1e-10, max_vdim=256, qnum_sector=0, rng_seed=42)

Evaluate results and compare with exact diagonalization

# virtual bond dimensions of optimized MPS
psi.bond_dims

[1, 4, 16, 30, 16, 4, 1]

# splitting entropies for each bond
entropy

[1.0572974362766947,
 1.0833061756872846,
 1.1618575985077486,
 1.0833061757321103,
 1.0572974362767673]

# represent 'psi' as vector
psi_vec = psi.to_statevector()
psi_vec.shape

(4096,)

# must be normalized
np.linalg.norm(psi_vec)

0.9999999999999998

# energy after each DMRG sweep
en_sweeps

[-18.48435890403331,
 -18.484358904033307,
 -18.484358904033314,
 -18.4843589040333]

# construct the (dense) matrix representation of the matrix product operator on the full Hilbert space
h_mat = hamiltonian.to_matrix()
# Hilbert space dimension is 4^nsites
h_mat.shape

(4096, 4096)

# check consistency with energy expectation value (difference should be numerically zero):
abs(np.vdot(psi_vec, h_mat @ psi_vec) - en_sweeps[-1])

3.552713678800501e-15

# reference eigenvalues (based on exact diagonalization of matrix representation)
w_ref = np.linalg.eigvalsh(h_mat)
# reference ground state energy
w_ref[0]

-18.484358962762183

# difference to reference ground state energy
en_sweeps - w_ref[0]

array([5.87288724e-08, 5.87288760e-08, 5.87288689e-08, 5.87288831e-08])