Severe Discrepancy between MPO and Trotter Time-Evolution?

Question

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

Answer 1

This is still an open issue.....

Answer 2

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

Answer 3

In addition to Matt's answer, which I also agree looks very likely to be the issue, I'd like to note that we have been finding that the toExpH approach to performing time evolution has not always been the most reliable unfortunately (though it depends on details of the Hamiltonian and time step size etc.). So in general you might not be able to get as good of results with that method as, say, the Trotter method for cases where both can be used.

In the near future, we are planning to offer alternative time evolution methods which will supersede the toExpH approach.

Miles