4th order Runge Kutta and TDVP

ITensor Julia Questions
1

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
2

Hi @SVGT,
It would really help us if you could provide a minimal working example.

3

Could you clarify what is meant by minimal working example? Would the code above not work?

Thanks,
SVGT

4

Please refer to this thread for a definition of MWE Please read: make it easier to help you - Meta Discussion - Julia Programming Language, in particular the 5-th point.

It is difficult for someone to fully debug the code of another person and it could take far more time then required.

That said, the first thing I’d advise you to do is to rewrite your code following the style guide Style Guide · The Julia Language. This itself can a lot of time avoid errors, and in case there are, it will be easier to debug it and more understandable for other people.

Other things related more directly to your question: Unfortunately there does not seem to be explicit errors, or at least none that I’m able to notice, so here there are some things that can help in debugging it:

  • Does your code work for non fermionic hamiltonians?
  • Does one of the 2 methods report the correct results in an analytical case?