spinless bosons

Hello!

How one can construct spinless boson basis?
Here there is only hard-core bosons http://itensor.org/docs.cgi?page=classes/spinless&vers=cppv2

And here there is only spinfull bosons http://itensor.org/docs.cgi?vers=cppv2&page=classes/hubbard

I want to consider 1D bose-hubbard model.

thank you!

+1 vote
edited by

Hi,

Spinless boson is tricky, but there are many ways to deal with it. See section V of U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).

The easiest thing one can do is just cutting the local Hilbert space to some max occupancy, e. g. double occupancy, as long as one is only interested in low energy physics. You just need to define your own siteset similar to the one already defined in the ITensor library, with the index changed to be e. g.

s = IQIndex{
nameint("site=",n),
Index(nameint("Emp",n),1,Site),QN("Nb=",0),
Index(nameint("N1",n),1,Site), QN("Nb=",1),
Index(nameint("N2",n),1,Site), QN("Nb=",2)
};


and changes the definition of the site operator accordingly.

commented by (640 points)
thx. could u please write more in detail "...changes the definition of the site operator accordingly."
commented by (240 points)
edited by
You just need to write a file similar to e. g. spinless.h, which is already defined in the ITensor library.
What you need to change:
1. Those "s = IQIndex{ ... }" lines.
2. The function "IQIndexVal state ...": now you have three cases in the if...else..., i.e. state == "0", state == "1", state == "2".
3. The function "IQTensor op...": since you are doing bosonic system, you no longer need those "C", "Cdag", "F"... and you need to change all other definitions of the bosonic operators. For example, "A" is changed to be:
if(opname == "A")
{
Op.set(N1,EmpP, 1);
Op.set(N2,N1P, sqrt(2));
}
with EmpP=sP(1), N1=s(2), N1P=sP(2), N2=s(3).
commented by (640 points)
Could you please see my modification of spinless.h file?

Could someone see my modification of spinless.h file to simulate spinless bosons in 1D Bose-Hubbard model?

#ifndef __ITENSOR_SPINLESS_H
#define __ITENSOR_SPINLESS_H
#include "itensor/mps/siteset.h"

namespace itensor {

class SpinlessBosons;

using SpinlessBos = BasicSiteSet<SpinlessBosons>;

class SpinlessBosons
{
IQIndex s;
public:

SpinlessBosons() { }
SpinlessBosons(IQIndex I) : s(I) { }
SpinlessBosons(int n, Args const& args = Args::global())
{
s = IQIndex{
nameint("SpinlessBos",n),
Index(nameint("Emp",n),1,Site),QN("Nb=",0),
Index(nameint("N1",n),1,Site), QN("Nb=",1),
Index(nameint("N2",n),1,Site), QN("Nb=",2)
};

}

IQIndex
index() const { return s; }

IQIndexVal
state(std::string const& state)
{

if(state == "Emp" || state == "0")
{
return s(1);
}
else
if(state == "N1" || state == "1")
{
return s(2);
}
else
if(state == "N2" || state == "2")
{
return s(3);
}
else
{
Error("State " + state + " not recognized");
}
return IQIndexVal{};
}

IQTensor
op(std::string const& opname,
Args const& args) const
{
auto sP = prime(s);
auto Emp  = s(1);
auto EmpP = sP(1);
auto N1  = s(2);
auto N1P = sP(2);
auto N2 = s(3);
auto N2P = sP(3);

auto Op = IQTensor(dag(s),sP);

if(opname == "N" || opname == "n")
{
Op.set(N1,N1P,1);
Op.set(N2,N2P,1);
}
else
if(opname =="A")
{
Op.set(N1, EmpP,1);
Op.set(N2,N1P,1);
}
else
{
Op.set(Emp, N1P,1);
Op.set(N1,N2P,1);
}

else
{
Error("Operator \"" + opname + "\" name not recognized");
}

return Op;
}
};

} //namespace itensor

#endif


the same code with comments (i.e. the code in spinless.h is commented out ) for comparison with the file spinless.h

#ifndef __ITENSOR_SPINLESS_H
#define __ITENSOR_SPINLESS_H
#include "itensor/mps/siteset.h"

