# Error related to errorMPOProd, and time-evolution (TEBD)

Hello Everyone,

We are facing two issues.

1. We were checking the result of the function errorMPOProd for a two-site operator (A) on a ground state (|G\rangle) obtained using DMRG, where the Hamiltonian (H) does not conserve any quantity (“ConserveQNs=”,false).
We expected the result (\big|\big| |y\rangle−A|G\rangle \big|\big| ^2) to be close to zero.
But we get a large error on sites close to the center of the chain and zero close to the ends of the chain.
Why is this so? Could this have any relation to the MPS being in a mixed canonical state?
How can one resolve this error?

2. We identified this error from another error.
We were checking the expectation value of the Hamiltonian (i) in the ground state (\langle G(0) |H|G(0) \rangle) and (ii) in the time-evolved ground state (\langle G(t) |H|G(t) \rangle) . Both should be the same. But they were not. Instead, on evolving in time using TEBD, we got an oscillation of the time-evolved expectation value, which suggested that there is either an error in the time-evolution or in the starting state itself.

We would appreciate any help or advice.

Subhajyoti

Here I am attaching our code too

H= \sum_{i\in even} JxS^x_iS^x_{i+1}+ \sum_{i\in odd} JyS^y_iS^y_{i+1}

Example operator A=S^+_i

#ifndef MASTER_HPP_INCLUDED
#define MASTER_HPP_INCLUDED
#include "itensor/all.h"
#include <iomanip> // std::setprecision
#include "values.hpp"
#include "input/input.hpp"
#include "input/pob.hpp"
//#include <omp.h>
#include <complex>
#include <vector>
using namespace std;
using namespace Eigen;
using namespace itensor;

