Discrepancy in correlation functions calculated from ED and DMRG for free fermions

Question

I was wondering if you could help me with the following issue: For the sake of benchmarking my code, I am considering free fermions hopping on a 1D lattice, and set the filling to 2 fermions only. I furthermore set the spins of both electrons to be up. I use DMRG to solve this problem.
On the other hand, I solve it by diagonalizing the Hamiltonian matrix.

Then I look at correlation functions such as $\langle \psi | c^\daggeri cj | \psi \rangle$ in the above two cases, but the off-diagonal elements of such matrix differ in the two cases.

Could you please check and let me know how I can resolve this problem? I provided a minimal code below, sorry it got long. I used your own example of DMRG for the extended Hubbard model.

Thank you

using ITensors
using LinearAlgebra


N = 10
Npart = 2
t1 = 1.0

sites = siteinds("Electron", N;
               conserve_qns = true)

ampo = AutoMPO()
for b=1:N-1
ampo += -t1, "Cdagup", b,    "Cup", b+1
ampo += -t1, "Cdagup", b+1,  "Cup", b
ampo += -t1, "Cdagdn", b,    "Cdn", b+1
ampo += -t1, "Cdagdn", b+1,  "Cdn", b
end
H = MPO(ampo,sites)

sweeps = Sweeps(6)
maxdim!(sweeps,50,100,200,400,800,800)
cutoff!(sweeps,1E-12)

state = ["Emp" for n=1:N]
for i=1:Npart
    state[i] = "Up"
end

psi0 = randomMPS(sites, state, min(10, floor(Int,N/2) ))

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy, psi = dmrg(H, psi0, sweeps)



println("\nGround State Energy = $energy")


#Below, correlation functions using the DMRG solution are calculated:


cor_mat_mps = zeros(N,N);

for i = 1:N
    orthogonalize!(psi,i)

    s_ind_i = siteind(psi,i);

    for j =1:N
        if j < i
            continue
        end

        s_ind_j = siteind(psi,j);

        if i == j
            cor_mat_mps[i,i] = scalar( 
                psi[i] * prime( dag(psi[i]) , s_ind_i ) * op(s_ind_i,"Nup")
                );
        else

            prod_aux = psi[i];

            for n = i+1:j
                prod_aux *= psi[n];
            end

            cor_mat_mps[i,j]=   scalar(
                prod_aux * prime( dag(prod_aux) , s_ind_i, s_ind_j) *
                   op(s_ind_i,"Cdagup")*op(s_ind_j,"Cup")
                ); 

            cor_mat_mps[j,i]= scalar(
                prod_aux * prime( dag(prod_aux) , s_ind_i, s_ind_j) *
                   op(s_ind_j,"Cdagup")*op(s_ind_i,"Cup")
                ) ; 

        end
    end
end

# The same Hamiltonian is diagonalized below and correlation functions are calculated.

h_nonint = zeros(N,N);

for n = 1:N
    for m = 1:N
#         h_nonint[n,m] = 0.01*rand();
        if abs(m-n) == 1
            h_nonint[n,m] = -1. ;
        end

    end
end


diagonalization_energy = sum( eigen(Hermitian(h_nonint)).values[1:2] );

println("\nGround State Energy from diagonalization = $diagonalization_energy")


vecs = eigen(Hermitian(h_nonint)).vectors;

cor_direct = zeros( N,N );


for m = 1:N
    for n = 1:N
        cor_direct[m,n] = conj(vecs[m,1]) * vecs[n,1];
        cor_direct[m,n] += conj(vecs[m,2]) * vecs[n,2];
    end
end


println()

# below, the difference between correaltion functions calcualted in the two ways is shown, first only the diagonal terms and then all the correlation matrix.
println(  norm(diag(cor_direct-cor_mat_mps))  )
println(  norm(cor_direct-cor_mat_mps)  )

Answer 1

