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[1]);
double miu=atof(argv[2]);
double U=atof(argv[3]);
double V=atof(argv[4]);
int L = atof(argv[5]);
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_i = op(sites,"Adag",i);
auto corfop_j = op(sites,"A",j);
psi.position(i);
auto psidag=dag(psi);
psidag.prime("Link");
auto l_i = leftLinkIndex(psi,i);
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);
}
auto lj = rightLinkIndex(psi,j);
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;
}
Looking forward to your help.