# different results for Hamiltonian written in different form

+1 vote
edited

Here I have a Bose-Hubbard model (with self-defined siteset file) and the density-density interaction can be written in two identical forms as follows

ampo += U, "Nb", i, "Nb", i;
ampo += -U, "Nb", i;

and

this two forms of "ampo" should be the same because @@N_b=A^{\dagger}A@@. However, the final results of dmrg is quite different for this two forms. In fact, the first form gives the correct answer while the second one does not and it just converge very slowly.

Here I presents the siteset file for Bose Hubbard model

#ifndef __ITENSOR_BH_H
#define __ITENSOR_BH_H
#include "itensor/mps/siteset.h"
#include <cmath>

namespace itensor
{
class BoseHubbard;
using BH=BasicSiteSet<BoseHubbard>;

auto Sqrt3=1.7320508075688772;

class BoseHubbard
{
IQIndex s;

public:
BoseHubbard() {}
BoseHubbard(IQIndex I): s(I) {}
BoseHubbard(int n, Args const& args=Args::global())
{
// spinless boson (3 boson/site at most)
s=IQIndex(nameint("site=",n),
Index(nameint("Emp",n),1,Site),QN({0,1}),
Index(nameint("b1",n),1,Site),QN({1,1}),
Index(nameint("b2",n),1,Site),QN({2,1}),
Index(nameint("b3",n),1,Site),QN({3,1}));
}

IQIndex index() const {return s;}

IQIndexVal state(std::string const& state)
{
if (state=="Emp")
{
return s(1);
}
else if (state=="b1")
{
return s(2);
}
else if (state=="b2")
{
return s(3);
}
else if (state=="b3")
{
return s(4);
}
}

IQTensor op(std::string const& opname, Args const& args) const
{
auto sP=prime(s);

IQIndexVal
Em(s(1)),EmP(sP(1)),
b1(s(2)),b1P(sP(2)),
b2(s(3)),b2P(sP(3)),
b3(s(4)),b3P(sP(4));

IQTensor Op(dag(s),sP);
// boson single-site operator
if (opname=="Nb")
{
Op.set(b1,b1P,1);
Op.set(b2,b2P,2);
Op.set(b3,b3P,3);
}
else if (opname=="Ab")
{
Op.set(Em,b1P,1);
Op.set(b1,b2P,std::sqrt(2.));
Op.set(b2,b3P,std::sqrt(3.));
}
{
Op.set(b1,EmP,1);
Op.set(b2,b1P,std::sqrt(2.));
Op.set(b3,b2P,std::sqrt(3.));
}
else if (opname=="Id")
{
Op.set(Em,EmP,1);
Op.set(b1,b1P,1);
Op.set(b2,b2P,1);
Op.set(b3,b3P,1);
}
else
{
Error("Operator " + opname + " name not recognized !");
}

return Op;
}

};
}
#endif

Here I restricts that at most 3 bosons can lie in the same site. The main file is prsented as follows

// this file tests spinless Bose-Hubbard model
#include "itensor/all.h"
#include <typeinfo>
#include <iostream>
#include <fstream>
#include <vector>
#include <stdlib.h>
#include <cmath>

using namespace itensor;
using namespace std;

int main()
{
auto N = 20;
auto J = 1.0;
auto U = 0.5;
auto V = 1.0;

auto sites = BH(N);
auto ampo = AutoMPO(sites);

// site index must start from 1
for (int i = 1; i < N; i++) {
// hopping term
ampo += -J, "Adagb", i, "Ab", i+1;
ampo += -J, "Adagb", i+1, "Ab", i;
// onsite interaction
// ampo += U/2.0, "Nb", i, "Nb", i;
// ampo += -U/2.0, "Nb", i;
ampo += -U/2.0, "Adagb", i, "Ab", i;
// NN interaction
// ampo += V, "Nb", i, "Nb", i+1;
}
// ampo += U/2.0, "Nb", N, "Nb", N;
// ampo += -U/2.0, "Nb", N;
ampo += -U/2.0, "Adagb", N, "Ab", N;

auto Hamil = IQMPO(ampo);
auto state = InitState(sites);
for (int i = 1; i <= N; i++)
{
if (i%2 == 0){
state.set(i, "b2");
}
else {
state.set(i, "Emp");
}

}
auto psi=IQMPS(state);

// DMRG parameter
auto sweeps = Sweeps(10);
sweeps.maxm() = 160;
sweeps.cutoff() = 1E-12;
// perform DMRG algorithm
auto energy=dmrg(psi,Hamil,sweeps,{"Quiet",true});
println("Ground state energy = ",energy);
}
commented by (460 points)
Hi. I was just looking at it and, maybe is it because that for the first form the single particle state is a eigenstate, while the second is not?

@@H_{1}\left\vert 1\right\rangle =0\left\vert 1\right\rangle @@
@@H_{2}\left\vert 1\right\rangle =0 @@

Best,
Yixuan
commented by (70.1k points)
Hi Junjie,
I have a similar question to Yixuan above. While your two definitions look ok to me for the case of 0, 1, and 2 particles on a single site, are they still the same for 3 or more particles on a single site? Does your Hilbert space allow more than 2 particles on the same site?
commented by (650 points)
thanks, does that mean that I should rearrange the order of creation and annihilation operators ?
commented by (650 points)
edited
I have posted my siteset file of bose Hubbard model. I have checked it and it should be ok. Here I restricted that at single site, there should be no more than 3 bosons. As suggested by yixuan, I change the order of creation and annihilation operator to ensure that the single particle state is a eigenstate for both forms. However,  the problem remains.
commented by (460 points)
Hi Junjie,
I think if you rearranged the order of creation and annihilation then the eigenvalues would be different, like this (if you meant switching the middle two operators)

$$H_{1}\left\vert 1\right\rangle=(N_{b}N_{b}-N_{b})\left\vert 1\right\rangle =0\left\vert 1\right\rangle$$
$$H_{2}\left\vert 1\right\rangle=A_{b}^{\dagger }A_{b}A_{b}^{\dagger }A_{b}\left\vert 1\right\rangle =1\left\vert 1\right\rangle$$

Even the single particle state is a eigenstate for both forms, they are still different Hamiltonion.

Best,
Yixuan
commented by (650 points)
hi,yixuan. As you see in my main file, I have subtracted Nb so the total Hamiltonian should be the same.
commented by (70.1k points)
Hi Junjie,
That's helpful to see more of your code. What about the following: it looks like you might have defined "Ab" and "Adagb" backwards in your site set file. The convention about operators in ITensor is that the index with primelevel=0 is the one acting on the initial state, and the one with primelevel=1 is the index corresponding to the final state.
commented by (650 points)
thanks,miles. You are right ! It's my fault ... Now it is fixed.
commented by (220 points)
reshown
Can anyone elaborate this a bit please? What exactly need to change in order to correct the site set file for BH model?
Thanks.
commented by (70.1k points)
Adding onto this comment, Junjie - would you be willing to share your boson site set? I would be happy to add it to ITensor officially.