Implementing Jaynes–Cummings–Hubbard model using Julia

Question

Hello,

I am working with Julia version of ITensor and would like to know if Jaynes-Cummings-Hubbard (JCH) model can be implemented in it.

Cheers,
Manpreet

Edit : The JCH model is given, for example, in this article https://arxiv.org/abs/1611.06404v2

Comments

Hi, could you please provide a reference with the definition of the exact Hamiltonian you wish to implement? (Please see also the question guidelines here: http://itensor.org/support/2003/guidelines-for-asking-questions)

Thanks,
Miles

---

Hi,
I have edited my question to add a reference for example. To begin with, I wish to implement it for a linear chain and without nearest-neighbour interaction.

Regards-
Manpreet

Answer 1

Hi, so yes you can definitely implement this model, but with one major limitation, which is that for any degree of freedom which is a boson (such as a photon mode), you'll have to limit the maximum number of states of the Hilbert space of that boson. This is just because ITensor and tensor network methods generally (barring some very specialized approaches) are formulated in tensor products of finite-dimensional vector spaces.

We have not yet defined a boson degree of freedom that gets included with ITensor, but we have made it straightforward to define custom site types or Hilbert spaces. For a detailed walkthrough of how to do this, you can see the following two example pages:
http://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_basic
http://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_qns

Once you define overloads of the op function which specify operator names associated to your custom site type or Hilbert space, you can use these within the AutoMPO feature of ITensor to include your custom operator names, similar to other examples of using AutoMPO where there are operator names like "Sz", "S+", "S-" for the case of spin Hilbert spaces.

I'd recommend testing everything out on the smallest possible system and comparing to (numerically) exact calculations to check that it's all working the way you expect.

If you do end up making a boson Hilbert space definition, please consider contributing it to our library by sending us a pull request on Github. We can work with you to make it ready to include into the library if you think it would be a good idea.

Best regards,
Miles

Comments

The recent ITensor paper also has some information about AutoMPO and making custom Hilbert spaces if you find it helpful:
https://arxiv.org/abs/2007.14822

---

Hello Miles,

Thank you for the reply, the instructions were really helpful.
To get more familiar with how ITensor (and Julia) work, I am trying to write a code for simple Bose-Hubbard model first. I followed the instructions given on the example pages and from the exthubbard example file. However I am getting some errors while running the code, which seem to be related to how randomMPS is initialized. Can you please help me with it ?

Thanks!

Best regards,
Manpreet


Code:

using ITensors

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

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

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

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

s = siteind("S=0"; conserve_qns=true)

nb = op("nb",s)
#@show nb

adop = op("adop",s)
#@show adop

aop = op("aop",s)
#@show aop

let
  N      = 20
  Nb     = 5
  maxocc = 4
  tb1    = 1.0
  tb2    = 0.0
  Ub     = 1.0

  sites = siteinds("S=0",N; conserve_qns=true)
  ampo = AutoMPO()
  for i=1:N
    ampo += +Ub/2.0,"nb*nb",i
    ampo += -Ub/2.0,"nb",i
  end

  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

  H = MPO(ampo,sites)

  sweeps = Sweeps(6)
  maxdim!(sweeps,100,200,400,800)
  cutoff!(sweeps,1E-12)
  @show sweeps

  state = ["Emp" for n=1:N]
  p = Nb
  state = [0 for n=1:N]
  for i=N:-1:1
   if p >= maxocc
      state[i] = maxocc
      p -= maxocc
   elseif p > 0
      state[i] = p
      p = 0
   end
  end
  @show state

  psi0 = randomMPS(sites,state,1)
#  @show flux(psi0)
#  energy,psi = dmrg(H,psi0,sweeps)
#  println("\nGround State Energy = $energy")

end

---

Hi, good work on this. Definitely make sure to test your DMRG results out on small systems or compare to exact results when you get everything running.

The error I got when running your code had to do with how the array `state` was defined. For your custom site type, there are no mappings of string state names to Index values (so mapping of names like "Emp" to an empty site or similar). But this is not a problem since you can make the `state` array just be an array of integers, as you are doing.

The error is coming from the fact that you are starting the `state` integers from zero. Because Julia is a1-indexed language, a state of 1 means the first setting of an Index, which here you have defined (in your `space` function) to a state with 0 bosons. So state[i] = 1 means an empty site. Then state[i] = 2 means one boson, and state[i] = 3 means two bosons etc. You'll need to fill your state array that way, writing the loop in such a way that you get the correct number of bosons you want.

Optionally, once you have the `state` array correct, try raising the bond dimension argument at the end of `randomMPS` from 1 to another value such as 2 or 8. This will scramble up the state which can help to make it a better initial state for DMRG. (A very specfic low-entangled state, like a product state, can make DMRG hard to converge since it has a poor basis and can trap DMRG in a local minimum.)

Best,
Miles

---

Hello Miles,

Many thanks for your help and valuable suggestions. I was able to debug the code and now it seems to work properly (attached). I think the selection/design of initial state can still be improved and I am working on it. Please do let me know if you have any suggestions in this regard or in general for overall improvement.

Best regards-
Manpreet

 using ITensors

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

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

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

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

s = siteind("S=0"; conserve_qns=true)
#@show s

nb = op("nb",s)
#@show nb

adop = op("adop",s)
#@show adop

aop = op("aop",s)
#@show aop

let
  N      = 20
  Nb     = 10
  maxocc = 4
  mltplt = maxocc+1
  tb1    = 1.0
  tb2    = 0.0
  Ub     = 1.0

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

  ampo = AutoMPO()
  for i=1:N
    ampo += +Ub/2.0,"nb*nb",i
    ampo += -Ub/2.0,"nb",i
  end

  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

  H = MPO(ampo,sites)

  sweeps = Sweeps(10)
  maxdim!(sweeps,20,200,400,800,1200)
  cutoff!(sweeps,1E-12)
  @show sweeps

  p = Nb
  state = [1 for n=1:N]
  for i=N:-1:1
   if p >= maxocc
      state[i] = mltplt
      p -= maxocc
   elseif p > 0
      state[i] = p+1
      p = 0
   end
  end

  psi0 = randomMPS(sites,state,20)
  @show flux(psi0)
  energy,psi = dmrg(H,psi0,sweeps)
  println("\nGround State Energy = $energy")

end

---

Your code looks great at a quick look. So is it giving results which match to some other exact results that you know of?

My one comment would be that you technically don't have to name your sites "S=0" unless that is a name you actually prefer. You could change the name from "S=0" to "Boson" instead, say, or any other string (up to 8 characters) that you like.

Miles

---

Hello

Yes, I checked the results it with another code and found that ground state energy matches well (although it takes more sweeps, maybe initial state is  not good?). I still need to check local densities, correlations, etc.
Thanks for the suggestion. I will change the naming convention, I don't have any particular preference for it.

Manpreet

---

Hello again Miles

I have another question regarding the Bose-Hubbard model.
I made some changes to the previous version of the code to include an on-site
Wb * (\delta(ni,3)) kind of interaction in the Hamiltonian. That is to say, this interaction has a finite value only if there are 3 atoms at a particular site.  Additionally I wish to keep a maximum of 6 bosons at a site, with particle number conserved.
After running some tests I find that I am getting wrong results (Energy).
Below is how I implemented. It will be great if you can let me know where am I making a mistake.
PS: Also, the keyword niter does not get recognized for some reason.

Thanks!

Regards-
Manpreet

using ITensors

function ITensors.space(::SiteType"Boson";
                        conserve_qns=true)
  if conserve_qns
    return [QN("nb",0)=>1,QN("nb",1)=>1,
            QN("nb",2)=>1,QN("nb",3)=>1,
            QN("nb",4)=>1,QN("nb",5)=>1,
            QN("nb",6)=>1]
  end
  return 5
end

function ITensors.op!(Op::ITensor,
                      ::OpName"nb",
                      ::SiteType"Boson",
                      s::Index)
  Op[s'=>2,s=>2] = 1
  Op[s'=>3,s=>3] = 2
  Op[s'=>4,s=>4] = 3
  Op[s'=>5,s=>5] = 4
  Op[s'=>6,s=>6] = 5
  Op[s'=>7,s=>7] = 6
end

function ITensors.op!(Op::ITensor,
                      ::OpName"dnb",
                      ::SiteType"Boson",
                      s::Index)
  Op[s'=>4,s=>4] = 1
end

function ITensors.op!(Op::ITensor,
                      ::OpName"aop",
                      ::SiteType"Boson",
                      s::Index)
  Op[s'=>1,s=>2] = sqrt(1)
  Op[s'=>2,s=>3] = sqrt(2)
  Op[s'=>3,s=>4] = sqrt(3)
  Op[s'=>4,s=>5] = sqrt(4)
  Op[s'=>5,s=>6] = sqrt(5)
  Op[s'=>6,s=>7] = sqrt(6)
end

function ITensors.op!(Op::ITensor,
                      ::OpName"adop",
                      ::SiteType"Boson",
                      s::Index)
  Op[s'=>2,s=>1] = sqrt(1)
  Op[s'=>3,s=>2] = sqrt(2)
  Op[s'=>4,s=>3] = sqrt(3)
  Op[s'=>5,s=>4] = sqrt(4)
  Op[s'=>6,s=>5] = sqrt(5)
  Op[s'=>7,s=>6] = sqrt(6)
end

s = siteind("Boson"; conserve_qns=true)
#@show s

nb = op("nb",s)
#@show nb

dnb = op("dnb",s)
#@show dnb

adop = op("adop",s)
#@show adop

aop = op("aop",s)
#@show aop

let
  N      = 20
  Nb     = 25
  maxocc = 6
  mltplt = maxocc+1
  tb1    = 0.053
  tb2    = 0.0
  Ub     = 1.0
  Wb     = -1.5

  sites = siteinds("Boson",N; conserve_qns=true)

  ampo = AutoMPO()
  for i=1:N
    ampo += +Ub/2.0,"nb*nb",i
    ampo += -Ub/2.0,"nb",i
    ampo += Wb,"dnb",i
  end

  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

  H = MPO(ampo,sites)

  sweeps = Sweeps(20, [
        "maxdim" "mindim" "cutoff" "noise"
        10       10       1e-12    1E-7
        10       10       1e-12    1E-7
        20       10       1e-12    1E-7
        20       10       1e-12    1E-7
        50       10       1e-12    1E-7
        50       10       1e-12    1E-7
        100      20       1e-12    1E-8
        100      20       1e-12    1E-8
        200      20       1e-12    1E-10
        200      20       1e-12    1E-10
        400      20       1e-12    1E-10
        400      20       1e-12    1E-10
        800      20       1e-12    0
        800      20       1e-12    0
    ])
  @show sweeps

  p = Nb
  state = [0 for n=1:N]
  q = div(p,N)
  if q > maxocc
    println("\nError: MaxOcc too small.")
  end
  if q >= 1
     state = [q for n=1:N]
     p = p - q * N
  end
  i = div(N - p, 2)
  while p >= 1
    state[i+1] = state[i+1] + 1
    p = p - 1
    i = i + 1
  end
  @show state

#  psi0 = randomMPS(sites,state,20)
  psi0 = productMPS(sites,state)
  @show flux(psi0)
  energy, psi = dmrg(H, psi0, sweeps)
  println("\nGround State Energy = $energy")

end

---

Hi, since this is not so related to the comments aboe, could you please post it as a new question?

Thanks, Miles