Hey ITensor team and users,
I am currently working on simulating the dynamics of the Bose-Hubbard model with nearest-neighbor interactions using ITensor in Julia. I have successfully implemented the time-independent version, but I am facing challenges in adapting it for the time-dependent case. I followed the tutorial provided on the ITensor website (MPSTimeEvolution) and incorporated other relevant snippets from discussions in this group.
As a test case, I am trying to evolve a density wave initial state (hardcore bosons) towards a superfluid state by applying only the hopping term, -t * (aᵢ† aⱼ + aⱼ† aᵢ), to this initial state. However, the energy values I obtain suggest that my code isn’t behaving as expected.
Since I am relatively new to Julia and the dynamics part, I would greatly appreciate any guidance or hints. I have attached my code below for reference.
Additionally, I have a couple of specific questions:
- Does the term
c = div(N,2)from the template handle systems with an odd number of sites properly? - Between TEBD and TDVP, which method would be more suitable for my case? I am currently using TEBD but am unsure if switching to TDVP would improve results.
Thank you in advance for any help!
Best regards,
Aman
Code :
using ITensors, ITensorMPS
function ITensors.space(::SiteType"S=0";
conserve_qns=false)
if conserve_qns
return [QN("nb",0)=>1,QN("nb",1)=>1]
end
return 5
end
function ITensors.op!(Op::ITensor,
::OpName"nb",
::SiteType"S=0",
s::Index)
Op[s'=>2,s=>2] = 1
end
function ITensors.op!(Op::ITensor,
::OpName"aop",
::SiteType"S=0",
s::Index)
Op[s'=>1,s=>2] = sqrt(1)
end
function ITensors.op!(Op::ITensor,
::OpName"adop",
::SiteType"S=0",
s::Index)
Op[s'=>2,s=>1] = sqrt(1)
end
s = siteind("S=0"; conserve_qns=true)
nb = op("nb",s)
adop = op("adop",s)
aop = op("aop",s)
let
N = 11
Nb = ((N-1)/2)+1
maxocc = 1
mltplt = maxocc+1
tb1 = 1.0
tb2 = 0.0
Vb = 0.0
sites = siteinds("S=0",N; conserve_qns=true)
ampo = AutoMPO()
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
for i=1:N-1
ampo += +Vb,"nb",i,"nb",i+1
end
H = MPO(ampo,sites)
p = Nb
state = [1 for n=1:N]
for i=N:-1:1
if p >=0
if isodd(i)
state[i] = 2
p -= 1
else
state[i] = 1
end
end
end
psi0 = productMPS(sites,state)
@show flux(psi0)
psi = psi0
energy = 0.0
println("\nGround State Energy = $energy")
nlocal = expect(psi,"nb")
@show nlocal
nncorr = correlation_matrix(psi,"nb","nb")
println(nlocal)
tau = 0.05
ttotal = 20.0
cutoff = 1E-8
# Make gates (1,2),(2,3),(3,4),...
gates = ITensor[]
for j in 1:(N - 1)
n1 = sites[j]
n2 = sites[j + 1]
hj = -tb1 * (op("adop", n1) * op("aop", n2) + op("aop", n1) * op("adop", n2))
Gj = exp(-1.0im * tau / 2 * hj)
push!(gates, Gj)
end
# Include gates in reverse order too
# (N,N-1),(N-1,N-2),...
append!(gates, reverse(gates))
c = div(N, 2) # center site
# Compute and print energy at each time step
# then apply the gates to go to the next time
for t in 0.0:tau:ttotal
energy_dyn = inner(psi', H, psi)
println("$t $energy_dyn")
energy_squared = inner(H, psi, H, psi)
variance = energy_squared - energy^2
println("\nEnergy variance = $variance")
t≈ttotal && break
psi = apply(gates, psi; cutoff)
normalize!(psi)
end
nlocal = expect(psi,"nb")
@show nlocal
return
end