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";

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 = productMPS(sites,state)

@show flux(psi0)

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

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

end