ITensorTDVP.jl v0.4

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensorTDVP"
 uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org> and contributors"]
-version = "0.3.1"
+version = "0.4.0"
 
 [deps]
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"

examples/01_tdvp.jl

@@ -1,4 +1,5 @@
-using ITensorMPS: MPO, OpSum, dmrg, inner, randomMPS, siteinds, tdvp
+using ITensors: MPO, OpSum, dmrg, inner, randomMPS, siteinds
+using ITensorTDVP: tdvp
 
 function main()
   n = 10
@@ -21,9 +22,9 @@ function main()
 
   ϕ = tdvp(
     H,
-    -1.0,
+    -20.0,
     ψ;
-    nsweeps=20,
+    time_step=-1.0,
     reverse_step=false,
     normalize=true,
     maxdim=30,

examples/02_dmrg-x.jl

@@ -1,4 +1,5 @@
-using ITensorMPS: MPO, MPS, OpSum, dmrg_x, inner, siteinds
+using ITensors: MPO, MPS, OpSum, inner, siteinds
+using ITensorTDVP: dmrg_x
 using Random: Random
 
 function main()
@@ -31,9 +32,7 @@ function main()
   initstate = rand(["↑", "↓"], n)
   ψ = MPS(s, initstate)
 
-  dmrg_x_kwargs = (
-    nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
-  )
+  dmrg_x_kwargs = (nsweeps=10, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1)
 
   e, ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
 

examples/03_models.jl

@@ -1,4 +1,4 @@
-using ITensorMPS: OpSum
+using ITensors: OpSum
 
 function heisenberg(n; J=1.0, J2=0.0)
   ℋ = OpSum()

examples/03_tdvp_time_dependent.jl

@@ -1,10 +1,11 @@
-using ITensors: @disable_warn_order
-using ITensorMPS: MPO, MPS, inner, randomMPS, siteinds, tdvp
+using ITensors: @disable_warn_order, contract
+using ITensors: MPO, MPS, inner, randomMPS, siteinds
+using ITensorTDVP: tdvp
 using LinearAlgebra: norm
 using Random: Random
 
 include("03_models.jl")
-include("03_solvers.jl")
+include("03_updaters.jl")
 
 function main()
   Random.seed!(1234)
@@ -20,7 +21,7 @@ function main()
   # How much information to output from TDVP
   # Set to 2 to get information about each bond/site
   # evolution, and 3 to get information about the
-  # solver.
+  # updater.
   outputlevel = 3
 
   # Frequency of time dependent terms
@@ -50,34 +51,24 @@ function main()
 
   tol = 1e-15
 
-  # ODE solver parameters
-  ode_alg = Tsit5()
-  ode_kwargs = (; reltol=tol, abstol=tol)
-
-  # Krylov solver parameters
-  krylov_kwargs = (; tol=tol, eager=true)
-
   @show n
   @show ω₁, ω₂
   @show J₁, J₂
   @show maxdim, cutoff, nsite
   @show start_linkdim
   @show time_step, time_stop
-  @show ode_alg
-  @show ode_kwargs
-  @show krylov_kwargs
 
-  ω⃗ = [ω₁, ω₂]
-  f⃗ = [t -> cos(ω * t) for ω in ω⃗]
+  ω⃗ = (ω₁, ω₂)
+  f⃗ = map(ω -> (t -> cos(ω * t)), ω⃗)
 
   # H₀ = H(0) = H₁(0) + H₂(0) + …
   ℋ₁₀ = heisenberg(n; J=J₁, J2=0.0)
   ℋ₂₀ = heisenberg(n; J=0.0, J2=J₂)
-  ℋ⃗₀ = [ℋ₁₀, ℋ₂₀]
+  ℋ⃗₀ = (ℋ₁₀, ℋ₂₀)
 
   s = siteinds("S=1/2", n)
 
-  H⃗₀ = [MPO(ℋ₀, s) for ℋ₀ in ℋ⃗₀]
+  H⃗₀ = map(ℋ₀ -> MPO(ℋ₀, s), ℋ⃗₀)
 
   # Initial state, ψ₀ = ψ(0)
   # Initialize as complex since that is what OrdinaryDiffEq.jl/DifferentialEquations.jl
@@ -88,60 +79,75 @@ function main()
 
   println()
   println("#"^100)
-  println("Running TDVP with ODE solver")
+  println("Running TDVP with ODE updater")
   println("#"^100)
   println()
 
-  function ode_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
-    return ode_solver(f⃗, H⃗₀, time_step, ψ₀; solver_alg=ode_alg, ode_kwargs..., kwargs...)
-  end
-
   ψₜ_ode = tdvp(
-    ode_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel
+    -im * TimeDependentSum(f⃗, H⃗₀),
+    time_stop,
+    ψ₀;
+    updater=ode_updater,
+    updater_kwargs=(; reltol=tol, abstol=tol),
+    time_step,
+    maxdim,
+    cutoff,
+    nsite,
+    outputlevel,
   )
 
   println()
-  println("Finished running TDVP with ODE solver")
+  println("Finished running TDVP with ODE updater")
   println()
 
   println()
   println("#"^100)
-  println("Running TDVP with Krylov solver")
+  println("Running TDVP with Krylov updater")
   println("#"^100)
   println()
 
-  function krylov_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
-    return krylov_solver(f⃗, H⃗₀, time_step, ψ₀; krylov_kwargs..., kwargs...)
-  end
-
   ψₜ_krylov = tdvp(
-    krylov_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel
+    -im * TimeDependentSum(f⃗, H⃗₀),
+    time_stop,
+    ψ₀;
+    updater=krylov_updater,
+    updater_kwargs=(; tol, eager=true),
+    time_step,
+    cutoff,
+    nsite,
+    outputlevel,
   )
 
   println()
-  println("Finished running TDVP with Krylov solver")
+  println("Finished running TDVP with Krylov updater")
   println()
 
   println()
   println("#"^100)
-  println("Running full state evolution with ODE solver")
+  println("Running full state evolution with ODE updater")
   println("#"^100)
   println()
 
   @disable_warn_order begin
-    ψₜ_full, _ = ode_solver(f⃗, prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel)
+    ψₜ_full, _ = ode_updater(
+      -im * TimeDependentSum(f⃗, contract.(H⃗₀)),
+      contract(ψ₀);
+      internal_kwargs=(; time_step=time_stop, outputlevel),
+      reltol=tol,
+      abstol=tol,
+    )
   end
 
   println()
-  println("Finished full state evolution with ODE solver")
+  println("Finished full state evolution with ODE updater")
   println()
 
   @show norm(ψₜ_ode)
   @show norm(ψₜ_krylov)
   @show norm(ψₜ_full)
 
-  @show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full))
-  @show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full))
+  @show 1 - abs(inner(contract(ψₜ_ode), ψₜ_full))
+  @show 1 - abs(inner(contract(ψₜ_krylov), ψₜ_full))
   return nothing
 end
 