namespace itensor {

//class SpinlessSite;
class SpinlessBosons;
//using Spinless = BasicSiteSet<SpinlessSite>;
using SpinlessBos = BasicSiteSet<SpinlessBosons>;

//class SpinlessSite
class SpinlessBosons
{
IQIndex s;
public:

//    SpinlessSite() { }
SpinlessBosons() { }
//    SpinlessSite(IQIndex I) : s(I) { }
SpinlessBosons(IQIndex I) : s(I) { }
//    SpinlessSite(int n, Args const& args = Args::global())
SpinlessBosons(int n, Args const& args = Args::global())
{
s = IQIndex{
nameint("SpinlessBos",n),
Index(nameint("Emp",n),1,Site),QN("Nb=",0),
Index(nameint("N1",n),1,Site), QN("Nb=",1),
Index(nameint("N2",n),1,Site), QN("Nb=",2)
};
/*
auto conserve_Nf = args.getBool("ConserveNf",true);
auto oddevenupdown = args.getBool("OddEvenUpDown",false);

if(!oddevenupdown) //usual case
{
auto q_occ = QN("Nf=",1);
if(not conserve_Nf) q_occ = QN("Pf=",1);

s = IQIndex{nameint("Spinless ",n),
Index(nameint("Emp ",n),1,Site),QN(),
Index(nameint("Occ ",n),1,Site),q_occ};

}
else
{
QN q_occ;
if(n%2==1) q_occ = QN("Sz",+1,"Nf=",1);
else       q_occ = QN("Sz",-1,"Nf=",1);
s = IQIndex{nameint("Spinless ",n),
Index(nameint("Emp ",n),1,Site),QN(),
Index(nameint("Occ ",n),1,Site),q_occ};

}*/
}

IQIndex
index() const { return s; }

IQIndexVal
state(std::string const& state)
{

if(state == "Emp" || state == "0")
{
return s(1);
}
else
if(state == "N1" || state == "1")
{
return s(2);
}
else
if(state == "N2" || state == "2")
{
return s(3);
}
else
{
Error("State " + state + " not recognized");
}
/*
if(state == "Emp" || state == "0")
{
return s(1);
}
else
if(state == "Occ" || state == "1")
{
return s(2);
}
else
{
Error("State " + state + " not recognized");
}*/
return IQIndexVal{};
}

IQTensor
op(std::string const& opname,
Args const& args) const
{
auto sP = prime(s);
auto Emp  = s(1);
auto EmpP = sP(1);
auto N1  = s(2);
auto N1P = sP(2);
auto N2 = s(3);
auto N2P = sP(3);

//auto Occ  = s(2);
//auto OccP = sP(2);

auto Op = IQTensor(dag(s),sP);

/*
if(opname == "N" || opname == "n")
{
Op.set(Occ,OccP,1);
}
else
if(opname == "C")
{
Op.set(Occ,EmpP,1);
}
else
if(opname == "Cdag")
{
Op.set(Emp,OccP,1);
}
else
*/
if(opname == "N" || opname == "n")
{
Op.set(N1,N1P,1);
Op.set(N2,N2P,1);
}
else
if(opname =="A")
{
Op.set(N1, EmpP,1);
Op.set(N2,N1P,1);
}
else
{
Op.set(Emp, N1P,1);
Op.set(N1,N2P,1);
}
/*
if(opname == "A")
{
Op.set(Occ,EmpP,1);
}
else
{
Op.set(Emp,OccP,1);
}
else
if(opname == "F" || opname == "FermiPhase")
{
Op.set(Emp,EmpP,1);
Op.set(Occ,OccP,-1);
}
else
if(opname == "projEmp")
{
Op.set(Emp,EmpP,1);
}
else
if(opname == "projOcc")
{
Op.set(Occ,OccP,1);
}
*/
else
{
Error("Operator \"" + opname + "\" name not recognized");
}

return Op;
}
};

} //namespace itensor

#endif

commented by (640 points)
edited

auto conserve_Nb = args.getBool("ConserveNb",true);

should I add something like that for conserving the number of particles?
Will the number of particles be conserved in the code written above?

2) I don't really understand what is

auto oddevenupdown = args.getBool("OddEvenUpDown",false);

and should I insert something similar in my file

3) Is it Ok to write
Op.set(N1,N2P,1);
the action of Adag on the sate with 1 particle should give the state with 2 particles multiplied by sqrt(2).
commented by (240 points)
edited by
Hi,
You can just ignore 1) and 2). For 3), You are right. I have corrected this typo in my answer.
Your first code looks good after changing "Op.set(N1,N2P,1)" to "Op.set(N1,N2P,std::sqrt(2))" and the corresponding one in opname ==  "A" (the correct is "Op.set(N2,N1P,std::sqrt(2))") and also in opname == "N" (the correct is "Op.set(N2,N2P,2)").
As Miles suggested below, you can test with max occupancy defined to be 2, 3, 4... to see if that makes any difference in your result.
+1 vote

This is a site set type we should offer. We're working on it!

But in the meantime, yes, you can create your own kind of site set for this case by copying one of the existing ones (such as "Spinless") and creating a new file and a new site set (perhaps boson.h and "class BosonSite" and "class Boson").

Though Mingru is correct that sometimes just having a maximum site occupancy of 1 is enough, in general that's not the case, very much depending on the particular Hamiltonian.

But for many cases, you can truncate the maximum occupancy to, say, 3 or 4 without incurring much error. You should raise and lower this maximum occupancy and check that your results are converged with respect to it.

For the operators you'll need to define, here is an example for maximum occupancy of 3:

ITensor
op(std::string const& opname,
Args const& args) const
{
auto sP = prime(s); //here "s" is the site index stored in this BosonSite object

IQIndexVal Emp(s(1)),
EmpP(sP(1)),
One(s(2)),
OneP(sP(2)),
Two(s(3)),
TwoP(sP(3)),
Three(s(4)),
ThreeP(sP(4));

if(opname == "N" || opname == "n")
{
Op(One,OneP) = 1;
Op(Two,TwoP) = 2;
Op(Three,ThreeP) = 3;
}
else
if(opname == "A")
{
Op(One,EmpP) = 1;
Op(Two,OneP) = std::sqrt(2);
Op(Three,TwoP) = std::sqrt(3);
}
else