Add solver examples and export contract

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensorMPS"
 uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"

examples/01_tdvp.jl

@@ -0,0 +1,41 @@
+using ITensorMPS: MPO, OpSum, dmrg, inner, random_mps, siteinds, tdvp
+
+function main()
+  n = 10
+  s = siteinds("S=1/2", n)
+
+  function heisenberg(n)
+    os = OpSum()
+    for j in 1:(n - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    return os
+  end
+
+  H = MPO(heisenberg(n), s)
+  ψ = random_mps(s, "↑"; linkdims=10)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
+
+  ϕ = tdvp(
+    H,
+    -20.0,
+    ψ;
+    time_step=-1.0,
+    maxdim=30,
+    cutoff=1e-10,
+    normalize=true,
+    reverse_step=false,
+    outputlevel=1,
+  )
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+
+  e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
+  @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2), e2
+
+  return nothing
+end
+
+main()

examples/02_dmrg-x.jl

@@ -0,0 +1,42 @@
+using ITensorMPS: MPO, MPS, OpSum, dmrg_x, inner, siteinds
+using Random: Random
+
+function main()
+  function heisenberg(n; h=zeros(n))
+    os = OpSum()
+    for j in 1:(n - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    for j in 1:n
+      if h[j] ≠ 0
+        os -= h[j], "Sz", j
+      end
+    end
+    return os
+  end
+
+  n = 10
+  s = siteinds("S=1/2", n)
+
+  Random.seed!(12)
+
+  # MBL when W > 3.5-4
+  W = 12
+  # Random fields h ∈ [-W, W]
+  h = W * (2 * rand(n) .- 1)
+  H = MPO(heisenberg(n; h), s)
+
+  initstate = rand(["↑", "↓"], n)
+  ψ = MPS(s, initstate)
+  e, ϕ = dmrg_x(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, normalize=true, outputlevel=1)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
+  @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ), e
+  @show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
+  return nothing
+end
+
+main()

examples/03_models.jl

@@ -0,0 +1,20 @@
+using ITensorMPS: OpSum
+
+function heisenberg(n; J=1.0, J2=0.0)
+  ℋ = OpSum()
+  if !iszero(J)
+    for j in 1:(n - 1)
+      ℋ += J / 2, "S+", j, "S-", j + 1
+      ℋ += J / 2, "S-", j, "S+", j + 1
+      ℋ += J, "Sz", j, "Sz", j + 1
+    end
+  end
+  if !iszero(J2)
+    for j in 1:(n - 2)
+      ℋ += J2 / 2, "S+", j, "S-", j + 2
+      ℋ += J2 / 2, "S-", j, "S+", j + 2
+      ℋ += J2, "Sz", j, "Sz", j + 2
+    end
+  end
+  return ℋ
+end

examples/03_tdvp_time_dependent.jl

