Customizing Hubbard model

Question

Dear Miles and ITensor Users,

I wish to implement the Bose-Hubbard model with a custom on-site interaction,
apart from the usual hopping and two-body interaction terms. The custom on-site interaction term is of the form Wb * (\delta(ni,3)), that is, 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. However I get wrong results (Energy) after running some tests. Below is how I implemented.
PS: Also, the keyword niter does not get recognized for some reason.

Any help will be really appreciated.

Thanks!

Regards-
Manpreet

using ITensors

function ITensors.space(::SiteType"Boson";
conserveqns=true)
if conserveqns
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 = productMPS(sites,state)
@show flux(psi0)
energy, psi = dmrg(H, psi0, sweeps)
println("\nGround State Energy = $energy")

end

Answer 1

Hi Manpreet,
The implementation of "dnb" looks correct to me. I think you are now just at the stage of doing various checks and debugging. Here all I can think to do to help you is to suggest some general debugging strategies I would try when studying a new model. Here are a few:

1. to test your Hamiltonian MPO, you can prepare some product states (using productMPS) and "sandwich" H between them using the overlap function to check whether H has the matrix elements you expect from your paper definition. Putting H between the same two states will check your interaction terms. Putting H between two different states will check your hopping terms.

2. for a small system size, you could convert your MPO to a dense matrix and diagonalize it to see if you get the result you expect

3. it may be that your model and Hamiltonian are correct, but that your DMRG run is failing to converge. I see you are using the noise term, which is good and may help a lot with this. You may also then need to carefully prepare your initial state to be ideally fairly close to the ground state. But this may not be the problem.

4. Are there some limits where you know the exact answer and can tune to these?

5. Does everything work correctly when you turn the interactions off? When you only have the more usual kind of interactions?

Miles

Comments

Hi Miles

Thank you for your suggestions. I will certainly try out those given in points 1, 2, 4.
Regarding 3:
I wanted to experiment with 'niter' as well but for some reason this keyword is not recognized.
Regarding 5:
With usual interactions I did not have much problem, although in that case I considered upto 4-5 bosons at a site.
I also notice that the value of flux remains unchanged irrespective of the changes I make in occupancy, state, etc.