Dear ITensor Team,

Thanks for your reply. Recently, I want to repeat entanglement entropy calculation as Fig.1 in Phys. Rev. Lett. 110, 115701 by using ITensor v3, which is just simulate 1/4 doped hubbard chain. I followed the tutorial http://itensor.org/docs.cgi?vers=cppv3&page=formulas/entanglement_mps

But i find my results(entropy and energy) are not converge ( we get different results every time )when i increase bond dimension. Is there any bug in my naive code? Please help me. Many thanks!

There is my naive code:

int main()

{

int Nx = 24;

auto N = Nx;

auto sites = Electron(N);

auto t = 1.0;

auto ampo = AutoMPO(sites);

for(int j=1; j<N; ++j)

{

ampo += -t, "Cdagup",j,"Cup",j+1;

ampo += -t, "Cdagdn",j,"Cdn",j+1;

}

for(int j=1; j<=N;++j)

{

ampo +=UU, "Nupdn",j;

}

auto H = toMPO(ampo);

auto state = InitState(sites,"Emp");

// for this part, naively, i want to doped my hubbard chain, is that correct?

for(int i=1; i<=Nx/4; ++i)

{

state.set(i,"Up");

}

for(int i=Nx/4+1; i<=Nx/2; ++i)

{

state.set(i,"Dn");

}

auto sweeps = Sweeps(10);

sweeps.maxdim() = 10, 20, 50,500;

sweeps.noise() = 1E-7, 1E-8, 1E-10, 0;

sweeps.cutoff() = 1E-10;

PrintData(sweeps);

auto psi0 = MPS(state);

auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});

PrintData(Nx);

PrintData(UU);

PrintData(t);

PrintData(totalQN(psi));

PrintData(maxLinkDim(psi));

PrintData(energy);

//Given an MPS called "psi",

//and some particular bond "b" (1 <= b < length(psi))

//across which we want to compute the von Neumann entanglement

auto b = Nx/2;

//"Gauge" the MPS to site b

psi.position(b);

//SVD this wavefunction to get the spectrum

//of density-matrix eigenvalues

auto l = leftLinkIndex(psi,b);

auto s = siteIndex(psi,b);

auto [U,S,V] = svd(psi(b),{l,s});

auto u = commonIndex(U,S);

//Apply von Neumann formula

//to the squares of the singular values

Real SvN = 0.;

for(auto n : range1(dim(u)))

{

auto Sn = elt(S,n,n);

auto p = sqr(Sn);

if(p > 1E-12) SvN += -p*log(p);
}
printfln("Across bond b=%d, SvN = %.10f",b,SvN-2*log(Nx)/3.0);

printfln("Across bond b=%d, SvN = %.10f",b,SvN);

return 0;

}