From d830f371755b87245c72f78704c8d17d7301eb1e Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 16 Feb 2024 10:03:03 -0500 Subject: [PATCH] Update examples and run in tests --- README.md | 2 +- examples/01_tdvp.jl | 71 ++++++------ examples/02_dmrg-x.jl | 74 +++++++------ examples/03_tdvp_time_dependent.jl | 1 + examples/04_tdvp_observers.jl | 133 ++++++++++++----------- examples/05_tdvp_nonuniform_timesteps.jl | 84 +++++++------- examples/05_utils.jl | 44 ++++++-- src/tdvp_generic.jl | 1 + test/Project.toml | 1 + test/test_examples.jl | 19 ++++ test/test_tdvp.jl | 16 ++- 11 files changed, 264 insertions(+), 182 deletions(-) create mode 100644 test/test_examples.jl diff --git a/README.md b/README.md index 3837353..74c4866 100644 --- 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 index 0fad273..c386c0a 100644 --- 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 index 2d6ce60..6d9e670 100644 --- 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 index 124caf3..e012219 100644 --- 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 index 8ab3f8a..f2c056d 100644 --- 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() diff --git a/examples/05_tdvp_nonuniform_timesteps.jl b/examples/05_tdvp_nonuniform_timesteps.jl index bf6f67a..0a76206 100644 --- a/examples/05_tdvp_nonuniform_timesteps.jl +++ b/examples/05_tdvp_nonuniform_timesteps.jl @@ -1,46 +1,52 @@ -using ITensors -using ITensorTDVP +using ITensors: MPO, MPS, OpSum, ProjMPO, siteinds +using KrylovKit: exponentiate +using Observers: observer include("05_utils.jl") -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 -outputlevel = 1 -nsteps = 10 -time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps] - -obs = observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi) - -s = siteinds("S=1/2", N; conserve_qns=true) -H = MPO(heisenberg(N), s) - -psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn") -psi = tdvp_nonuniform_timesteps( - ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs -) - -times = obs.times -psis = obs.psis - -println("\nResults") -println("=======") -print("step = ", 0) -print(", time = ", zero(ComplexF64)) -print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3)) -println() -for n in 1:length(times) - print("step = ", n) - print(", time = ", round(times[n]; digits=3)) - print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3)) + N = 10 + cutoff = 1e-12 + outputlevel = 1 + nsteps = 10 + time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps] + + obs = observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi) + + s = siteinds("S=1/2", N; conserve_qns=true) + H = MPO(heisenberg(N), s) + + psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn") + psi = tdvp_nonuniform_timesteps( + ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs + ) + + times = obs.times + psis = obs.psis + + println("\nResults") + println("=======") + print("step = ", 0) + print(", time = ", zero(ComplexF64)) + print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3)) println() + for n in 1:length(times) + print("step = ", n) + print(", time = ", round(times[n]; digits=3)) + print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3)) + println() + end + return nothing end + +main() diff --git a/examples/05_utils.jl b/examples/05_utils.jl index b168b14..5d61efd 100644 --- a/examples/05_utils.jl +++ b/examples/05_utils.jl @@ -1,9 +1,7 @@ -using ITensors -using ITensorTDVP -using Observers -using Printf - -using ITensorTDVP: tdvp_solver, tdvp_step, process_sweeps, TDVPOrder +using ITensors: MPS, maxlinkdim +using ITensorTDVP: TDVPOrder, process_sweeps, tdvp_solver, tdvp_step, process_sweeps +using Observers: observer, update! +using Printf: @printf function tdvp_nonuniform_timesteps( solver, @@ -14,10 +12,15 @@ function tdvp_nonuniform_timesteps( time_start=0.0, order=2, (step_observer!)=observer(), + maxdim=ITensorTDVP.default_maxdim(), + mindim=ITensorTDVP.default_mindim(), + cutoff=ITensorTDVP.default_cutoff(), + noise=ITensorTDVP.default_noise(), + outputlevel=ITensorTDVP.default_outputlevel(), kwargs..., ) nsweeps = length(time_steps) - maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, kwargs...) + maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise) tdvp_order = TDVPOrder(order, Base.Forward) current_time = time_start for sw in 1:nsweeps @@ -55,6 +58,29 @@ function tdvp_nonuniform_timesteps( return psi end -function tdvp_nonuniform_timesteps(H, psi::MPS; kwargs...) - return tdvp_nonuniform_timesteps(tdvp_solver(; kwargs...), H, psi; kwargs...) +function tdvp_nonuniform_timesteps( + H, + psi::MPS; + ishermitian=ITensorTDVP.default_ishermitian(), + issymmetric=ITensorTDVP.default_issymmetric(), + solver_tol=ITensorTDVP.default_solver_tol(exponentiate), + solver_krylovdim=ITensorTDVP.default_solver_krylovdim(exponentiate), + solver_maxiter=ITensorTDVP.default_solver_maxiter(exponentiate), + solver_outputlevel=ITensorTDVP.default_solver_outputlevel(exponentiate), + kwargs..., +) + return tdvp_nonuniform_timesteps( + tdvp_solver( + exponentiate; + ishermitian, + issymmetric, + solver_tol, + solver_krylovdim, + solver_maxiter, + solver_outputlevel, + ), + H, + psi; + kwargs..., + ) end diff --git a/src/tdvp_generic.jl b/src/tdvp_generic.jl index 066e78f..6279ac8 100644 --- a/src/tdvp_generic.jl +++ b/src/tdvp_generic.jl @@ -7,6 +7,7 @@ using ITensors: checkdone!, disk, linkind, + maxlinkdim, permute function _tdvp_compute_nsweeps(t; time_step=default_time_step(t), nsweeps=default_nsweeps()) diff --git a/test/Project.toml b/test/Project.toml index 7f0787b..0cdcad6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,4 +5,5 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_examples.jl b/test/test_examples.jl new file mode 100644 index 0000000..9592425 --- /dev/null +++ b/test/test_examples.jl @@ -0,0 +1,19 @@ +@eval module $(gensym()) +using Suppressor: @suppress +using ITensorTDVP: ITensorTDVP +using Test: @testset +@testset "Run examples" begin + examples_files = [ + "01_tdvp.jl", + "02_dmrg-x.jl", + # "03_tdvp_time_dependent.jl", + "04_tdvp_observers.jl", + "05_tdvp_nonuniform_timesteps.jl", + ] + examples_path = joinpath(pkgdir(ITensorTDVP), "examples") + @testset "Running example file $f" for f in examples_files + println("Running example file $f") + @suppress include(joinpath(examples_path, f)) + end +end +end diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index 38bcb32..97c7146 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -1,6 +1,20 @@ @eval module $(gensym()) using ITensors: - ITensors, AbstractObserver, ITensor, MPO, MPS, OpSum, randomMPS, inner, op, siteinds + ITensors, + AbstractObserver, + ITensor, + MPO, + MPS, + OpSum, + apply, + dag, + expect, + inner, + noprime, + op, + prime, + randomMPS, + siteinds using ITensorTDVP: tdvp using KrylovKit: exponentiate using LinearAlgebra: norm