examples/03_updaters.jl

@@ -0,0 +1,22 @@
+using ITensors: ITensor, array, inds, itensor
+using ITensorTDVP: TimeDependentSum, to_vec
+using KrylovKit: exponentiate
+using OrdinaryDiffEq: ODEProblem, Tsit5, solve
+
+function ode_updater(operator, init; internal_kwargs, alg=Tsit5(), kwargs...)
+  (; current_time, time_step) = (; current_time=zero(Bool), internal_kwargs...)
+  time_span = typeof(time_step).((current_time, current_time + time_step))
+  init_vec, to_itensor = to_vec(init)
+  f(init::ITensor, p, t) = operator(t)(init)
+  f(init_vec::Vector, p, t) = to_vec(f(to_itensor(init_vec), p, t))[1]
+  prob = ODEProblem(f, init_vec, time_span)
+  sol = solve(prob, alg; kwargs...)
+  state_vec = sol.u[end]
+  return to_itensor(state_vec), (;)
+end
+
+function krylov_updater(operator, init; internal_kwargs, kwargs...)
+  (; current_time, time_step) = (; current_time=zero(Bool), internal_kwargs...)
+  state, info = exponentiate(operator(current_time), time_step, init; kwargs...)
+  return state, (; info)
+end

examples/04_tdvp_observers.jl

@@ -1,4 +1,5 @@
-using ITensorMPS: MPO, MPS, OpSum, expect, inner, siteinds, tdvp
+using ITensors: MPO, MPS, OpSum, expect, inner, siteinds
+using ITensorTDVP: tdvp
 using Observers: observer
 
 function main()

src/ITensorTDVP.jl

@@ -1,20 +1,21 @@
 module ITensorTDVP
 export TimeDependentSum, dmrg_x, linsolve, tdvp, to_vec
+include("ITensorsExtensions.jl")
+using .ITensorsExtensions: to_vec
 include("applyexp.jl")
 include("defaults.jl")
 include("update_observer.jl")
-include("solver_utils.jl")
+include("timedependentsum.jl")
 include("tdvporder.jl")
-include("tdvpinfo.jl")
 include("sweep_update.jl")
 include("alternating_update.jl")
 include("tdvp.jl")
 include("dmrg.jl")
 include("dmrg_x.jl")
-include("projmpo_apply.jl")
-include("contract_mpo_mps.jl")
-include("projmps2.jl")
-include("projmpo_mps2.jl")
+include("reducedcontractproblem.jl")
+include("contract.jl")
+include("reducedconstantterm.jl")
+include("reducedlinearproblem.jl")
 include("linsolve.jl")
 using PackageExtensionCompat: @require_extensions
 function __init__()

src/ITensorsExtensions.jl

@@ -0,0 +1,9 @@
+module ITensorsExtensions
+using ITensors: ITensor, array, inds, itensor
+function to_vec(x::ITensor)
+  function to_itensor(x_vec)
+    return itensor(x_vec, inds(x))
+  end
+  return vec(array(x)), to_itensor
+end
+end