Update examples and run in tests

ITensor · Feb 16, 2024 · d830f37 · d830f37
1 parent e3e2f5d
commit d830f37
Show file tree

Hide file tree

Showing 11 changed files with 264 additions and 182 deletions.
diff --git a/README.md b/README.md
@@ -32,6 +32,6 @@ pkg> add ITensorTDVP
 
 ## About
 
-This package is effectively a generalization of the DMRG code in [ITensors.jl](https://github.com/ITensor/ITensors.jl), using the MPS/MPO types from that package. It provides a general MPS "solver" interface which allows us to implement other MPS/MPO optimization/solver functionality like DMRG (`ITensorTDVP.dmrg` or `eigsolve`), TDVP (`tdvp` or `exponentiate`), linear solving (`linsolve`), DMRG-X (`dmrg_x`), etc. while sharing most of the code across those different functions. Therefore, it effectively supercedes the DMRG functionality in ITensors.jl (`dmrg`), and provides its own `ITensorTDVP.dmrg`/`eigsolve` function that is essentially the same as the `dmrg` function from ITensors.jl. This package is fairly stable and appropriate for general use. The primary missing feature is a lack of modern subspace expansion tools for methods like TDVP. However, 2-site TDVP or TEBD is often sufficient for performing subspace expansion (except when [it's not](https://arxiv.org/abs/2005.06104)).
+This package is effectively a generalization of the DMRG code in [ITensors.jl](https://github.com/ITensor/ITensors.jl), using the MPS/MPO types from that package. It provides a general MPS "solver" interface which allows us to implement other MPS/MPO optimization/solver functionality like DMRG (`ITensorTDVP.dmrg`), TDVP (`ITensorTDVP.tdvp`), linear solving (`ITensorTDVP.linsolve`/`KrylovKit.linsolve`), DMRG-X (`ITensorTDVP.dmrg_x`), etc. while sharing most of the code across those different functions. Therefore, it effectively supercedes the DMRG functionality in ITensors.jl (`dmrg`), and provides its own `ITensorTDVP.dmrg` function that is essentially the same as the `dmrg` function from ITensors.jl. This package is fairly stable and appropriate for general use. The primary missing feature is a lack of modern subspace expansion tools for methods like TDVP. However, 2-site TDVP or TEBD is often sufficient for performing subspace expansion (except when [it's not](https://arxiv.org/abs/2005.06104)).
 
 However, note that future developments, including modern subspace expansion tools, are being developed in our next-generation tensor network library [ITensorNetworks.jl](https://github.com/mtfishman/ITensorNetworks.jl). The goal of that package is to provide contraction, optimization, and evolution tools for general tensor networks, as well as methods like DMRG, TDVP, and linear solving for tree tensor networks, and the eventual goal is to replace this package which is limited to solvers for just MPS/MPO (linear/path graph) tensor networks. However, ITensorNetworks.jl is under heavy development and is _not_ meant for general usage at the moment, except for those who are brave enough to handle missing features and breaking interfaces. Basically, for the average user who wants stable and reliable code, if you need to use MPS-based TDVP or linear solving, you should use this package for the time being.
diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl
@@ -1,42 +1,47 @@
-using ITensors
-using ITensorTDVP
-
-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
+using ITensors: ITensors, MPO, OpSum, inner, randomMPS, siteinds
+using ITensorTDVP: ITensorTDVP, 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
-  return os
-end
 
-H = MPO(heisenberg(n), s)
-ψ = randomMPS(s, "↑"; linkdims=10)
+  H = MPO(heisenberg(n), s)
+  ψ = randomMPS(s, "↑"; linkdims=10)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
 
-@show inner(ψ', H, ψ) / inner(ψ, ψ)
+  ϕ = tdvp(
+    H,
+    -1.0,
+    ψ;
+    nsweeps=20,
+    reverse_step=false,
+    normalize=true,
+    maxdim=30,
+    cutoff=1e-10,
+    outputlevel=1,
+  )
 
-ϕ = tdvp(
-  H,
-  -1.0,
-  ψ;
-  nsweeps=20,
-  reverse_step=false,
-  normalize=true,
-  maxdim=30,
-  cutoff=1e-10,
-  outputlevel=1,
-)
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
 
-@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+  e2, ϕ2 = ITensors.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
 
-e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
+  @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)
 
-@show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)
+  ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
 
-ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
+  @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)
+  return nothing
+end
 
-@show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)
+main()
diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl
@@ -1,44 +1,48 @@
-using ITensors
-using ITensorTDVP
-using LinearAlgebra
-
-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
+using ITensors: MPO, MPS, OpSum, inner, siteinds
+using ITensorTDVP: dmrg_x
+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
-  return os
-end
 
