TDVP vs TEBD results different on a small Heisenberg chain

[ZhenZhao]

Hi there,

Many thanks for the ITensorTDVP.jl v0.4! I was testing it after its recent release, using an AFM Heisenberg chain with 4 sites, and time-dependent B field on the first site,

H(t)=\sum_{i=1}^3\mathbf{S}_i\cdot\mathbf{S}_{i+1}-B(t)S_1^z.

At t=0, I set B=3, and use dmrg to get a ground state MPS \psi(0). For t>0, B=0, and I use a)TDVP and b) TEBD to evolve \psi(0) with time_step=0.05 until t=1. For each time step, I measure \langle S_1^z\rangle. My TEBD results agree with ED benchmark, but TDVP results are pretty different.

I guess I might input arguments in the function tdvp in a wrong way. However, if I use tdvp and H(t=0) to evolve \psi(0), \psi(0) stays unchanged, which seems to be correct. Does this mean that, since \psi(0) has small maxlinkdim, I should evolve it with TEBD for several steps then use TDVP?

Below I put my codes. Thanks again!

Best,
Zhen

using ITensors, ITensorMPS
using Observers: observer

function spMPS()
# heisenberg chain
  N = 4
  spsites = siteinds("S=1/2",N;conserve_qns=true,conserve_sz=false)
  cutoff = 1E-8
  B = 3.0 # only on site 1
   
  # H0 with B
  os = OpSum()
  os -= B, "Sz", 1

  for i in 1:(N-1)
    os += 1.0/2, "S+", i, "S-", i+1
    os += 1.0/2, "S-", i, "S+", i+1
    os += 1.0, "Sz", i, "Sz", i+1
  end
  H = MPO(os,spsites)
  # H1 with B=0
  os1 = OpSum()
  for i in 1:(N-1)
    os1 += 1.0/2, "S+", i, "S-", i+1
    os1 += 1.0/2, "S-", i, "S+", i+1
    os1 += 1.0, "Sz", i, "Sz", i+1
  end
  H1 = MPO(os1,spsites)
  
  # dmrg
  state = [isodd(i) ? "Up" : "Dn" for i in 1:N] 
  psi0 = randomMPS(spsites, state)
  obs = DMRGObserver(; energy_tol=1e-8)
  
  nsweeps = 50
  maxdim = [32, 50, 100, 200, 400]#, 800]
  noise = [1E-3, 1E-4, 1E-5, 1E-6, 1E-7, 1E-8, 1E-10, 0.0]
  
  energy, psi = dmrg(H, psi0; 
               nsweeps, maxdim, cutoff,noise,observer=obs,eigsolve_krylovdim=8)
  @show energy
  # tdvp, with observer to measure <S1z(t)>
  step(; sweep) = sweep
  current_time(; current_time) = current_time
  return_state(; state) = state
  measure_sz(; state) = expect(state, "Sz"; sites=1)
  tdvpobs = observer(
    "steps" => step, "times" => current_time, 
    "states" => return_state, "Sz" => measure_sz
  )
  init = copy(psi)
  state = tdvp(
    H1, 1.0, init; time_step=0.05, cutoff=1e-12, (step_observer!)=tdvpobs,
    reverse_step=true, normalize=true,maxdim=100, 
    nsite = 2, updater_backend="exponentiate",outputlevel=1
  )
  
  println("\nResults")
  println("=======")
  for n in 1:length(tdvpobs.steps)
    print("step = ", tdvpobs.steps[n])
    print(", time = ", round(tdvpobs.times[n]; digits=3))
    print(", norm = ", round(abs(inner(tdvpobs.states[n], tdvpobs.states[n])); digits=3))
    print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(tdvpobs.states[n], init)); digits=3))
    print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(tdvpobs.states[n], state)); digits=3))
    print(", ⟨S1z⟩ = ", round(tdvpobs.Sz[n]; digits=3))
    println()
  end
return

spMPS()

[yangxusolidstate]

Hey Zhen, I think in tdvp the time \tau is translated to e^{ \tau*H} so you should use -im*1.0 for your real time evolution.

[mtfishman]

Also the time step should be -im * 0.05, again if you are doing real time evolution.

[ZhenZhao]

Thanks for the answer.

\Zhen