@@ -0,0 +1,153 @@
+using ITensors: @disable_warn_order
+using ITensorMPS: MPO, MPS, contract, inner, random_mps, siteinds, tdvp
+using LinearAlgebra: norm
+using Random: Random
+
+include("03_models.jl")
+include("03_updaters.jl")
+
+function main()
+  Random.seed!(1234)
+
+  # Time dependent Hamiltonian is:
+  # H(t) = H₁(t) + H₂(t) + …
+  #      = f₁(t) H₁(0) + f₂(t) H₂(0) + …
+  #      = cos(ω₁t) H₁(0) + cos(ω₂t) H₂(0) + …
+
+  # Number of sites
+  n = 6
+
+  # 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
+  # updater.
+  outputlevel = 3
+
+  # Frequency of time dependent terms
+  ω₁ = 0.1
+  ω₂ = 0.2
+
+  # Nearest and next-nearest neighbor
+  # Heisenberg couplings.
+  J₁ = 1.0
+  J₂ = 1.0
+
+  time_step = 0.1
+  time_stop = 1.0
+
+  # nsite-update TDVP
+  nsite = 2
+
+  # Starting state bond/link dimension.
+  # A product state starting state can
+  # cause issues for TDVP without
+  # subspace expansion.
+  start_linkdim = 4
+
+  # TDVP truncation parameters
+  maxdim = 100
+  cutoff = 1e-8
+
+  tol = 1e-15
+
+  @show n
+  @show ω₁, ω₂
+  @show J₁, J₂
+  @show maxdim, cutoff, nsite
+  @show start_linkdim
+  @show time_step, time_stop
+
+  ω⃗ = (ω₁, ω₂)
+  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⃗₀ = map(ℋ₀ -> MPO(ℋ₀, s), ℋ⃗₀)
+
+  # Initial state, ψ₀ = ψ(0)
+  # Initialize as complex since that is what OrdinaryDiffEq.jl/DifferentialEquations.jl
+  # expects.
+  ψ₀ = complex.(random_mps(s, j -> isodd(j) ? "↑" : "↓"; linkdims=start_linkdim))
+
+  @show norm(ψ₀)
+
+  println()
+  println("#"^100)
+  println("Running TDVP with ODE updater")
+  println("#"^100)
+  println()
+
+  ψₜ_ode = tdvp(
+    -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 updater")
+  println()
+
+  println()
+  println("#"^100)
+  println("Running TDVP with Krylov updater")
+  println("#"^100)
+  println()
+
+  ψₜ_krylov = tdvp(
+    -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 updater")
+  println()
+
+  println()
+  println("#"^100)
+  println("Running full state evolution with ODE updater")
+  println("#"^100)
+  println()
+
+  @disable_warn_order begin
+    ψₜ_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 updater")
+  println()
+
+  @show norm(ψₜ_ode)
+  @show norm(ψₜ_krylov)
+  @show norm(ψₜ_full)
+
+  @show 1 - abs(inner(contract(ψₜ_ode), ψₜ_full))
+  @show 1 - abs(inner(contract(ψₜ_krylov), ψₜ_full))
+  return nothing
+end
+
+main()

examples/03_updaters.jl

@@ -0,0 +1,23 @@
+using Compat: @compat
+using ITensors: ITensor, array, inds, itensor
+using ITensorMPS: 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

examples/04_tdvp_observers.jl

@@ -0,0 +1,45 @@
+using ITensorMPS: MPO, MPS, OpSum, expect, inner, siteinds, tdvp
+using Observers: observer
+
+function main()
+  function heisenberg(N)
+    os = OpSum()
+    for j in 1:(N - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    return os
+  end
+
+  N = 10
+  s = siteinds("S=1/2", N; conserve_qns=true)
+  H = MPO(heisenberg(N), s)
+
+  step(; sweep) = sweep
+  current_time(; current_time) = current_time
+  return_state(; state) = state
+  measure_sz(; state) = expect(state, "Sz"; sites=length(state) ÷ 2)
+  obs = observer(
+    "steps" => step, "times" => current_time, "states" => return_state, "sz" => measure_sz
+  )
+
+  init = MPS(s, n -> isodd(n) ? "Up" : "Dn")
+  state = tdvp(
+    H, -1.0im, init; time_step=-0.1im, cutoff=1e-12, (step_observer!)=obs, outputlevel=1
+  )
+
+  println("\nResults")
+  println("=======")
+  for n in 1:length(obs.steps)
+    print("step = ", obs.steps[n])
+    print(", time = ", round(obs.times[n]; digits=3))
+    print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(obs.states[n], init)); digits=3))
+    print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(obs.states[n], state)); digits=3))
+    print(", ⟨Sᶻ⟩ = ", round(obs.sz[n]; digits=3))
+    println()
+  end
+  return nothing
+end
+
+main()

examples/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[compat]
+ITensors = "0.6.7"

src/ITensorMPS.jl

@@ -10,6 +10,7 @@ using .Experimental: Experimental
 include("Deprecated.jl")
 using .Deprecated: Deprecated, dmrg
 export dmrg
+@reexport using ITensors: contract
 @reexport using ITensors.ITensorMPS:
   @OpName_str,
   @SiteType_str,

test/Project.toml

@@ -1,5 +1,12 @@
 [deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"