Problems with on site interactions Bose-Hubbard Julia

DMRG and Numerical Methods
,
1

Hello,

I’m trying to implement a dmrg code for the computation of the GS of the Bose Hubbard model in Julia. The code that I’m using is the following one

using ITensors

let
  L = 4
  Np = 2
  miu = 1.0
  U = 1.0
  J = 1.0
  sites = siteinds("Boson", L; conserve_qns=true, conserve_number=true)


  os = OpSum()


 
  os = OpSum()
  for b in 1:(L - 1)
    os += -miu - U*0.5, "N", b
    os += U*0.5, "N", b, "N", b
    os += -J, "A", b, "Adag", b + 1
    os += -J, "Adag", b, "A", b + 1
  end
  os += -miu - U*0.5, "N", L
  os += U*0.5, "N", L, "N", L
  os += -J, "A", 1, "Adag", L
  os += -J, "Adag", 1, "A", L


  println(os)
  
  H = MPO(os, sites)

  states = ["0" for j=1:L]
    
    for j in 1:Np
        states[j] = "1"
    end
    
    psi0 = randomMPS(sites,states)



  sweeps = Sweeps(5)
  setmaxdim!(sweeps, 100,100,100,200,200)
  setcutoff!(sweeps, 1E-13)
  setnoise!(sweeps,1E-7,1E-10,1E-5,0,1E-8)

  energy, psi = dmrg(H,psi0, sweeps)

  println("\nGround State Energy = $energy")



  return
end

The problem is that my results seems to be independent on the choice of the interaction strength U. If I change the other parameters the energy changes, but if change U no. So, my impression is that for the computer “U0.5, “N”, b" and "U0.5, “N”, b, “N”, b” are the same, thus in the Hamiltonian they cancel each other and there isn’t the dependence on the interaction.

Anyway, this type of formalism (“U*0.5, “N”, b, “N”, b”) is the one that is typically used to write the on-site interaction at least with the C++ version. Does the Julia version needs something different to write properly the on site interaction? Or am I just doing something wrong somewhere?

Kind regards,
Francesco

2

Hi Francesco,
I think the issue here is what you said, about the “N” operator in this case squaring to 1. I don’t think it’s a Julia vs. C++ issue, but just a question about the mathematics of boson operators.

The key thing here is that “Boson” sites by default are “hard core” bosons, meaning they have a maximum occupancy of 1. You can check this by printing out one of your site indices:

julia> sites[1]
(dim=2|id=306|"Boson,Site,n=1") <Out>
 1: QN("Number",0) => 1
 2: QN("Number",1) => 1

Given your Hamiltonian definition, I’m guessing you are intending to study “soft core” bosons where theoretically an infinite number can occupy the same site. In practice this means you should work with a finite maximum occupancy and check that your results are converged in the limit of increasing this number.

To set a larger maximum occupancy, you can pass the dim keyword like this:

sites = siteinds("Boson", L; conserve_qns=true,dim=4)

(You don’t need to also pass conserve_number=true if you’ve already set conserve_qns=true, because conserve_qns=true means “conserve all possible QN’s” which includes number for bosons.)

3

Dear Miles,

Thank you, I didn’t realize that the dim=2 is set by default. Now it works.

Thank you,
Francesco