Hi, so I see in your code that you are not including any Jordan-Wigner “F” string operators in your MPS correlation function calculations. This is necessary when computing pairs of operators such as C^dagger and C which have an odd fermion parity. (As a side comment, we are nearly done with a new feature that will remove this requirement, using a technique that makes tensor indices anticommutative, but it’s technical and still in the debugging phase.)

For some resources on how to do this and why, here is a sample code for the C++ version showing how to compute fermion correlators properly for the case of spinless fermions:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/spinless_correlator_mps

For a discussion about Jordan-Wigner string and mappings to bosonic operators (all ITensors are bosonic in that sense), see this tutorial page: http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/fermions

Finally, the approach you are taking of making prod_aux where you combine many MPS tensors together looks like it will have an exponential cost, if I am understanding your code correctly there. This page has an explanation of the approach that is efficient for calculating MPS correlators: http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations

We can discuss some more if you have questions -

Best regards,
Miles

Comments

Thank you for your elaborate response. I see the problem now; I was assuming that the Jordan-Wigner string is taken into account automatically when using C operators, since things worked properly with AutoMPO. But it seems it's not the case. Looking forward to working with this new feature that you mentioned.
Also, thanks for providing the more efficient approach to calculate the correlation functions, I was aware of this other approach but used the one above for clarity.
Thank you!

---

I was wondering if you could help me with how I can adopt the code sample you provided above to the spinful case?
For the spinful case, when I have both spin up and down particles, the code does not give correct correlation function. Please see a sample code below.
Thanks

===============================
    using ITensors

    N = 10
    Npart = 2
    t1 = 1.0

    sites = siteinds("Electron", N;
                   conserve_qns = true)

    ampo = AutoMPO()
    for b=1:N-1
    ampo += -t1, "Cdagup", b,    "Cup", b+1
    ampo += -t1, "Cdagup", b+1,  "Cup", b
    ampo += -t1, "Cdagdn", b,    "Cdn", b+1
    ampo += -t1, "Cdagdn", b+1,  "Cdn", b
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(6)
    maxdim!(sweeps,50,100,200,400,800,800)
    cutoff!(sweeps,1E-12)

    state = ["Emp" for n=1:N]
    for i=1: floor(Int,Npart/2)
        state[i] = "Up"
        state[floor(Int,Npart/2)+i] = "Dn"
    end


    psi0 = randomMPS(sites, state, min(10, floor(Int,N/2) ))

    # Check total number of particles:
    @show flux(psi0)

    # Start DMRG calculation:
    energy, psi = dmrg(H, psi0, sweeps)



    println("\nGround State Energy = $energy")


    #Below, correlation functions using the DMRG solution are calculated:


    cor_mat_mps = zeros(N,N);

    for i = 1:N
        orthogonalize!(psi,i)

        s_ind_i = siteind(psi,i);
        l_ind_i = linkind(psi,i);

        for j =1:N
            if j < i
                continue
            end


            if i == j
                cor_mat_mps[i,i] = scalar(
                    psi[i] * prime( dag(psi[i]) , s_ind_i ) * op(s_ind_i,"Nup")
                    );
            else

                s_ind_j = siteind(psi,j);
                l_ind_j_left = linkind(psi,j-1);


                prod_aux = psi[i] * prime(dag(psi[i]), s_ind_i, l_ind_i);



                for n = i+1:j - 1
                    prod_aux *= psi[n] ;
                    prod_aux *= op(sites,"F",n) ;
                    prod_aux *= prime( dag( psi[n] ) ) ;
                end



                cor_mat_mps[i,j]=   scalar(
                    prod_aux * psi[j] * prime( dag( psi[j] ), s_ind_j, l_ind_j_left) *
                       op(s_ind_i,"Cdagup")*op(s_ind_j,"Cup")
                    );


                cor_mat_mps[j,i]= cor_mat_mps[i,j];

            end
        end
    end