0 votes
asked by (570 points)

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!

3 Answers

+1 vote
answered by (240 points)
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 (570 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 (570 points)
Could you please see my modification of spinless.h file?
0 votes
answered by (570 points)

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
        if(opname == "Adag")
            {
            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
        if(opname == "Adag")
            {
            Op.set(Emp, N1P,1); 
            Op.set(N1,N2P,1);   
            }
        /*
        if(opname == "A")
            {
            Op.set(Occ,EmpP,1);
            }   
        else
        if(opname == "Adag")
            {
            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 (570 points)
edited by
1) I wonder about

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);   
after    if(opname == "Adag")?
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
answered by (32.3k points)

Hi AbraDabra,
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
    if(opname == "Adag")
        {
        Op(Emp,OneP) = 1;
        Op(One,TwoP) = std::sqrt(2);
        Op(Two,ThreeP) = std::sqrt(3);
        }

If you'd like to post your code when you're done with a first draft or have some more specific questions please just post a comment or a new question -

Best regards,
Miles

Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...