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^\dagger*i c*j | \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) )
```