ITensorTDVP.jl v0.4

Project.toml

@@ -1,9 +1,10 @@
 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]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -18,6 +19,7 @@ Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
 ITensorTDVPObserversExt = "Observers"
 
 [compat]
+Compat = "4"
 ITensors = "0.3.58, 0.4, 0.5, 0.6"
 KrylovKit = "0.6, 0.7"
 Observers = "0.2"

README.md

@@ -22,6 +22,95 @@ However, as noted above we now recommend installing and loading `ITensorMPS` ins
 
 ## News
 
+### ITensorTDVP.jl v0.4 Release Notes
+
+#### Breaking changes
+
+- When calling `tdvp(operator, t, init; kwargs...)`, `t` is now interpreted as the total evolution time, while in `ITensorTDVP.jl` v0.3 and below it was being interpreted as the time step. To upgrade, change code like:
+```julia
+tdvp(operator, -0.1im, init; nsweeps=10, kwargs...)
+```
+to:
+```julia
+tdvp(operator, -1.0im, init; nsweeps=10, kwargs...)
+```
+to evolve to time `1.0` over `10` steps with time steps of `0.1`. Also note that `ITensorTDVP.jl` v0.4 introduces `nsteps` as an alias for `nsweeps` in the `tdvp` function. `nsteps` is now the preferred syntax for specifying the number of time steps, and `nsweeps` may be deprecated in `tdvp` in the future. For example:
+```julia
+t = -1.0im
+nsteps = 10
+tdvp(operator, t, init; nsteps, kwargs...)
+```
+or specify the time steps instead of the number of steps:
+```julia
+t = -1.0im
+time_step = -0.1im
+tdvp(operator, t, init; time_step, kwargs...)
+```
+- In `tdvp`, a custom local updater/solver must now be passed as a keyword argument `updater`, as opposed to as the first argument which was the syntax in `ITensorTDVP.jl` v0.3 and below. So code like:
+```julia
+tdvp(custom_updater, operator, t, init; kwargs...)
+```
+must be changed to:
+```julia
+tdvp(operator, t, init; updater=custom_updater, kwargs...)
+```
+- The keyword argument `psi` that was being passed to observers in `tdvp`, `linsolve`, etc. which stored the current state has been renamed to `state`. Change code like:
+```julia
+measure_sz(; psi) = expect(psi, "Sz")
+obs = observer("Sz" => measure_sz)
+tdvp(operator, t, init; (observer!)=obs, kwargs...)
+```
+to:
+```julia
+measure_sz(; state) = expect(state, "Sz")
+obs = observer("Sz" => measure_sz)
+tdvp(operator, t, init; (observer!)=obs, kwargs...)
+```
+- Only the argument ordering `tdvp(operator, t, init; kwargs...)` is now supported. `tdvp(t, operator, init; kwargs...)` and `tdvp(operator, init, t; kwargs...)` have been removed.
+- In `tdvp`, the keyword argument `solver_backend` has been renamed to `updater_backend`. Change code like:
+```julia
+tdvp(operator, t, init; solver_backend="applyexp", kwargs...)
+```
+to:
+```julia
+tdvp(operator, t, init; updater_backend="applyexp", kwargs...)
+```
+- In `tdvp` and `ITensorTDVP.dmrg`, keyword arguments passed to the local solver should now be passed in a NamedTuple in the `updater_kwargs` keyword argument, such as `updater_kwargs=(; tol=1e-5, krylovdim=20)`, instead of as keyword arguments `solver_tol`, `solver_krylovdim`, etc. Change code like:
+```julia
+tdvp(operator, t, init; solver_tol=1e-5, krylovdim=20, kwargs...)
+ITensorTDVP.dmrg(operator, init; solver_tol=1e-5, solver_krylovdim=20, kwargs...)
+```
+to:
+```julia
+tdvp(operator, t, init; updater_kwargs=(; tol=1e-5, krylovdim=20), kwargs...)
+ITensorTDVP.dmrg(operator, init; updater_kwargs=(; tol=1e-5, krylovdim=20), kwargs...)
+```
+- In `linsolve`, the keyword argument `solver_kwargs` has been renamed to `updater_kwargs`.
+- In `ITensorTDVP.dmrg`, `dmrg_x`, and `linsolve`, the keyword argument `step_observer!` has been renamed to `sweep_observer!`. Either name is allowed in `tdvp` but `step_observer!` is preferred and the name `sweep_observer!` may be deprecated in `tdvp` in future versions.
+- Support for `ITensors.AbstractObserver`-based observers has been removed, use `Observers.observer` instead.
+- In `contract(operator::MPO, state::MPS; alg="fit", kwargs...)`, and `apply(operator::MPO, state::MPS; kwargs...)`, the keyword argument for specifying an initial guess for the result is now called `init` instead of `init_mps`. Additionally, in `contract`, `init` should have primed site indices, or more generally should have site indices which are those that are not shared by the input operator and state.
+- In custom local updaters/solvers, the keyword arguments `time_step`, `current_time`, and `outputlevel` are now being passed as a NamedTuple in a new keyword argument `internal_kwargs`. Change local updaters/solvers from:
+```julia
+function custom_updater(operator, init; time_step, current_time, outputlevel, kwargs...)
+  ## updater code
+end
+```
+to:
+```julia
+function custom_updater(operator, init; internal_kwargs, kwargs...)
+  (; time_step, current_time, outputlevel) = internal_kwargs
+  # Or if the updater only needs to use a subset of them:
+  # (; time_step) = internal_kwargs
+  ## updater code
+end
+```
+
+#### New features
+
+- `nsteps` is now an alias for the `nsweeps` keyword argument in `tdvp` and the preferred syntax for setting the number of time steps of TDVP.
+- `TimeDependentSum` now accepts coefficients and terms that are Tuples, along with the previous interface which accepted Vectors.
+- Custom local updaters/solvers can be passed in a keyword argument `updater` to `ITensorTDVP.dmrg`, `dmrg_x`, and `linsolve`, which is consistent with the new syntax for `tdvp`.
+
 ### ITensorTDVP.jl v0.3 Release Notes
 
 #### Breaking changes

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()
@@ -30,12 +31,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
-  )
-
-  e, ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+  e, ϕ = dmrg_x(H, ψ; nsweeps=10, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1)
 
   @show inner(ψ', H, ψ) / inner(ψ, ψ)
   @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2

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,23 @@
+using Compat: @compat
+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...)
+  @compat (; 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...)
+  @compat (; current_time, time_step) = (; current_time=zero(Bool), internal_kwargs...)
+  state, info = exponentiate(operator(current_time), time_step, init; kwargs...)
+  return state, (; info)
+end