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

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 :

//
// 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);
}
{
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);
}
{
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

}


+1 vote
answered Feb 18 by (18,230 points)

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

commented Feb 18 by (340 points)
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
commented Feb 19 by (340 points)
reshown Feb 20 by miles
By the way, I try your siteset file  ("nesites.h"). It has the same probelms as mine.
commented Feb 20 by (18,230 points)
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 -
commented Feb 20 by (340 points)
edited Feb 21
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.
commented Feb 22 by (18,230 points)
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.