I'm following an article: arXiv:0808.1675. When I calculate the spin and charge gap of Fig.2, I can't get the correct answer,my codes is:

```
#include "itensor/all.h"
```

using namespace itensor;

int

main()

{

using Holstein =MixedSiteSet<ElectronSite,BosonSite>;

auto N=60;

auto U=1.0;

auto t=1.0;

auto eplison=1.25;

auto omega=5.0;

auto sites=Holstein(N,{"ConserveNf",true,"ConserveNb",false,"MaxOcc=",4});

```
auto ampo =AutoMPO(sites);
for(int j=1;j<=N-2;j=j+2)
{
ampo += -t,"Cdagup",j+2,"Cup",j;
ampo += -t,"Cdagdn",j+2,"Cdn",j;
ampo += -t,"Cdagup",j,"Cup",j+2;
ampo += -t,"Cdagdn",j,"Cdn",j+2;
}
for(int j=1;j<=N;j=j+2)
{
ampo += U,"Nup",j,"Ndn",j;
ampo += -sqrt(eplison*omega),"Nup",j,"A",j+1;
ampo += -sqrt(eplison*omega),"Nup",j,"Adag",j+1;
ampo += -sqrt(eplison*omega),"Ndn",j,"A",j+1;
ampo += -sqrt(eplison*omega),"Ndn",j,"Adag",j+1;
ampo += omega,"Adag",j+1,"A",j+1;
}
auto H=toMPO(ampo);
auto sweeps =Sweeps(20);
sweeps.noise()=1E-6,1E-6,1E-8,1E-12;
sweeps.maxdim()=10,20,100,100,200,400,800;
sweeps.cutoff()=1E-10;
auto state =InitState(sites);
for(auto j:range1(N))
{
if(j%4==1)state.set(j,"Dn");
else if(j%4==3)state.set(j,"Up");
else state.set(j,"0");
}
auto psi0=randomMPS(state);
auto [energy,psi]=dmrg(H,psi0,sweeps,{"Quiet=",true});
auto state1=InitState(sites);
for(auto j:range1(N))
{
if(j==1)state1.set(j,"UpDn");
else if(j%4==1&&j!=1)state1.set(j,"Dn");
else if(j%4==3)state1.set(j,"Up");
else state1.set(j,"0");
}
auto psi1=randomMPS(state1);
auto [energyadd,psi2]=dmrg(H,psi1,sweeps,{"Quiet=",true});
auto state2=InitState(sites);
for(auto j:range1(N))
{
if(j==3)state2.set(j,"Emp");
else if(j%4==1)state2.set(j,"Dn");
else if(j%4==3&&j!=3)state2.set(j,"Up");
else state2.set(j,"0");
}
auto psi3=randomMPS(state2);
auto [energyminus,psi4]=dmrg(H,psi3,sweeps,{"Quiet=",true});
auto state3=InitState(sites);
for(auto j:range1(N))
{
if(j==1)state3.set(j,"Up");
else if(j%4==1&&j!=1)state3.set(j,"Dn");
else if(j%4==3)state3.set(j,"Up");
else state3.set(j,"0");
}
auto psi5=randomMPS(state3);
auto [energyspin,psi6]=dmrg(H,psi5,sweeps,{"Quiet=",true});
printfln("spin gap=",energyspin-energy);
printfln("charge excitation=", energyadd+energyminus-2*energy);
return 0;
}
```

When I calculate N=60 in the code(actually the system is 30) , lambda=0.625, boson is 4. the spin gap is 0.351311, the charge gap is 0.442756. They are not correct. I try the other ways, but it doesn't work. Whether I ignore something?