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();