Problems about 1-D Bose-Fermi mixture when using DMRG

Question

I have written a SiteSet file to describe a 1-D Bose-Fermi mixture system. Boson lies on even site while fermion lies on odd site. Boson is spinless and fermion is spin polarized. Boson is regular soft and we restrict that no more than 3 boson can lie on the same site. (truncation is 3) All works fine, including AutoMPO and IQMPO, and dmrg function works without faults. The initial state is set to be Mott-insulator for boson and half-filling for fermion. But the result of DMRG is obviously wrong including the following questions:

(1). No matter how many the number of site is, dmrg works so quickly. (less than 0.1 second) It seems that no calculation was done in dmrg;
(2). The vN Entropy at the center bond is always zero;

PS: the Hamiltonian is presented as follows

The SiteSet file I defined is presented as follows :

//
// Distributed under ITensor Library License, Version 2.x
// Author: Junjie Chen
// Date: 2017/02/14
//
#ifndef __ITENSOR_BOSEFERMIMIX__H
#define __ITENSOR_BOSEFERMIMIX__H
#include "itensor/mps/siteset.h"
#include <cmath>

namespace itensor
{
    class BoseFermiMix;
    using BF=BasicSiteSet<BoseFermiMix>;

class BoseFermiMix
{
    IQIndex s;

public:
    BoseFermiMix() {}
    BoseFermiMix(IQIndex I): s(I) {}
    BoseFermiMix(int n, Args const& args=Args::global())
    {
        if(n%2==1) // spinless fermion
        {
            s=IQIndex(nameint("site=",n),
                Index(nameint("Em",n),1,Site), QN("Sz=",0,"Nf=",0), // vaccum
                Index(nameint("f1",n),1,Site), QN("Sz=",0,"Nf=",1) // 1 fermion
                );  
        }
        else // spinless boson (3 boson/site at most)
        {
            s=IQIndex(nameint("site=",n),
            Index(nameint("Em",n),1,Site),QN("Sz=",0,"Nb=",0),
            Index(nameint("b1",n),1,Site),QN("Sz=",0,"Nb=",1),
            Index(nameint("b2",n),1,Site),QN("Sz=",0,"Nb=",2),
            Index(nameint("b3",n),1,Site),QN("Sz=",0,"Nb=",3)
            );
        }
    }

    IQIndex index() const {return s;}

    IQIndexVal state(std::string const& state)
    {
        if (state=="Em")
        {
            return s(1);
        }
        else if (state=="f1" || 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)),
        f1(s(2)),f1P(sP(2)),
        b1(s(2)),b1P(sP(2)),
        b2(s(3)),b2P(sP(3)),
        b3(s(4)),b3P(sP(4));

        IQTensor Op(dag(s),sP); // operator with 1 arrow in and 1 arrow out

        // Fermion single-site operator
        if (opname=="N" || opname=="n")
        {
            Op.set(f1,f1P,1);
        }
        else if (opname=="C")   
        {
            Op.set(Em,f1P,1);
        }
        else if (opname=="Cdag")
        {
            Op.set(f1,EmP,1);
        }
        else if (opname=="A")
        {
            Op.set(Em,f1P,1);
        }
        else if (opname=="Adag")
        {
            Op.set(f1,EmP,1);
        }
        else if (opname=="F" || opname=="FermiPhase")
        {
            Op.set(Em,EmP,1);
            Op.set(f1,f1P,-1);
        }
        else if (opname=="projEmp")
        {
            Op.set(Em,EmP,1);   
        }
        else if (opname=="projf1")
        {
            Op.set(f1,f1P,-1);
        }
        // Boson single-site operator
        else 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,Sqrt2);
            Op.set(b2,b3P,Sqrt3);
        }
        else if (opname=="Adagb")
        {
            Op.set(b1,EmP,1);
            Op.set(b2,b1P,Sqrt2);
            Op.set(b3,b2P,Sqrt3);
        }
        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

The main file I wrote is presented as follows :

#include "itensor/all.h"

using namespace itensor;

