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