void input::Master_print(int ii,int jj,int order)
{
cout<< "======================================================" <<endl;
cout<< "length of the interaction: "<< longrange <<endl;
cout<< "======================================================" <<endl; auto sites    =SpinHalf(N,{"ConserveQNs=",false});
auto ampo = AutoMPO(sites);
//-------------------------------------------- Hamiltonian formation --------------------------------//
//-------------------------------------------- Hamiltonian formation --------------------------------//
for(int i = 1; i <= N; i++)
{
for(int j = i+1; j <= N; j++)
{

if(abs(JJx(i-1,j-1))>0.0 && i%2!=0)
{
//cout<< i <<"    "<<j<<endl;
//---------------------------------------------------------------//
ampo += JJx(i-1,j-1)*0.25,"S+",i,"S-",j; // Heisenberg exchange term
ampo += JJx(i-1,j-1)*0.25,"S-",i,"S+",j; // Heisenberg exchange term

ampo += JJx(i-1,j-1)*0.25,"S+",i,"S+",j; // Heisenberg exchange term
ampo += JJx(i-1,j-1)*0.25,"S-",i,"S-",j; // Heisenberg exchange term

//---------------------------------------------------------------//
//---------------------------------------------------------------//
}
if(abs(JJy(i-1,j-1))>0.0 && i%2==0)
{
//cout<< i <<"    "<<j<<endl;
//---------------------------------------------------------------//
ampo += JJy(i-1,j-1)*0.25,"S+",i,"S-",j; // Heisenberg exchange term
ampo += JJy(i-1,j-1)*0.25,"S-",i,"S+",j; // Heisenberg exchange term

ampo += -JJy(i-1,j-1)*0.25,"S+",i,"S+",j; // Heisenberg exchange term
ampo += -JJy(i-1,j-1)*0.25,"S-",i,"S-",j; // Heisenberg exchange term

//---------------------------------------------------------------//
//---------------------------------------------------------------//
}
}
}
//--------------------------------------------------//
//--------------------------------------------------//
auto H0 = toMPO(ampo);
//-------------------------------------------- Hamiltonian formation --------------------------------//
//-------------------------------------------- Hamiltonian formation --------------------------------//
//-------------------------------------------- Initial MPS --------------------------------//
//-------------------------------------------- Initial MPS --------------------------------//
auto state = InitState(sites);
auto psi0 =randomMPS(sites);// MPS(state);
//-------------------------------------------- Initial MPS --------------------------------//
//-------------------------------------------- Initial MPS --------------------------------//

//-------------------------------------------- GS Calculation --------------------------------//
//-------------------------------------------- GS Calculation --------------------------------//
auto sweeps = Sweeps(20);
sweeps.noise() = 1E-2,1E-6,1E-8,1E-10,1E-12;
sweeps.maxdim() = 50,100,200,500,2000;
sweeps.cutoff() = 1E-10;
auto [energy,psi] = dmrg(H0,psi0,sweeps,{"Quiet=",true});
//-------------------------------------------- GS Calculation --------------------------------//
//-------------------------------------------- GS Calculation --------------------------------//

//------------------------- Print the final energies reported by DMRG -------------------------------------//
printfln("\nGround State Energy = %.16f",energy);
printfln("\nGround State Energy error= %.16f",energy-inner(psi,H0,psi));
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//============ Time Evolution =====================//
//------------------------------------------------------------------------------------------------------------------------------------------

//---------------------------------------------------------------------//
//--------------- creat tievolution file -------------------//
//---------------------------------------------------------------------//
corrsSzfile << "#ttotal <Sz(SL,ttotal).Sz(SR,0)> t(ms) SR   SL"<<endl;
//---------------------------------------------------------------------//
corrsSzfile << std::fixed << std::setprecision(16);
//---------------------------------------------------------------------//
//---------------------------------------------------------------------//

//-------------------------------------------- (Step -2)----------------------------------------------------------------//
//======================================================================================================================//
//==================== Defining MPO ====================//
//======================================================================================================================//
int ncent1 = ii; //------- Central Site
int ncent2 = jj; //------- Next site

//-------------------------------------------------------------------------------------------------------------//
//--------------------- Defining O_1 ---------------------//
//-------------------------------------------------------------------------------------------------------------//
auto ampo1 = AutoMPO(sites);
for(int i = ii-longrange; i <= ii+longrange; i++)
{
if(order==0)
{
if(i==ii-longrange)
ampo1 += 0.5,"S+",ii;
}
else if (order==1)
{
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
if(i>=1 && abs(JJx(ii-1,i-1)+JJy(ii-1,i-1))>0.0 && i<=N)
{
//---------------------------------------------------------------//
ampo1 += (JJx(ii-1,i-1)+JJy(ii-1,i-1))*0.25,"S+",ii,"S+",i;
//---------------------------------------------------------------//
}
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
}
else if (order==2)
{
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
for(int kpp = ii-longrange; kpp <= ii+longrange; kpp++)
{
if (kpp>=1 && kpp<=N)
{
if(i>=1 && abs(JJ(ii-1,i-1)*JJ(ii-1,kpp-1))>0.0 && i<=N)
{
//---------------------------------------------------------------//
ampo1 += JJ(ii-1,i-1)*JJ(ii-1,kpp-1),"Sz",i,"Sz",kpp;
//---------------------------------------------------------------//
}

}
}
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
}

}
//-------------------------------------------------------------------------------------------------------------//
//--------------------- Defining O_2 ---------------------//
//-------------------------------------------------------------------------------------------------------------//
auto ampo2 = AutoMPO(sites);
for(int i = jj-longrange; i <= jj+longrange; i++)
{
if(order==0)
{
if(i==jj-longrange)
ampo2 += 0.5,"S+",jj;
}
else if (order==1)
{
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
if(i>=1 && abs(JJx(jj-1,i-1)+JJy(jj-1,i-1))>0.0 && i<=N)
{
//---------------------------------------------------------------//
ampo2 += (JJx(jj-1,i-1)+JJy(jj-1,i-1))*0.25,"S+",jj,"S+",i;
//---------------------------------------------------------------//
}
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
}
else if (order==2)
{
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
for(int kpp = jj-longrange; kpp <= jj+longrange; kpp++)
{
if (kpp>=1 && kpp<=N)
{
if(i>=1 && abs(JJ(jj-1,i-1)*JJ(jj-1,kpp-1))>0.0 && i<=N)
{
//---------------------------------------------------------------//
ampo2 += JJ(jj-1,i-1)*JJ(jj-1,kpp-1),"Sz",i,"Sz",kpp;
//---------------------------------------------------------------//
}

}
}
//----------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------//
}

}
//-------------------------------------------------------------------------------------------------------------//
//-------------------------------------------------------------------------------------------------------------//
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 psi000 = psi;
//-------------------------- Guess States ---------------------//
auto x0 = randomMPS(sites);//MPS(state);
// //---------------------------------------- ---------------------//
// //---------------------------------------------------------------------------------------------------------------//
auto y2 = applyMPO(O_2,psi,x0,{"MaxDim=",1000,"Cutoff=",1E-12,"Nsweep=",20}); // O_jj|g⟩= |jj⟩,"Method=","Fit",
auto y1 = applyMPO(O_1,psi,x0,{"MaxDim=",1000,"Cutoff=",1E-12,"Nsweep=",20}); // O_ii|g⟩= |ii⟩,"Method=","Fit",
y2.noPrime();
y1.noPrime();
//---------------------------------------------------------------------------------------------------------------//

//---------------------------------------------------------------------------------------------------------------//
printfln("\nError= %.16f",errorMPOProd(y2,O_2,psi)); //should be close to zero |||jj⟩−O|g⟩||^2
printfln("Error= %.16f\n",errorMPOProd(y1,O_1,psi000)); //should be close to zero |||ii⟩−O|g⟩||^2
//---------------------------------------------------------------------------------------------------------------//

//-------------------------------------------- End of (Step -2)---------------------------------------------------------//

//-------------------------------------------- (Step -3)----------------------------------------------------------------//
Real tstep0 = 0.02;//---time step (smaller is generally more accurate)
Real tstep = tstep0;//---time step (smaller is generally more accurate)
Real ttotal_ = 10.0;//---total time to evolve
int t_iters = static_cast<int>(ttotal_/tstep);
std::cout << "t_iters = " << t_iters << '\n';
Real cutoff = 1E-12; //--truncation error cutoff when restoring MPS form
//-------------------------------------------- End of (Step -3) ---------------------------------------------------------//

//-------------------------------------------- (Step -4)----------------------------------------------------------------//
auto gates = vector<BondGate>(); //-- Create a std::vector (dynamically sizeable array) //-- to hold the Trotter gates
//-----------------Create the gates exp(-i*tstep/2*hterm)----------------------------------
//---------------------------- and add them to gates --------------------------------------
for(int b = 1; b <= N-1; ++b)
{
if(abs(JJx(b-1,b))>0.0 && b%2!=0)
{
//---------------------------------------------------------------//
auto hterm = JJx(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S-",b+1);
hterm += JJx(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S+",b+1);

hterm += JJx(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S+",b+1);
hterm += JJx(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S-",b+1);

//---------------------------------------------------------------//
//---------------------------------------------------------------//
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}
if(abs(JJy(b-1,b))>0.0 && b%2==0)
{
//---------------------------------------------------------------//
auto hterm = JJy(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S-",b+1);
hterm += JJy(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S+",b+1);

hterm += JJy(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S+",b+1);
hterm += JJy(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S-",b+1);

//---------------------------------------------------------------//
//---------------------------------------------------------------//
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}

}
//---------------- Create the gates exp(-i*tstep/2*hterm) in reverse -------------------------
//---------------- order (to get a second order Trotter breakup which ------------------------
//------------------ does a time step of "tstep") and add them to gates ----------------------
for(int b = N-1; b >= 1; --b)
{
if(abs(JJx(b-1,b))>0.0 && b%2!=0)
{
//---------------------------------------------------------------//
auto hterm = JJx(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S-",b+1);
hterm += JJx(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S+",b+1);

hterm += JJx(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S+",b+1);
hterm += JJx(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S-",b+1);

//---------------------------------------------------------------//
//---------------------------------------------------------------//
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}
if(abs(JJy(b-1,b))>0.0 && b%2==0)
{
//---------------------------------------------------------------//
auto hterm = JJy(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S-",b+1);
hterm += JJy(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S+",b+1);

hterm += JJy(b-1,b)*0.25*op(sites,"S+",b)*op(sites,"S+",b+1);
hterm += JJy(b-1,b)*0.25*op(sites,"S-",b)*op(sites,"S-",b+1);

//---------------------------------------------------------------//
//---------------------------------------------------------------//
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}
}
//-------------------------------------------- End of (Step -4)----------------------------------------------------------------//

//-------------------------------------------- (Step -5)----------------------------------------------------------------//
//-------------------------------------------- Save initial state ------------------------------------------------------
//-------------------------------------------- auto psi00 = psi --------------------------------------------------------
//-------------------------------------------- Calculate correlations at different ttotal's ----------------------------
//-------------------------------------------- for (int n = 1; n <= t_iters; ++n)---------------------------------------
Real ttotal = 0.0;
//------------------------------------ Calculate <Sz_2(t=0) Sz_1(t=0)> and write to file ------------------------------------
auto t_start0 = std::chrono::high_resolution_clock::now();
//---------------------------------Initialize starting state to be modified -------------------------------------------------

y1 = psi;
y2 = psi;
auto zz = innerC(y1,H0,y2);// <ii|jj> // At time =0
cout << ttotal << " " << zz.real() << " " << zz.imag() << "  " << ncent2 << " " << ncent1 << endl;

auto t_end0 = std::chrono::high_resolution_clock::now(); //
//---------------------------------------------------------------------------------------------------------------------------------//
double elapsed_time_ms0 = std::chrono::duration<double, std::milli>(t_end0-t_start0).count();
corrsSzfile << ttotal << " " << zz.real() << " " << zz.imag() << " " << elapsed_time_ms0 << " " << ncent2 << " " << ncent1 << endl;
std::cout << "ttotal = " << ttotal << ", <O_"<<order<<"("<<ncent1<<").O_"<<order<<"("<<ncent2<<")> = " << zz << endl;
//-------------------------------------------- End of (Step -5)--------------------------------------------------------------------//

//--------------------------------------------Time Evolution ------------------------------------------------------------------//
//--------------------------------------------------- (Step -6)----------------------------------------------------------------//
//--------------------------- Calculate and write for nonzero ttotal's in the following loop --------------------------------
//--------------------------------------------- while (ttotal < ttotal_+tstep) ----------------------------------------------
while (ttotal < ttotal_)
{
auto t_start = std::chrono::high_resolution_clock::now();
//---------------- Real ttotal = n*t_iters---------------//
ttotal += tstep0;
//-------------- Initialize starting states -------------//
auto psi00 = y2;
auto phi0 = y1;
//-------------------------------------------------------//
//-------------------------------------------------------//
//--------------Time-evolve, overwriting psi00 when done--------------//
gateTEvol(gates,tstep0,tstep,y2,{"Cutoff=",cutoff,"Verbose=",true});
auto zz = innerC(y2,H0,y2); // At time = ttotal;
//---------------------------------------------------------------------------------------------------------------------------------------//
auto t_end = std::chrono::high_resolution_clock::now();
double elapsed_time_ms = std::chrono::duration<double, std::milli>(t_end-t_start).count();
corrsSzfile << ttotal << " " << 2.0*zz.real()/2.0 << " " << 2*zz.imag()/2.0 << " " << elapsed_time_ms << " " << ncent2 << " " << ncent1 << '\n';
std::cout << "ttotal = " << ttotal << ", <O_"<<order<<"("<<ncent1<<").O_"<<order<<"("<<ncent2<<")> = " << zz << endl;
//---------------------------------------------------------------------------------------------------------------------------------------//
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
cout<< "======================================================" <<endl;
auto diff=end-start;
cout<<"Time for single point calculation calculation:  "<< chrono::duration <double,milli>(diff).count()<<"   "<<endl;
cout<< "======================================================" <<endl;
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------


}
#endif

We mentioned here step by step The above 2 issues where they are occurring. if we have codded anything wrong then mention that also.
We would appreciate any help or advice.