# Implementing Jaynes–Cummings–Hubbard model using Julia

+1 vote
edited

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

commented by (70.1k points)

Thanks,
Miles
commented by (160 points)
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

+1 vote

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

commented by (70.1k points)
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
commented by (160 points)
Hello Miles,

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,
::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

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
end

for b=1:N-2
if b % 2 == 1
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 commented by (70.1k points) 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 commented by (160 points) edited by 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
commented by (70.1k points)
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
commented by (160 points)
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
commented by (160 points)
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,
::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

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
end

for b=1:N-2
if b % 2 == 1
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
commented by (70.1k points)
Hi, since this is not so related to the comments aboe, could you please post it as a new question?

Thanks, Miles