Hello! I’ve been trying to implement a 4th order runge-kutta integrator for dynamics and compare with the ITensorTDVP package. Here’s the code for a 6-site system, I’m seeing small discrepancies with certain observables (electron density on site 3, for example) and am wondering if there is something in the code setup that I am not understanding properly or could improve.
Is my code for Runge-Kutta setup properly?
Thanks
using ITensors, ITensorMPS
let
N = 6
NL = N/3
Ndot = N/3
LastDot = NL + Ndot
t_ham = 1
U = 0
U_end = 1
nprint = 20
Bfield = 50
sites = siteinds("Electron",N; conserve_qns=false)
state = ["Up" for i in 1:N]
os = OpSum()
for l=NL+1:LastDot
l = Int(l)
os += -Bfield, "Sz", l
end
for j=1:N-1
os += -t_ham, "Cdagup", j, "Cup", j+1
os += -t_ham, "Cdagup", j+1, "Cup", j
os += -t_ham, "Cdagdn", j+1, "Cdn", j
os += -t_ham, "Cdagdn", j, "Cdn", j+1
end
H = MPO(os,sites)
nsweeps = 6 # number of sweeps is 6
maxdim = [10,20,100,100,200, 500] # gradually increase states kept
cutoff = [1E-10] # desired truncation error
ϕ = random_mps(sites)
energy,ϕ = dmrg(H,ϕ;nsweeps,maxdim,cutoff)
println("Static Energy: ", energy)
# DYNAMICS ----------------------------------------------------------
os = OpSum()
for l=NL+1:LastDot
l = Int(l)
os += Bfield, "Sz", l
end
for j=1:N-1
os += -t_ham, "Cdagup", j, "Cup", j+1
os += -t_ham, "Cdagup", j+1, "Cup", j
os += -t_ham, "Cdagdn", j+1, "Cdn", j
os += -t_ham, "Cdagdn", j, "Cdn", j+1
end
H_dyn = MPO(os, sites)
tau_real = 0.001
tau = -0.001 * 1im
total_steps = 5001
currstep = 0
ntot = expect(ϕ,"Ntot", sites=1:N)
println(" ")
println("edens at t=0: ", ntot)
ntot = expect(ϕ,"Sz", sites=1:N)
println(" ")
println("spinz at t=0: ", ntot)
ntot = expect(ϕ,"S-", sites=1:N)
# PRESERVE PRE-DYNAMICS PSI
ϕ1 = ϕ
for i in 1:total_steps
k1 = tau * contract(NDTensors.Algorithm(:densitymatrix), H_dyn, ϕ; cutoff = 1e-10, maxdim = 500)
noprime!(k1)
temp1 = add(ϕ,0.5*k1; cutoff=1e-10, maxdim=1000)
k2 = tau * contract(NDTensors.Algorithm(:densitymatrix), H_dyn, temp1; cutoff = 1e-10, maxdim = 500)
noprime!(k2)
temp2 = add(ϕ,0.5*k2; cutoff=1e-10, maxdim=1000)
k3 = tau * contract(NDTensors.Algorithm(:densitymatrix), H_dyn, temp2; cutoff = 1e-10, maxdim = 500)
noprime!(k3)
temp3 = add(ϕ,k3; cutoff=1e-10, maxdim=1000)
k4 = tau * contract(NDTensors.Algorithm(:densitymatrix), H_dyn, temp3; cutoff = 1e-10, maxdim = 500)
noprime!(k4)
t1 = (1/6) * k1
t2 = (1/3) * k2
t3 = (1/3) * k3
t4 = (1/6) * k4
ϕ = add(t1, t2, t3, t4, ϕ; alg = "densitymatrix", cutoff=1e-10, maxdim=1000)
noprime!(ϕ)
normalize!(ϕ)
currstep += 1
if currstep % nprint == 0
println("Step: ", currstep * tau_real)
flush(stdout)
ntot = expect(ϕ,"Ntot", sites=1:N)
println(" ")
println("edens: ", ntot)
flush(stdout)
ntot = expect(ϕ,"Sz", sites=1:N)
println(" ")
println("Spin Z: ", ntot)
flush(stdout)
end
end
println("END OF RK4 ############################### BEGIN TDVP")
flush(stdout)
# REASSIGN PSI TO INITIAL CONDITIONS
ϕ = ϕ1
nprint = 2
tau_real = 0.01
tau = -0.01 * 1im
total_steps = 500
currstep = 0
for i in 1:total_steps
ϕ = tdvp(
H_dyn,
tau,
ϕ,
time_step = tau,
normalize = true,
maxdim = 1000,
cutoff = 1e-7,
outputlevel = 1,
)
currstep += 1
if currstep % nprint == 0
println("Step: ", currstep * tau_real)
flush(stdout)
ntot = expect(ϕ,"Ntot", sites=1:N)
println(" ")
println("edens: ", ntot)
ntot = expect(ϕ,"Sz", sites=1:N)
println(" ")
println("Spin Z: ", ntot)
flush(stdout)
end
end
flush(stdout)
end