Inconsistency between result of dmrgpy and ITensor.jl

[HamidArianZad]

Dear Jose,

I would like to report an inconsistency between results obtained for the ground-state energy of a four-site cluster Hubbard chain using dmrgpy and Itensor.jl. The Hubbard model includes two different hopping term t=1 (intra-chain), and tp=0.06t (inter-chain).

The code that I used is as follows:

Nfst = 6 # number of four-site cells

t = 1             # intra-chain hopping term: t=1 -> energy in units of t
tp = 0.06     # inter-chain hopping term
U = 1             # Hubbard term
mu = 0.5      # Chemical potential

B = 0.1      # Zeeman term 

Ne = 4 * Nfst   # total number of sites

fc = fermionchain.Spinful_Fermionic_Chain(Ne)  # create the chain

h = 0    # Hamiltonian 

for unit in range(Nfst):   # intra-chain hopping terms
    h -= t * fc.Cdagup[4 * unit + 2] * fc.Cup[4 * unit]
    h -= t * fc.Cdagdn[4 * unit + 2] * fc.Cdn[4 * unit]
    h -= t * fc.Cdagup[4 * unit + 2] * fc.Cup[4 * unit + 1]
    h -= t * fc.Cdagdn[4 * unit + 2] * fc.Cdn[4 * unit + 1]
    h -= t * fc.Cdagup[4 * unit + 2] * fc.Cup[4 * unit + 3]
    h -= t * fc.Cdagdn[4 * unit + 2] * fc.Cdn[4 * unit + 3]

for unit in range(Nfst-1):    # inter-chain hopping term
    h -= tp * fc.Cdagup[4 * unit + 3] * fc.Cup[4 * unit + 4]
    h -= tp * fc.Cdagdn[4 * unit + 3] * fc.Cdn[4 * unit + 4]
    h += mu * (fc.Nup[4 * unit + 3] + fc.Ndn[4 * unit + 3])   
    h -= mu * (fc.Nup[4 * unit + 4] + fc.Ndn[4 * unit + 4])

h = h + h.get_dagger()

for i in range(4*Nfst):
    h -= B * fc.Sz[i]                    # Field term
    h += U*(fc.Nup[i]-0.5)*(fc.Ndn[i]-0.5)     # Hubbard term

fc.set_hamiltonian(h)                       # set the Hamiltonian to the object
wf = fc.get_gs()                                  # compute ground state
Mz = [fc.vev(Szi) for Szi in fc.Sz]      # magnetization per site

Sz_cluster = np.zeros(Nfst)              # total magnetization per unit block (four-site clusters)

nbt = [-2, -1, 0, +1]                           # sites per cluster
for unit in np.arange(0, Nfst):
    for i in nbt:
        Sz_cluster[unit] += wf.dot(fc.Sz[4 * unit + 2 + i]*wf).real

After spending much time, I understood that When I consider tp <0.35, I get a large discrepancy for the results of Itensor.jl. To remove all degeneracies I added a finite Zeeman term and also a chemical potential to fix the total quantum number as I wish. In fact, I need to have ground state |S, S_zT> = |1, +1> in half-filling. Everything seems to be fine when tp>0.35, but when I take tp<0.35 I loss the consistency between dmrgpy and Itensor.jl results.

Have you any idea to address this problem?

[joselado]

Dear HamidArianZad,

thank you for the message. I have just checked and the ground state energies for a small system with the different methods are consistent, so I would say that there is in principle no issue.

The difference may stem from some additional degeneracy, or a different convergence with the bond dimension. As a test, perhaps the way to proceed is to see if you can reproduce this difference in a small system that you can also solve exactly with exact diagonalization or analytically.

Best regards, Jose

[HamidArianZad]

Dear Jose,

Thank you very much for your re-check the results in both dmrgpy and Itensor. In this regard, I see some differences between results obtained from QuSpin package and dmrgpy. I am just curious to realize that what is the difference between creation and annihilation operators and normalization in QuSpin package and dmrgpy? For example, I used the below Hamiltonian in both packages to gain the ground-state energy and some low-lying excited states and their energies. But I get completely different results for the ground and excited states and their energies.

Nfst = 2     # Number of unit blocks 
Nsites = 4 * Nfst
basis = spinful_fermion_basis_1d(Nsites, Nf=(Ne_up, Ne_dn))

CdagsCs = []  # s=up,dn
Nupdn = []

 # hopping term
    if t != 0.0:
        nbt = [-2, -1, +1]
        CdagsCs += [[-t, 4 * unit + 2, 4 * unit + 2 + j] for unit in range(Nfst)
                    for j in nbt]
        # hc
        CdagsCs += [[-t, 4 * unit + 2 + j, 4 * unit + 2] for unit in range(Nfst)
                    for j in nbt]

    # tp: hopping between nearest neighbor sites of each adjacent block
    if tp != 0.0:
        CdagsCs += [[-tp, 4 * unit + 3, 4 * unit + 3 + 1] for unit in range(Nfst - 1)]
        # hc
        CdagsCs += [[-tp, 4 * unit + 3 + 1, 4 * unit + 3] for unit in range(Nfst - 1)]

    # U
    Nupdn += [[U, i, i] for i in range(4 * Nfst)]

    static = [
        ['+-|', CdagsCs],  # CdagupCup
        ['|+-', CdagsCs],  # CdagdnCdn
        ['n|n', Nupdn]
    ]
    dynamic = []

    no_checks = dict(check_pcon=False, check_symm=False, check_herm=False)
    H = hamiltonian(static, dynamic, basis=basis, dtype=np.float64, **no_checks)
   En, psin = H.eigsh(k=3, which='SA')

