+1 vote
asked by (180 points)

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.

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");

auto psi0 = randomMPS(state);
auto [energy, psi] = dmrg(H, psi0, sweeps, {"Quiet=",true});
//fprintf(fp,"the number of on-site particles\n");
float Numall = 0;
for( auto i : range1(L))//local density

  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;
//calculating the expectation of correlation function
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);

auto psidag=dag(psi);
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);

return 0;

Looking forward to your help.

1 Answer

+1 vote
answered by (70.1k points)

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.

Hope that answers your question! If not please post a comment below.


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.
Your answers just solve the problem bothering me for days.
Thank you very much, Miles!
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.