# How to set total number of particles in bose-hubbard model?

+1 vote

Hello everyone!

I'd like to ask a question about how to set the total number of particles in bose-hubbard model without controlling the particle number of every site by using state.set(i,"number").

It's Hamiltonian is
$\hat{H}=-t\sum{i=1}^L(\hat{b}i^\dagger \hat{b}{i+1}+H.c.)+\frac{U}{2}\sum{i=1}^L\hat{n}i(\hat{n}i-1)-\mu\sum{i=1}^L\hat{n}i$

I want to calculate its phase diagram through its ground-state energy difference when one particle was added or removed.
$\mu(N)_+=E(N+1)-E(N)$
$\mu(N)_-=E(N)-E(N-1)$

Here is my code by C++ version. I know how to control the particle number on every site but I just want to control total number of particles not every site.

  double J=atof(argv);
double miu=atof(argv);
double U=atof(argv);
double V=atof(argv);
int L = atof(argv);

auto sites=Boson(L,{"MaxOcc=",5,"ConserveQNs",true,"ConserveNb",false});

auto sweeps = Sweeps(8);
sweeps.maxdim() = 40,80,400,400,800,800,1000,1000;
sweeps.cutoff() = 1E-16;

auto ampo = AutoMPO(sites);
for (int i = 1;i < L; i += 1)
{
ampo += -miu - U/2, "N", i;
ampo += U/2, "N", i,"N",i;
ampo += -J, "A", i, "Adag", i + 1;
ampo += -J, "Adag", i, "A", i + 1;
ampo += V, "N", i, "N", i + 1;
}
ampo += -miu - U/2, "N", L;
ampo += U/2, "N", L,"N",L;
/*ampo += -J, "A", 1, "Adag", L;
ampo += -J, "Adag", 1, "A", L;
ampo += V, "N", 1, "N", L;
*/
auto H = toMPO(ampo);

auto state = InitState(sites);
/*for (int i : range1(L))
{
state.set(i, "3");
}
*/
PrintData(sweeps);

auto psi0 = randomMPS(state);
auto [energy, psi] = dmrg(H, psi0, sweeps, {"Quiet=",true});
fprintf(fp,"%f,%f\n",miu,energy);
fclose(fp);
fp=fopen(file_name1.str().c_str(),"w");
//fprintf(fp,"the number of on-site particles\n");
float Numall = 0;
for( auto i : range1(L))//local density
{
psi.position(i);

auto ket = psi(i);
auto bra = dag(prime(ket,"Site"));

auto Numop = op(sites,"N",i);
auto Numi = elt(bra*Numop*ket);
Numall = Numall + Numi;
fprintf(fp,"%i,%.12f\n",i,Numi);
}
fclose(fp);
//calculating the expectation of correlation function
fp=fopen(file_name2.str().c_str(),"w");
fprintf(fp,"i,j,lgj,corf\n");
auto i = 1;
for(int j = 2;j < L+1; ++j)//measure a correlation function
//at site i and j
{
//make the operators
auto corfop_j = op(sites,"A",j);

psi.position(i);
auto psidag=dag(psi);
auto C = prime(psi(i),l_i)*corfop_i;
C *=prime(psidag(i),"Site");
for(int k = i+1; k<j; ++k)
{
C *= psi(k);
C *= psidag(k);
}
C *= prime(psi(j),lj)*corfop_j;
C *= prime(psidag(j),"Site");
auto corf_1j=elt(C);
fprintf(fp,"%i,%i,%.6f,%.12f\n",i,j,log10(j),corf_1j);
}

fclose(fp);
return 0;
}


+1 vote

Hi, so ITensor DMRG should already be doing what you want, namely making a state which has a fixed total particle number, not just a fixed number per site. (If the number of particles on a site never changed, actually, then ITensor would give incorrect results almost always since most ground states do not have that property.)

When you use the InitState feature to make a state that does indeed have a fixed number of particles on each site, that is just a special case because it is a product state. But when you input this product state into the dmrg function, because of the hopping terms in the Hamiltonian, the particles will start to move to other sites. After the DMRG algorithm has converged, you can check that there will no longer be a fixed number of particles on each site by measuring the density operator on each site.

What the DMRG code does conserve is the total number of particles you set. So even though they move around, the total number won't change.

If you want to find the ground state with one fewer or one more particle, the way to do that is to make your initial state have one fewer or one more particle, then input this into DMRG and it will again keep the total number of particles fixed throughout.

Miles

commented by (180 points)
Oh, I get it!
So, you mean if I want to set the total number of particles, say 64 particles in 32 sites, just fix 2 particles per site by using 'state.set(i,2)'.
And if I want to add one additional particle or remove one particle, just make particle number on one of sites plus 1 or minus 1, like 'state.set(i,3)' or 'state.set(i,1)'.
One more question is that does the site I chose to  add or remove one particle would effect the final result?
Thanks a lot!
commented by (70.1k points)
Great - yes you have definitely got it.

Your question about which site to add or remove from is a good question. Ideally it should not affect the answer, but DMRG results can sometimes depend on the initial state that is used, meaning either that it can get "stuck" into a state which is not the absolute ground state of that particle number sector or that it can take a longer time to converge so that one might need to do more sweeps or other tricks to help the system converge.

So I would just try adding or removing the particle near the *center* of the system, then do as many sweeps as you can reasonably afford. Plot the density of the +1 and -1 particle states to see if they look fairly left-right symmetric. Then if the energy seems well converged with the number of sweeps and the densities look converged you should be in good shape. Of course if you know some exact limits where you have benchmark answers it is good to check against those.
commented by (180 points)
OK，I know what you mean.