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)
//and add them to gates
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