# Discrepancy in Calculating Correlation Function Using iTensor TDVP

**URL:** <https://itensor.discourse.group/t/discrepancy-in-calculating-correlation-function-using-itensor-tdvp/1705>\
**Category:** ITensor C++ Questions\
**Tags:** cpp, tdvp\
**Created:** [May 10, 2024, 3:14pm UTC](https://itensor.discourse.group/t/discrepancy-in-calculating-correlation-function-using-itensor-tdvp/1705 "2024-05-10T15:14:43Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![subhajyoti.pal](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@subhajyoti.pal](https://itensor.discourse.group/u/subhajyoti.pal)\
**Post date:** [May 10, 2024, 3:14pm UTC](https://itensor.discourse.group/t/discrepancy-in-calculating-correlation-function-using-itensor-tdvp/1705/1 "2024-05-10T15:14:43Z")

</div>

Hi,

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⟩
y2.noPrime();
y1.noPrime();

psi1=y1; 
psi2=y2;

//------------------------- 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
println("----------------------------------------GSE-TDVP---------------------------------------");

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.noise()=1E-1;
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};
            addBasis(y1,H,epsilonK,{"Cutoff",1E-12,
                                    "Method","Fit",
                                    "KrylovOrd",3,
                                    "DoNormalize",true,
                                    "Quiet",true});

            // 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,
                                        "DoNormalize",false,
                                        "Quiet",true,
                                        "NumCenter",1,
                                        "ErrGoal",1E-12});

        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;
 }
file0.close();

```

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [May 19, 2024, 9:43pm UTC](https://itensor.discourse.group/t/discrepancy-in-calculating-correlation-function-using-itensor-tdvp/1705/2 "2024-05-19T21:43:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![subhajyoti.pal](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@subhajyoti.pal](https://itensor.discourse.group/u/subhajyoti.pal)\
**Post date:** [June 9, 2024, 6:45am UTC](https://itensor.discourse.group/t/discrepancy-in-calculating-correlation-function-using-itensor-tdvp/1705/3 "2024-06-09T06:45:18Z")

</div>

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