# Severe Discrepancy between MPO and Trotter Time-Evolution?

edited

I am comparing the two methods for time-evolution.

First, I write the Hamiltonian in MPO form using autoMPO, and then time evolve using toExpH and applyMPO.

Secondly, I implement a TEBD type algorithm by creating a sequence of two-site gates using BondGate and time evolving with gateTEvol.

However, both give very different results for the dynamics for a Fermi-Hubbard model. I have tested this for a spin-1/2 chain and this works fine, I am wondering if the trotter gate method does not properly take into account the fermionic statistics?

Here is sample code:

#include "itensor/all.h"

using namespace itensor;

int main(int argc, char* argv[])
{
int BD = 1024;

int t1 = 1;

int N = 20;
auto sites = Electron(N);
auto state = InitState(sites);
int p=1;
for(int i = N; i >= 1; --i)
{
if(p == 1)
{
state.set(i,"Up");
p = 2;
}
else
{
state.set(i,"Dn");
p = 1;
}
}
auto psi = MPS(state);

//
// Create the Hamiltonian using AutoMPO
//
auto ampo = AutoMPO(sites);
for(int b = 1; b < N; ++b)
{
ampo += -t1,"Cdagup",b,"Cup",b+1;
ampo += -t1,"Cup",b,"Cdagup",b+1;
ampo += -t1,"Cdagdn",b,"Cdn",b+1;
ampo += -t1,"Cdn",b,"Cdagdn",b+1;
}
auto H = toMPO(ampo);
auto tau = 0.01;
auto expH = toExpH(ampo,tau*Cplx_i);

auto ttotal = 1.0;
auto nt = int(ttotal/tau+(1e-9*(ttotal/tau)));

std::ofstream ofGup("DATA/MPO_TE_Test.txt");

for(int n = 1; n <= nt; ++n)
{
cpu_time sw_time;

psi = applyMPO(expH,psi,{"Method=","DensityMatrix","MaxDim=",BD,"Cutoff=",1E-15});
psi.noPrime().normalize();

double updo;

for(int j = 1; j <= N; ++j)
{
psi.position(j);
updo = (eltC(dag(prime(psi(j),"Site"))*op(sites,"Nup",j)*psi(j))).real();
printfln(ofGup,"%d ",updo);
}

printfln("MPO: Time %d",n*tau);
printfln("Maximum MPS bond dimension is %d",maxLinkDim(psi));

auto sm = sw_time.sincemark();
printfln("CPU time = %s ",showtime(sm.time));
}

auto sites2 = Electron(N);
auto state2 = InitState(sites2);
p=1;
for(int i = N; i >= 1; --i)
{
if(p == 1)
{
state2.set(i,"Up");
p = 2;
}
else
{
state2.set(i,"Dn");
p = 1;
}
}
auto psi_TEBD = MPS(state2);

Print(totalQN(psi_TEBD));

auto gates = vector<BondGate>();

//Create the gates exp(-i*tstep/2*hterm)
for(int b = 1; b <= N-1; ++b)
{
auto hterm = -t1*op(sites2,"Cdagup",b)*op(sites2,"Cup",b+1);
hterm += -t1*op(sites2,"Cup",b)*op(sites2,"Cdagup",b+1);
hterm += -t1*op(sites2,"Cdagdn",b)*op(sites2,"Cdn",b+1);
hterm += -t1*op(sites2,"Cdn",b)*op(sites2,"Cdagdn",b+1);

auto g = BondGate(sites2,b,b+1,BondGate::tReal,tau/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)
{
auto hterm = -t1*op(sites2,"Cdagup",b)*op(sites2,"Cup",b+1);
hterm += -t1*op(sites2,"Cup",b)*op(sites2,"Cdagup",b+1);
hterm += -t1*op(sites2,"Cdagdn",b)*op(sites2,"Cdn",b+1);
hterm += -t1*op(sites2,"Cdn",b)*op(sites2,"Cdagdn",b+1);

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

Real cutoff = 1E-15; //truncation error cutoff when restoring MPS form
std::ofstream ofGup2("DATA/TEBD_TE_Test.txt");

for(int n = 1; n <= nt; ++n)
{
cpu_time sw_time;

//Time evolve, overwriting psi when done
gateTEvol(gates,tau,tau,psi_TEBD,{"Cutoff=",cutoff,"Verbose=",true,"MaxDim=",BD});
psi_TEBD.noPrime().normalize();

double updo;

for(int j = 1; j <= N; ++j)
{
psi_TEBD.position(j);
updo = (eltC(dag(prime(psi_TEBD(j),"Site"))*op(sites2,"Nup",j)*psi_TEBD(j))).real();
printfln(ofGup2,"%d ",updo);
}

printfln("Time %d",n*tau);
printfln("Maximum MPS bond dimension is %d",maxLinkDim(psi_TEBD));

auto sm = sw_time.sincemark();
printfln("CPU time = %s ",showtime(sm.time));
}

return 0;
}


What have I done wrong??
Note that this is not an issue with the time-step or bond dimension.

Thanks,

Stuart

+1 vote
edited

This is still an open issue.....

Hi Stuart,

I believe the issue is that you should define your Hamiltonian as follows:

for(int b = 1; b < N; ++b)
{
ampo += -t1,"Cdagup",b,"Cup",b+1;
ampo += -t1,"Cdagup",b+1,"Cup",b;
ampo += -t1,"Cdagdn",b,"Cdn",b+1;
ampo += -t1,"Cdagdn",b+1,"Cdn",b;
}


to correctly account for the Fermionic sign (take a look at: https://itensor.org/docs.cgi?page=tutorials/fermions for more information).

Let me know if that helps.

-Matt