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
}
```