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.));
}
else if (opname=="Adagb")
{
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, "Adagb", i, "Ab", i;
ampo += -U/2.0, "Adagb", i, "Ab", i;
// NN interaction
// ampo += V, "Nb", i, "Nb", i+1;
ampo += V, "Adagb", i, "Ab", i, "Adagb", i+1, "Ab", i+1;
}
// ampo += U/2.0, "Nb", N, "Nb", N;
// ampo += -U/2.0, "Nb", N;
ampo += U/2.0, "Adagb", N, "Ab", N, "Adagb", N, "Ab", 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);
}