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

subhajyoti.pal · May 16, 2024, 8:53am

Hello Everyone,

We are facing two issues.

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?
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.
Thanks in advance.

Subhajyoti

subhajyoti.pal · May 17, 2024, 8:07am

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)
{
  auto start =chrono::steady_clock::now();
  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 -------------------// 
//---------------------------------------------------------------------//
std::ofstream corrsSzfile;      
corrsSzfile.open (title1);
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;
    //---------------------------------------------------------------------------------------------------------------------------------------//
}
corrsSzfile.close();
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
cout<< "======================================================" <<endl; 
auto end =chrono::steady_clock::now();
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.
Thanks in advance.

Subhajyoti

Topic		Replies	Views
TEBD for electrons with hopping and onsite energy ITensor Julia Questions julia , mps , fermions	6	365	October 4, 2024
Errors in TEBD on a Mixed-Site Model DMRG and Numerical Methods julia , mps	1	255	July 31, 2023
Unexpected oscillations in long-time TEBD results ITensor Julia Questions julia , mps	7	112	June 21, 2024
TEBD on 2D fermionic system ITensor Julia Questions	3	123	September 12, 2024
std::bad_alloc when performing dmrg ITensor C++ Questions	5	261	August 3, 2022

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

Related topics