Discrepancy in Calculating Correlation Function Using iTensor TDVP


I have been utilizing the iTensor C++ version for the past three months, working on solving a t-J model with next-nearest-neighbor interaction. Specifically, the Hamiltonian I am dealing with is given by:

H = \sum_{i} J_1 S_i S_{i+1} + \sum_{i} J_2 S_i S_{i+2} + \sum_{i} (t_1 c^\dagger_i c_{i+1} + t_2 c^\dagger_i c_{i+2}) + \text{h.c.}

Now, my objective is to calculate the correlation function

\langle G | O_i(t) O_j(0) | G \rangle

where |G\rangle represents the ground state of the t-J Hamiltonian obtained using DMRG in iTensor,

O_i and O_j are some spin operators (e.g S^Z_i) represented as MPO.

I have written a sample code using the time-dependent variational principle (TDVP) for the nearest-neighbor interaction to achieve this. However, I have encountered a discrepancy between the results obtained from my TDVP code and those from the TEBD calculation. I have cross-validated the TEBD results with exact diagonalization for verification.

Attached to this is my TDVP code for calculating the correlation function. I would greatly appreciate your assistance in identifying any potential issues or discrepancies in my implementation.

Thank you very much for your attention to this matter, and I look forward to your guidance and insights.

Best regards,

Subhajyoti Pal

My code

char title1[50];
sprintf(title1,"SZ.dat"); //| ==== Name | I; 
std::ofstream file0;      
file0.open (title1);
auto O_1 = toMPO(ampo1); //op(sites,"Sz",ncent1);//--automatically returns tensor on site
auto O_2 = toMPO(ampo2);// op(sites,"Sz",ncent2);//--similarly ncent2 = (N/2)+1
auto psi1=psi //psi is calculated from DMRG (Ground state)
auto psi2=psi
state = InitState(sites);
for(int j = 1; j <= N; j++)
    if(j%2==0)  state.set(j,"Up");
    else if(j%2!=0)  state.set(j,"Dn");
auto x0 = MPS(state);
auto y2 = applyMPO(O_2,psi1,x0,{"Method=","Fit","MaxDim=",500,"Cutoff=",1E-12,"Nsweep=",8}); // O_jj|g⟩= |jj⟩
auto y1 = applyMPO(O_1,psi2,x0,{"Method=","Fit","MaxDim=",500,"Cutoff=",1E-12,"Nsweep=",8}); // O_ii|g⟩= |ii⟩


//------------------------- Print the final energies reported by DMRG -------------------------------------//
//printfln("\nGround State Energy = %.16f\t%d\t%d",energy,ncent1,order);
printfln("\nGround State Energy error= %.16f",energy-inner(psi,H,psi));
file0<<t*n<<"   "<< real(innerC(y2,y1)*exp_prefac)<<"   "<<imag(innerC(y2,y1)*exp_prefac)<<endl; // <g|O_ii(0)O_jj(0)|g⟩

// start TDVP, either one site or two site algorithm can be used by adjusting the "NumCenter" argument

double energy0 = real(innerC(psi2,H,psi1));
printfln("Initial energy = %.5f", energy);
printfln("Using overlap = %.10f", real(innerC(y2,y1)));

auto sweeps = Sweeps(1);
sweeps.maxdim() =2000;
sweeps.cutoff() = 1E-12;
sweeps.niter() = 20;

for(int n = 1; n <= nsw; n+=1)

        if(n < 5)
            //Global subspace expansion
            std::vector<Real> epsilonK = {1E-12, 1E-12};

            //  epsilonK = {1E-12, 1E-12};
            // addBasis(y2,H,epsilonK,{"Cutoff",1E-12,
            //                         "Method","Fit",
            //                         "KrylovOrd",3,
            //                         "DoNormalize",true,
            //                         "Quiet",true});
            } //DensityMatrix
       // TDVP sweep
        energy = tdvp(y1,H,-1_i*t,sweeps,{"Truncate",true,

        complex<double> prefac(0,t*n*energy0);
        complex<double> exp_prefac = exp(prefac);
        cout<< exp_prefac*innerC(y2,y1)<< "   "<<innerC(y2,y1)<<"   "<<exp_prefac<<endl;
        #Writting the data of <g|O_ii(t*n)O_jj(0)|g⟩
        file0<<t*n<<"   "<< real(innerC(y2,y1)*exp_prefac)<<"   "<<imag(innerC(y2,y1)*exp_prefac)<<endl;

Hi Subhajyoti,
It’s difficult for me to be sure what is going wrong just from looking at code.

Here are some things I would suggest you may try to reveal why you are getting different results with the two methods:

  • please try running both on a very small system and comparing each one to a third calculation just obtaining the exact answer by exponentiating H as a matrix
  • try comparing a simpler property like \langle S^z_j \rangle to make sure that matches before comparing correlation functions
  • try comparing equal-time correlation functions before comparing time-dependent correlation functions
  • you might also be able to consider the J_1=J_2=0 case which is non-interacting fermions to test your code, though there the spin correlations are probably zero unless you add a chemical potential or an initial state that flips a spin etc.

I hope those suggestions can help you.

Hi miles
Thank you for your suggestion. I will follow all the steps and get back to you.