evolution with complicated Hamiltonian using Trotter gates

Question

Hello everyone!

How one can evolve the folloing hamiltonian using Trotter gates?

$$
H=\sum_{i=1}^{L-1}h_{i,i+1} +\sum_{i=1}^{L}T_i
$$

For instance, for this Hamiltonian
$$
H=\sum_{i=1}^{L-1}S^x _i S^x _{i+1} +\sum_{i=1}^{L}S^z _i
$$

using Gate = BondGate<ITensor>;

auto gates = vector<Gate>();


for(int b = 1; b <= N-1; ++b)
    {
    ITensor Sxb = sites.op("Sx",b); //Sx must be converted to an ITensor prior to usage
    ITensor Sxb1 = sites.op("Sx",b+1);
    auto hterm = Sxb*Sxb1;
    auto g = Gate(sites,b,b+1,Gate::tReal,tstep/2.,hterm);
    gates.push_back(g);
    }

if I write

hterm +=sites.op("Sz",b);

in the for-loop then I obtain the Hamiltonian without the last term @@S^z_{L}@@, but if I write

hterm +=0.5*sites.op("Sz",b)+0.5*sites.op("Sz",b+1);

then there is the discrepancy for the first and the last terms, i.e. @@S^z_{1}@@ and @@S^z_{L}@@.

Comments

Actually the first code,

hterm +=sites.op("Sz",b)

is not working. Error occured:
-------------------------------------------------
ITensorT::operator+=: different number of indices
ITensorT::operator+=: different number of indices
Aborted
----------------------------------------------------
whereas the second code

hterm +=0.5*sites.op("Sz",b)+0.5*sites.op("Sz",b+1);

gives also an error

------------------------
ITensorT::operator+=: different index structure

ITensorT::operator+=: different index structure
Aborted

---

The full code is attached if necessary (without additional term @@\sum S^z\_i@@)

#include "itensor/all.h"

using namespace itensor;
using std::vector;

int main()
{
int N = 10; //number of sites
Real tstep = 0.1; //time step (smaller is generally more accurate)
Real ttotal = 1.0; //total time to evolve
Real cutoff = 1E-5; //truncation error cutoff when restoring MPS form

//Define a site set object "sites" which lets us
//easily obtain Site indices defining our Hilbert space
//and S=1/2 single-site operators
auto sites = SpinHalf(N);
auto state = InitState(sites,"Dn");//initial state
auto psi = MPS(state);

//Define the type "Gate" as a shorthand for BondGate<ITensor>
using Gate = BondGate<ITensor>;

//Create a std::vector (dynamically sizeable array)
//to hold the Trotter gates
auto gates = vector<Gate>();

//Create the gates exp(-i*tstep/2*hterm)
//and add them to gates

for(int b = 1; b <= N-1; ++b)
    {
    ITensor Sxb = sites.op("Sx",b); //Sx must be converted to an ITensor prior to usage
    ITensor Sxb1 = sites.op("Sx",b+1);
    auto hterm = Sxb*Sxb1;
    hterm +=0.5*sites.op("Sz",b)+0.5*sites.op("Sz",b+1);
   
    auto g = Gate(sites,b,b+1,Gate::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)
    {
    ITensor Sxb = sites.op("Sx",b);
    ITensor Sxb1 = sites.op("Sx",b+1);
    auto hterm = Sxb*Sxb1;
    hterm +=0.5*sites.op("Sz",b)+0.5*sites.op("Sz",b+1);
   
    auto g = Gate(sites,b,b+1,Gate::tReal,tstep/2.,hterm);
    gates.push_back(g);
   
    }


//Time evolve, overwriting psi when done
gateTEvol(gates,ttotal,tstep,psi,{"Cutoff=",cutoff,"Verbose=",true});

printfln("Maximum MPS bond dimension after time evolution is %d",maxM(psi));

//define an operator to find expactation value
auto ampo = AutoMPO(sites);
for(int j = 1; j <= N; ++j)
    {
    ampo += "Sz",j;
    }
auto H1 = MPO(ampo);

//Print <psi_t|H1|psi_t>
//(Will be complex so using overlapC which can return complex);
Print(overlapC(psi,H1,psi).real());

return 0;
}

---

I tried the following and it works:
    hterm += sites.op("Sz",b)*sites.op("Id",b+1);// automatically convert to ITensor
    if(b==N-1) hterm += sites.op("Id",b)*sites.op("Sz",b+1);

Answer 1

You can just add a if sentence in your for loop.
Something like:
hterm +=sites.op("Sx",b);
if(b == N-1) hterm +=sites.op("Sx",b+1);

Comments

The snippet
------------------------
for(int b = 1; b <= N-1; ++b)
    {
    ITensor Sxb = sites.op("Sx",b); //Sx must be converted to an ITensor prior to usage
    ITensor Sxb1 = sites.op("Sx",b+1);
    auto hterm = Sxb*Sxb1;
    hterm +=sites.op("Sz",b);
    if(b==N-1)
    {
    hterm += sites.op("Sz",b+1);   
    }   
   
    auto g = Gate(sites,b,b+1,Gate::tReal,tstep/2.,hterm);
    gates.push_back(g);
    }
------------------------
 gives rise to the error

ITensorT::operator+=: different number of indices
ITensorT::operator+=: different number of indices
Aborted

---

Please print out the different ITensors you are adding together and inspect the results of the print out. This should enable you to see why you are getting tensors with different number of indices.

It's important to note that only ITensors with the same index structure (same set of indices) are allowed to be added together.

(At a quick glance, it looks like perhaps you are trying to add an operator which acts on two sites to one which only acts on one site. I think the solution is that you need to tensor-product the Sz operator on site b with an identity operator acting on site b+1 before adding it to hterm.)

---

Thank you very much, Miles!

hterm += sites.op("Sz",b)*sites.op("Id",b+1);

is working.