-n = 10
-s = siteinds("S=1/2", n)
+  n = 10
+  s = siteinds("S=1/2", n)
 
-using Random
-Random.seed!(12)
+  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)
+  # 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)
+  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, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
+  )
 
-ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+  ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
+  @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+  @show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
+  return nothing
+end
 
-@show inner(ψ', H, ψ) / inner(ψ, ψ)
-@show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
-@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
-@show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
+main()
diff --git a/examples/03_tdvp_time_dependent.jl b/examples/03_tdvp_time_dependent.jl
@@ -149,6 +149,7 @@ function main()
 
   @show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full))
   @show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full))
+  return nothing
 end
 
 main()
diff --git a/examples/04_tdvp_observers.jl b/examples/04_tdvp_observers.jl
@@ -1,81 +1,86 @@
-using ITensors
-using ITensorTDVP
-using Observers
+using ITensors: MPO, MPS, OpSum, expect, siteinds
+using ITensorTDVP: tdvp
+using Observers: observer
 
-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
+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
-  return os
-end
 
-N = 10
-cutoff = 1e-12
-tau = 0.1
-ttotal = 1.0
+  N = 10
+  cutoff = 1e-12
+  tau = 0.1
+  ttotal = 1.0
 
-s = siteinds("S=1/2", N; conserve_qns=true)
-H = MPO(heisenberg(N), s)
+  s = siteinds("S=1/2", N; conserve_qns=true)
+  H = MPO(heisenberg(N), s)
 
-function step(; sweep, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return sweep
+  function step(; sweep, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return sweep
+    end
+    return nothing
   end
-  return nothing
-end
 
-function current_time(; current_time, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return current_time
+  function current_time(; current_time, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return current_time
+    end
+    return nothing
   end
-  return nothing
-end
 
-function measure_sz(; psi, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return expect(psi, "Sz"; sites=N ÷ 2)
+  function measure_sz(; psi, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return expect(psi, "Sz"; sites=N ÷ 2)
+    end
+    return nothing
   end
-  return nothing
-end
 
-function return_state(; psi, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return psi
+  function return_state(; psi, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return psi
+    end
+    return nothing
   end
-  return nothing
-end
 
-obs = observer(
-  "steps" => step, "times" => current_time, "psis" => return_state, "Sz" => measure_sz
-)
+  obs = observer(
+    "steps" => step, "times" => current_time, "psis" => return_state, "Sz" => measure_sz
+  )
 
-psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
-psi_f = tdvp(
-  H,
-  -im * ttotal,
-  psi;
-  time_step=-im * tau,
-  cutoff,
-  outputlevel=1,
-  normalize=false,
-  (observer!)=obs,
-)
+  psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi_f = tdvp(
+    H,
+    -im * ttotal,
+    psi;
+    time_step=-im * tau,
+    cutoff,
+    outputlevel=1,
+    normalize=false,
+    (observer!)=obs,
+  )
 
-steps = obs.steps
-times = obs.times
-psis = obs.psis
-Sz = obs.Sz
+  steps = obs.steps
+  times = obs.times
+  psis = obs.psis
+  Sz = obs.Sz
 
-println("\nResults")
-println("=======")
-for n in 1:length(steps)
-  print("step = ", steps[n])
-  print(", time = ", round(times[n]; digits=3))
-  print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi)); digits=3))
-  print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
-  print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
-  println()
+  println("\nResults")
+  println("=======")
+  for n in 1:length(steps)
+    print("step = ", steps[n])
+    print(", time = ", round(times[n]; digits=3))
+    print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi)); digits=3))
+    print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
+    print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
+    println()
+  end
+  return nothing
 end
+
+main()