+1 vote
asked by (650 points)
edited by

Here I have a Bose-Hubbard model (with self-defined siteset file) and the density-density interaction can be written in two identical forms as follows

ampo += U, "Nb", i, "Nb", i;
ampo += -U, "Nb", i;

and

ampo += U, "Adagb", i, "Adagb", i, "Ab", i, "Ab", i;

this two forms of "ampo" should be the same because @@N_b=A^{\dagger}A@@. However, the final results of dmrg is quite different for this two forms. In fact, the first form gives the correct answer while the second one does not and it just converge very slowly.

Here I presents the siteset file for Bose Hubbard model

#ifndef __ITENSOR_BH_H
#define __ITENSOR_BH_H
#include "itensor/mps/siteset.h"    
#include <cmath>

namespace itensor
{
  class BoseHubbard;
  using BH=BasicSiteSet<BoseHubbard>;

  auto Sqrt3=1.7320508075688772;

  class BoseHubbard
   {
    IQIndex s;

  public:
    BoseHubbard() {}
    BoseHubbard(IQIndex I): s(I) {}
    BoseHubbard(int n, Args const& args=Args::global())
    {
      // spinless boson (3 boson/site at most)
      s=IQIndex(nameint("site=",n),
      Index(nameint("Emp",n),1,Site),QN({0,1}),
      Index(nameint("b1",n),1,Site),QN({1,1}),
      Index(nameint("b2",n),1,Site),QN({2,1}),
      Index(nameint("b3",n),1,Site),QN({3,1}));
    }

    IQIndex index() const {return s;}

    IQIndexVal state(std::string const& state)
    {
      if (state=="Emp")
      {
        return s(1);
      }
      else if (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)),
      b1(s(2)),b1P(sP(2)),
      b2(s(3)),b2P(sP(3)),
      b3(s(4)),b3P(sP(4));

      IQTensor Op(dag(s),sP);
      // boson single-site operator
      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,std::sqrt(2.));
        Op.set(b2,b3P,std::sqrt(3.));
      }
      else if (opname=="Adagb")
      {
        Op.set(b1,EmP,1);
        Op.set(b2,b1P,std::sqrt(2.));
        Op.set(b3,b2P,std::sqrt(3.));
      }
      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

Here I restricts that at most 3 bosons can lie in the same site. The main file is prsented as follows

// this file tests spinless Bose-Hubbard model
#include "itensor/all.h"
#include <typeinfo>
#include <iostream>
#include <fstream>
#include <vector>
#include <stdlib.h>
#include <cmath>

using namespace itensor;
using namespace std;

int main()
{
  auto N = 20;
  auto J = 1.0;
  auto U = 0.5;
  auto V = 1.0;

  auto sites = BH(N);
  auto ampo = AutoMPO(sites);

  // site index must start from 1
  for (int i = 1; i < N; i++) {
    // hopping term
    ampo += -J, "Adagb", i, "Ab", i+1;
    ampo += -J, "Adagb", i+1, "Ab", i;
    // onsite interaction
    // ampo += U/2.0, "Nb", i, "Nb", i;
    // ampo += -U/2.0, "Nb", i;
    ampo += U/2.0, "Adagb", i, "Ab", i, "Adagb", i, "Ab", i;
    ampo += -U/2.0, "Adagb", i, "Ab", i;
    // NN interaction 
    // ampo += V, "Nb", i, "Nb", i+1;
    ampo += V, "Adagb", i, "Ab", i, "Adagb", i+1, "Ab", i+1;
  }
  // ampo += U/2.0, "Nb", N, "Nb", N;
  // ampo += -U/2.0, "Nb", N;
  ampo += U/2.0, "Adagb", N, "Ab", N, "Adagb", N, "Ab", N;
  ampo += -U/2.0, "Adagb", N, "Ab", N;

  auto Hamil = IQMPO(ampo);
  auto state = InitState(sites);
  for (int i = 1; i <= N; i++)
  {
    if (i%2 == 0){
      state.set(i, "b2");
    }
    else {
      state.set(i, "Emp");
    }

  }
  auto psi=IQMPS(state);

  // DMRG parameter
  auto sweeps = Sweeps(10);
  sweeps.maxm() = 160;
  sweeps.cutoff() = 1E-12;
  // perform DMRG algorithm
  auto energy=dmrg(psi,Hamil,sweeps,{"Quiet",true});
  println("Ground state energy = ",energy);
}
commented by (460 points)
Hi. I was just looking at it and, maybe is it because that for the first form the single particle state is a eigenstate, while the second is not?

@@H_{1}\left\vert 1\right\rangle =0\left\vert 1\right\rangle @@
@@H_{2}\left\vert 1\right\rangle =0 @@

Best,
Yixuan
commented by (70.1k points)
Hi Junjie,
I have a similar question to Yixuan above. While your two definitions look ok to me for the case of 0, 1, and 2 particles on a single site, are they still the same for 3 or more particles on a single site? Does your Hilbert space allow more than 2 particles on the same site?
commented by (650 points)
thanks, does that mean that I should rearrange the order of creation and annihilation operators ?
commented by (650 points)
edited by
I have posted my siteset file of bose Hubbard model. I have checked it and it should be ok. Here I restricted that at single site, there should be no more than 3 bosons. As suggested by yixuan, I change the order of creation and annihilation operator to ensure that the single particle state is a eigenstate for both forms. However,  the problem remains.
commented by (460 points)
Hi Junjie,
I think if you rearranged the order of creation and annihilation then the eigenvalues would be different, like this (if you meant switching the middle two operators)

$$H_{1}\left\vert 1\right\rangle=(N_{b}N_{b}-N_{b})\left\vert 1\right\rangle =0\left\vert
1\right\rangle $$
$$H_{2}\left\vert 1\right\rangle=A_{b}^{\dagger }A_{b}A_{b}^{\dagger }A_{b}\left\vert 1\right\rangle
=1\left\vert 1\right\rangle $$

Even the single particle state is a eigenstate for both forms, they are still different Hamiltonion.

Best,
Yixuan
commented by (650 points)
hi,yixuan. As you see in my main file, I have subtracted Nb so the total Hamiltonian should be the same.
commented by (70.1k points)
Hi Junjie,
That's helpful to see more of your code. What about the following: it looks like you might have defined "Ab" and "Adagb" backwards in your site set file. The convention about operators in ITensor is that the index with primelevel=0 is the one acting on the initial state, and the one with primelevel=1 is the index corresponding to the final state.
commented by (650 points)
thanks,miles. You are right ! It's my fault ... Now it is fixed.
commented by (220 points)
reshown by
Can anyone elaborate this a bit please? What exactly need to change in order to correct the site set file for BH model?
Thanks.
commented by (70.1k points)
Adding onto this comment, Junjie - would you be willing to share your boson site set? I would be happy to add it to ITensor officially.

1 Answer

0 votes
answered by (70.1k points)

Hi Junjie, glad it is working now. (Answering here to make sure the post is officially answered. See discussion above.)

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.

Categories

...