# different results for Hamiltonian written in different form

+1 vote
edited Dec 16, 2017

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

ampo += U, "Adagb", i, "Adagb", i, "Ab", i, "Ab", i;


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 Dec 15, 2017 by (170 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 Dec 16, 2017 by (19,870 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 Dec 16, 2017 by (490 points)
thanks, does that mean that I should rearrange the order of creation and annihilation operators ?
commented Dec 16, 2017 by (490 points)
edited Dec 16, 2017
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 Dec 17, 2017 by (170 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 Dec 17, 2017 by (490 points)
hi,yixuan. As you see in my main file, I have subtracted Nb so the total Hamiltonian should be the same.
commented Dec 17, 2017 by (19,870 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 Dec 17, 2017 by (490 points)
thanks,miles. You are right ! It's my fault ... Now it is fixed.