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