3. DMRG algorithm¶

This tutorial demonstrates an application of the DMRG algorithm to find the ground state of the Fermi-Hubbard model on a one-dimensional lattice.

[1]:

# import NumPy and the PyTeNet package
import numpy as np
import pytenet as ptn

3.1. Construct the Fermi-Hubbard Hamiltonian¶

[2]:

# 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 a matrix product operator (MPO)
hamiltonian = ptn.fermi_hubbard_1d_mpo(nsites, t, u, mu)

[3]:

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

[3]:

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

[4]:

# virtual bond dimensions
hamiltonian.bond_dims

[4]:

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

3.2. Run two-site DMRG¶

[5]:

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

[6]:

# random initial MPS
rng = np.random.default_rng(42)
psi = ptn.MPS.construct_random(nsites, hamiltonian.qsite, qnum_sector, max_vdim=18, dtype='real', rng=rng)

# actually run DMRG; `psi` is updated in-place
numsweeps = 4
en_sweeps = ptn.dmrg_twosite(hamiltonian, psi, numsweeps, tol_split=1e-8)

# value after last iteration
e0 = en_sweeps[-1]
e0

C:\Code\pytenet\pytenet\krylov.py:45: RuntimeWarning: beta[0] ~= 0 encountered during Lanczos iteration.
  warnings.warn(

[6]:

np.float64(-18.48435890403327)

[7]:

psi.bond_dims

[7]:

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

3.3. Evaluate results and compare with exact diagonalization¶

[8]:

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

[8]:

(4096,)

[9]:

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

[9]:

np.float64(0.9999999999999994)

[10]:

# energy after each DMRG sweep
en_sweeps

[10]:

array([-18.4843589, -18.4843589, -18.4843589, -18.4843589])

[11]:

# 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

[11]:

(4096, 4096)

[12]:

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

[12]:

np.float64(1.4210854715202004e-14)

[13]:

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

[13]:

np.float64(-18.484358962762272)

[14]:

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

[14]:

array([5.87290430e-08, 5.87289968e-08, 5.87289932e-08, 5.87290003e-08])