int main()
{
 auto N=200;
 auto Jb=1.0;
 auto Jf=1.0;
 auto Ub=10.0;
 auto Vbf=0.2;

 auto sites=BF(N);
 auto ampo=AutoMPO(sites);
 // create 1-D Bose-Fermi mixture Hamiltonian
 for (int i=1; i<N-1; ++i)
 {
   if (i%2==1)
 {
    ampo += -Jf, "Cdag", i, "C", i+2;
    ampo += -Jf, "Cdag", i+2, "C", i;
  }
  else if (i%2==0)
  {
    ampo += -Jb, "Adagb", i, "Ab", i+2;
    ampo += -Jb, "Adagb", i+2, "Ab", i;
  }
}
 for (int i=1; i<N; ++i)
{
   if (i%2==1)
{
  ampo += Vbf, "N", i, "Nb", i+1;
}
 else
{
  ampo += Vbf, "Nb", i, "N", i+1;
}
}
  for (int i=1; i<=N; ++i)
 {
   if (i%2==0)   
 {
   ampo += Ub/2.0, "Nb*Nb", i;
   ampo += -Ub/2.0, "Nb", i;
 }
}

 auto Hamil=IQMPO(ampo);
 // initialize the state
 auto state=InitState(sites);
 for (int i=1; i<=N; i++)
 {
   if (i%4==1)
 {
  state.set(i,"f1");
 }
  else if (i%4==3)
{
  state.set(i,"Em");
}
  else if (i%4==0)
{
  state.set(i,"b1");
}
  else if (i%4==2)
{
  state.set(i,"b1");
}
}
 auto psi=IQMPS(state);
 // DMRG parameter
 auto sweeps=Sweeps(1);
 sweeps.maxm()=100;
 sweeps.cutoff()=1E-12;
 // perform DMRG algorithm
 auto energy=dmrg(psi,Hamil,sweeps,{"Quiet",true});
 println("Ground state energy = ",energy);

 // measurement (particle number on each site)
 for (int i=1; i<=N; ++i)
 {
   if (i%2==0)
 {
  psi.position(i);
  IQTensor ket=psi.A(i);
  IQTensor bra=dag(prime(ket,Site));
  IQTensor nbi=sites.op("Nb",i);
  auto numB=(bra*nbi*ket).real();
  println("The ",i,"-th Boson site particle numer = ",numB);
}
} // boson

 for (int i=1; i<=N; ++i)
 {
   if (i%2==1)
{
  psi.position(i);
  IQTensor ket=psi.A(i);
  IQTensor bra=dag(prime(ket,Site));
  IQTensor nfi=sites.op("N",i);
  auto numF=(bra*nfi*ket).real();
  println("The ",i,"-th fermion site particle numer = ",numF);
}
} // fermion

}

Answer 1

Hi Junjie,
Your code looks good (sorry for not answering already via email - it's been a busy week). I think posting on here is a good idea anyway so people can learn from your code.

There may be a few different bugs causing issues; it's especially tricky to make new fermionic site sets other than the standard ones we provide but I can help you look into that later.

For now I see at least one bug, which is that the way you define the "F" operator gives a minus sign if there is one particle on a site regardless of whether its a bosonic or fermionic site (even though you use f1, it's (correctly) defined identically to b1). So you should put a flag (like a boolean) into your definition of BoseFermiMix that 'remembers' whether that site is fermionic or bosonic, then check this flag when making the "F" operator and put either a -1 or +1 as appropriate.

Try that and let me know if it still doesn't work and we can keep looking at things.

Miles

Comments

thanks, miles. I shall try to fix it according to your nice answer. But another problem comes out : in the Hamiltonian, fermion hopping term is something as Cdag_{i} * C_{i+2}, i.e. not nearest-neighbor hopping. Dose AutoMPO handle it automatically ? thanks

---

By the way, I try your siteset file  ("nesites.h"). It has the same probelms as mine.

---

Yes, AutoMPO does automatically handle second-neighbor and further-neighbor hopping terms, and other similar terms. You have to put in the minus signs just as you would on paper, as I see you did correctly in your main file.

About nesites.h, what problem are you finding that it has? For that site set, both types of sites are fermionic. We have been using it successfully on our project. Thanks -

---

Hi, miles. I'm sorry that I wrote a wrong Hamiltonian of spin-1/2 two leg ladder model. Now I fix it and your "netsite.h" works pretty fine. But I still have no idea about my problem. In fact, I don't know where and how to add some flags into my code to distinguish fermion site and boson site. I guess the flag should be passed by using "Args" objects but I find it hard to pass it to "BasicSiteSet" class which is important to build a whole sites from a single site class . Should I modify the "BasicSiteSet" class or do something else ? I may post my modified code in a new question. Thank you so much.

---

Hi Junjie, I see that you figured it out in the latest code you sent, but just for anyone else reading these comments I wanted to post the answer to your last question.

As you noticed, the BasicSiteSet class template passes an integer (the site number) to the constructor of the local site object (the object you named BoseFermiMix). So you can use this site number to determine fermion and boson sites and set a boolean flag (which you have to define) inside the BoseFermiMix object as a private internal variable.