State evolution in Bose-Hubbard model - Julia

ITensor Julia Questions
,
1

Hey ITensor team and users,

I am currently working on simulating the dynamics of the Bose-Hubbard model with nearest-neighbor interactions using ITensor in Julia. I have successfully implemented the time-independent version, but I am facing challenges in adapting it for the time-dependent case. I followed the tutorial provided on the ITensor website (MPSTimeEvolution) and incorporated other relevant snippets from discussions in this group.

As a test case, I am trying to evolve a density wave initial state (hardcore bosons) towards a superfluid state by applying only the hopping term, -t * (aįµ¢ā€  aā±¼ + aā±¼ā€  aįµ¢), to this initial state. However, the energy values I obtain suggest that my code isnā€™t behaving as expected.

Since I am relatively new to Julia and the dynamics part, I would greatly appreciate any guidance or hints. I have attached my code below for reference.

Additionally, I have a couple of specific questions:

  1. Does the term c = div(N,2) from the template handle systems with an odd number of sites properly?
  2. Between TEBD and TDVP, which method would be more suitable for my case? I am currently using TEBD but am unsure if switching to TDVP would improve results.

Thank you in advance for any help!

Best regards,
Aman

Code :

using ITensors, ITensorMPS

function ITensors.space(::SiteType"S=0";
                        conserve_qns=false)
  if conserve_qns
    return [QN("nb",0)=>1,QN("nb",1)=>1]
  end
  return 5
end

function ITensors.op!(Op::ITensor,
                      ::OpName"nb",
                      ::SiteType"S=0",
                      s::Index)
  Op[s'=>2,s=>2] = 1
end

function ITensors.op!(Op::ITensor,
                      ::OpName"aop",
                      ::SiteType"S=0",
                      s::Index)
  Op[s'=>1,s=>2] = sqrt(1)
end

function ITensors.op!(Op::ITensor,
                      ::OpName"adop",
                      ::SiteType"S=0",
                      s::Index)
  Op[s'=>2,s=>1] = sqrt(1)
end

s = siteind("S=0"; conserve_qns=true)
nb = op("nb",s)
adop = op("adop",s)
aop = op("aop",s)

let
  N      = 11
  Nb     = ((N-1)/2)+1
  maxocc = 1
  mltplt = maxocc+1
  tb1    = 1.0
  tb2    = 0.0
  Vb     = 0.0

  sites = siteinds("S=0",N; conserve_qns=true)

  ampo = AutoMPO()
  for b=1:N-1
    ampo += -tb1,"adop",b,"aop",b+1
    ampo += -tb1,"adop",b+1,"aop",b
  end

  for b=1:N-2
    if b % 2 == 1
      ampo += -tb2,"adop",b,"aop",b+2
      ampo += -tb2,"adop",b+2,"aop",b
    end
  end

  for i=1:N-1
    ampo += +Vb,"nb",i,"nb",i+1
  end

  H = MPO(ampo,sites)

  p = Nb
  state = [1 for n=1:N]
  for i=N:-1:1
   if p >=0
    if isodd(i)
       state[i] = 2
       p -= 1
    else
       state[i] = 1
    end
   end
  end

  psi0 = productMPS(sites,state)
  @show flux(psi0)
  psi = psi0
  energy = 0.0
  println("\nGround State Energy = $energy")
  nlocal = expect(psi,"nb")
  @show nlocal
  nncorr = correlation_matrix(psi,"nb","nb")
  println(nlocal)

  tau = 0.05
  ttotal = 20.0
  cutoff = 1E-8
  # Make gates (1,2),(2,3),(3,4),...
  gates = ITensor[]
  for j in 1:(N - 1)
    n1 = sites[j]
    n2 = sites[j + 1]
    hj = -tb1 * (op("adop", n1) * op("aop", n2) + op("aop", n1) * op("adop", n2))
    Gj = exp(-1.0im * tau / 2 * hj)
    push!(gates, Gj)
  end
  # Include gates in reverse order too
  # (N,N-1),(N-1,N-2),...
  append!(gates, reverse(gates))

  c = div(N, 2) # center site

  # Compute and print energy at each time step
  # then apply the gates to go to the next time
  for t in 0.0:tau:ttotal
    energy_dyn = inner(psi', H, psi)
    println("$t $energy_dyn")
    energy_squared = inner(H, psi, H, psi)
    variance = energy_squared - energy^2
    println("\nEnergy variance = $variance")

    tā‰ˆttotal && break

    psi = apply(gates, psi; cutoff)
    normalize!(psi)
  end
  nlocal = expect(psi,"nb")
  @show nlocal

  return
end
2

Hi Aman,
The code c = div(N,2) is the Julia integer division function. For even values of N, it returns N/2. For odd values of N, it returns (N-1)/2. I will let you decide if that behavior is what you need or not.

Note that there is no requirement that you use div(N,2) or even the variable c as we defined it. It was just shown there as an example of measuring a local property at a certain site near the center. In your modified code, it looks like you are not even using c, so you can just delete that line.

About TEBD versus TDVP, I would say the main reasons to use TDVP are when you are doing a system with complicated or long-range interaction such that you need an MPO, or you are doing a very expensive time evolution, where TDVP can offer the benefit of being more efficient. But if you can use TEBD, especially if you are doing a 1D system with short-range interactions, then I would say TEBD is the best ā€œsafeā€ choice because it automatically adapts the bond dimension very nicely after each step whereas with TDVP this can currently be a much trickier issue, and TDVP can sometimes give wrong answers unless you carefully check it and perform certain ā€œsubspace expansionā€ steps to help it to converge, which we currently only have partial support for.

So in short, I would recommend TEBD for your case, which I think is a 1D short-range case.

Finally, I am not sure why you are not getting the results you expect but please do carefully consider the role of finite-size and open-boundary effects, and also please explore various parameter settings like the maxdim and cutoff parameters controlling the truncation at each time step.

3

Hi Dr. Miles,
Could I ask a simple question about ITensor.jl for boson system under this post? Recently I switched my ITensor from C++ version to Julia version. In C++ version I defined the SiteType as sites = Boson(L, {"MaxOcc=", n, "ConserveQNs", true, "ConserveNb", false}) to make a cutoff on the maximum occupation to reduce the dimension. In Julia version I havenā€™t found such a setting for the cutoff from the guidance. I am wondering if there is a similar embedded commands for the maximum occupancy setting in Julia version? Or how could I set such a cutoff? Think you so much for your nice work.

4

Hi Ning,
Would you please post your question as a new topic / new post? I think it is a good question, but deserves its own separate post. Thanks.

Miles

1 Like