After running this code in QuSpin for one, two and three number of unit blocks where in each unit block we have an imbalance term |N_up-N_dn| = 2 such that we will have the ground state |S, SzT> = |1, +1>, I get following results:

N_up= 3 N_dn= 1
Energies =  [**-2.84821708** -1.17008649 -1.17008649 -0.73205081]  # -2.84821708 is the ground-state eigenenergy
N_up= 5 N_dn= 3
Energies =  [**-5.70230319** -5.69728019 -5.46228373 -5.4585275 ]   # -5.70230319 is the ground-state eigenenergy
Ne_up= 7 Ne_dn= 5
Energies =  [**-8.55877328** -8.55382438 -8.55376966 -8.55136592]  # -8.55877328 is the ground-state eigenenergy

while the same Hamiltonian (below code) in dmrgpy shows different values of the ground and low-lying eigenenergies when I fix the number of up and down particles of each block with the imbalance term |N_up-N_dn| = 2 to have the ground state |S, SzT> = |1, +1>.

Nfst = 1

t = 0.5*1  # t=1 -> energy in units of t
tp = 0.06
U = 1
mu = -0.45
B = -0.1

Ne = 4 * Nfst 

fc = fermionchain.Spinful_Fermionic_Chain(Ne)  # create the chain
h = 0
for unit in range(Nfst):
    if t != 0.0:
        nbt = [-2, -1, +1]
        for j in nbt:
            h -= t * fc.Cdagup[4 * unit + 2] * fc.Cup[4 * unit + 2 + j]
            h -= t * fc.Cdagdn[4 * unit + 2] * fc.Cdn[4 * unit + 2 + j]
            h += mu * (fc.Nup[4 * unit] + fc.Ndn[4 * unit])
            h += mu * (fc.Nup[4 * unit + j] + fc.Ndn[4 * unit + j])

    if tp != 0.0 and unit < Nfst-1:
        h -= tp * fc.Cdagup[4 * unit + 3] * fc.Cup[4 * unit + 3 + 1]
        h -= tp * fc.Cdagdn[4 * unit + 3] * fc.Cdn[4 * unit + 3 + 1]
        h += mu * (fc.Nup[4 * unit + 3] + fc.Ndn[4 * unit + 3])
        h -= mu * (fc.Nup[4 * unit + 4] + fc.Ndn[4 * unit + 4])

for i in range(4 * Nfst):
    h += U * (fc.Nup[i] - 0.5) * (fc.Ndn[i] - 0.5)
    h += B * fc.Sz[i]

h = h + h.get_dagger()

fc.set_hamiltonian(h)
Mz = [fc.vev(Szi) for Szi in fc.Sz]
es1, ws1 = fc.get_excited_states(n=3, mode="DMRG")
print("Exited Energy 0-2 = ", es1)

Please note that I could not fix the number of particles in half-filling mode such that N_up-N_dn = 2 that lead to SzT = +1, except by imposing a chemical potential mu = -0.45 for last site of a block and the first site of its next block with opposite sign, and multiplying hopping term t by the factor of 0.5. With this consideration I could produce the state |1, +1> in dmrgpy. Below are the dmrgpy results for the energy of the ground state |S, SzT> = |1, +1> with different number of unit blocks Nfst:

Nfst = 1 # for one block
GS Energy =  -2.2700864866260377
Particle number =  4.000000000000007
Sz_total =  1.0000000000000009
N_up = 2.9999999999999996
N_dn = 0.9999999999999999
N_total = 4.0

-----------------------------------------------------------------
Nfst = 2 # for two blocks
GS Energy =  -5.125375950332136
Particle number =  8.000000000000007
Sz_total =  0.9999999999819709
N_up = 4.9999999999999996
N_dn = 3.0000000000032516
N_total = 8.0
-----------------------------------------------------------
Nfst = 3 # for three blocks
GS Energy =  -7.995624536193942
Particle number = 12.000000000065112
Sz_total =  0.9999999999999878
N_up = 7.000000000032609
N_dn = 5.000000000032608
N_total = 12.000000000065214

where GS Energy is different from that of obtained from QuSpin. I also get different results for the eigenenergies of excited states.

1. Is there any difference between introduced creation and annahilation operators in QuSpin and dmrgpy?
2. Is there any difference between normalization of the Hamiltonian parameters in QuSpin and dmrgpy?
3. I get different results when I consider:

h = h + h.get_dagger()  

for i in range(4 * Nfst):
    h += U * (fc.Nup[i] - 0.5) * (fc.Ndn[i] - 0.5)
    h += B * fc.Sz[i]

and when

for i in range(4 * Nfst):
    h += U * (fc.Nup[i] - 0.5) * (fc.Ndn[i] - 0.5)
    h += B * fc.Sz[i]

h = h + h.get_dagger()  

How can I realize where is the correct place of the term: h = h + h.get_dagger() ?

If you let me know your valuable comments, it would be very helpful for me and I am very thankful to you for your responses.

Best regards, Hamid

[joselado]

Dear Hamid,

thanks for the follow up. Unfortunately I am not aware of the internal the details of QuSpin, so I would not know if there is any difference between dmrgpy and QuSpin.

Perhaps the simplest approach would be take a model where you know the analytic solution, and see if the codes are consistent with it.

Regarding the statement h = h + h.get_dagger() , this is just a minimal way of including all the hermitian conjugate terms. If you include explicitly all the terms in the Hamiltonian, this statement does not have to be included.

Best regards, Jose