From 3cc10c1405dbd7091f889765ad139b325047e8b9 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Fri, 24 May 2024 15:44:27 -0400 Subject: [PATCH] Add solver examples and export contract (#20) * Add solver examples and export contract * Bump to v0.2.3 --- Project.toml | 2 +- examples/01_tdvp.jl | 41 ++++++ examples/02_dmrg-x.jl | 42 ++++++ examples/03_models.jl | 20 +++ examples/03_tdvp_time_dependent.jl | 153 ++++++++++++++++++++++ examples/03_updaters.jl | 23 ++++ examples/04_tdvp_observers.jl | 45 +++++++ examples/Project.toml | 11 ++ src/ITensorMPS.jl | 1 + test/Project.toml | 7 + test/runtests.jl | 48 ++----- test/test_examples.jl | 15 +++ test/test_exports.jl | 46 +++++++ test/utils/TestITensorMPSExportedNames.jl | 1 + 14 files changed, 414 insertions(+), 41 deletions(-) create mode 100644 examples/01_tdvp.jl create mode 100644 examples/02_dmrg-x.jl create mode 100644 examples/03_models.jl create mode 100644 examples/03_tdvp_time_dependent.jl create mode 100644 examples/03_updaters.jl create mode 100644 examples/04_tdvp_observers.jl create mode 100644 examples/Project.toml create mode 100644 test/test_examples.jl create mode 100644 test/test_exports.jl diff --git a/Project.toml b/Project.toml index ae48ba2..71c8684 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorMPS" uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" authors = ["Matthew Fishman ", "Miles Stoudenmire "] -version = "0.2.2" +version = "0.2.3" [deps] ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342" diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl new file mode 100644 index 0000000..d4be3bc --- /dev/null +++ b/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() diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl new file mode 100644 index 0000000..4305136 --- /dev/null +++ b/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() diff --git a/examples/03_models.jl b/examples/03_models.jl new file mode 100644 index 0000000..51f5aee --- /dev/null +++ b/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 diff --git a/examples/03_tdvp_time_dependent.jl b/examples/03_tdvp_time_dependent.jl new file mode 100644 index 0000000..85b79a8 --- /dev/null +++ b/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() diff --git a/examples/03_updaters.jl b/examples/03_updaters.jl new file mode 100644 index 0000000..ee617b7 --- /dev/null +++ b/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 diff --git a/examples/04_tdvp_observers.jl b/examples/04_tdvp_observers.jl new file mode 100644 index 0000000..2d6b807 --- /dev/null +++ b/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() diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 0000000..558a964 --- /dev/null +++ b/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" diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 6932f32..3618787 100644 --- a/src/ITensorMPS.jl +++ b/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, diff --git a/test/Project.toml b/test/Project.toml index 72cc60a..69fd6a2 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/runtests.jl b/test/runtests.jl index 1424d2f..44e961e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,46 +1,14 @@ @eval module $(gensym()) +using Test: @testset using ITensorMPS: ITensorMPS -using ITensorTDVP: ITensorTDVP -using ITensors: ITensors -include("utils/TestITensorMPSExportedNames.jl") -using Test: @test, @test_broken, @testset +test_path = joinpath(pkgdir(ITensorMPS), "test") +test_files = filter( + file -> startswith(file, "test_") && endswith(file, ".jl"), readdir(test_path) +) @testset "ITensorMPS.jl" begin - @testset "exports" begin - @test issetequal( - names(ITensorMPS), - [ - [:ITensorMPS] - # ITensorTDVP reexports - [:TimeDependentSum, :dmrg_x, :expand, :linsolve, :tdvp, :to_vec] - # ITensors and ITensors.ITensorMPS reexports - TestITensorMPSExportedNames.ITENSORMPS_EXPORTED_NAMES - ], - ) - end - @testset "Aliases" begin - @test ITensorMPS.Experimental.dmrg === ITensorTDVP.dmrg - @test ITensorMPS.dmrg === ITensors.ITensorMPS.dmrg - end - @testset "Not exported" begin - @test ITensorMPS.sortmergeterms === ITensors.ITensorMPS.sortmergeterms - # Should we fix this in ITensors.jl by adding: - # ```julia - # using .ITensorMPS: sortmergeterms - # ``` - # ? - @test_broken ITensorMPS.sortmergeterms === ITensors.sortmergeterms - for f in [ - :AbstractProjMPO, - :AbstractMPS, - :ProjMPS, - :makeL!, - :makeR!, - :set_terms, - :sortmergeterms, - :terms, - ] - @test getfield(ITensorMPS, f) === getfield(ITensors.ITensorMPS, f) - end + @testset "$filename" for filename in test_files + println("Running $filename") + @time include(joinpath(test_path, filename)) end end end diff --git a/test/test_examples.jl b/test/test_examples.jl new file mode 100644 index 0000000..3e0f481 --- /dev/null +++ b/test/test_examples.jl @@ -0,0 +1,15 @@ +@eval module $(gensym()) +using Suppressor: @suppress +using ITensorMPS: ITensorMPS +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" + ] + examples_path = joinpath(pkgdir(ITensorMPS), "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_exports.jl b/test/test_exports.jl new file mode 100644 index 0000000..7515de7 --- /dev/null +++ b/test/test_exports.jl @@ -0,0 +1,46 @@ +@eval module $(gensym()) +using ITensorMPS: ITensorMPS +using ITensorTDVP: ITensorTDVP +using ITensors: ITensors +include("utils/TestITensorMPSExportedNames.jl") +using Test: @test, @test_broken, @testset +@testset "Exports and aliases" begin + @testset "Exports" begin + @test issetequal( + names(ITensorMPS), + [ + [:ITensorMPS] + # ITensorTDVP reexports + [:TimeDependentSum, :dmrg_x, :expand, :linsolve, :tdvp, :to_vec] + # ITensors and ITensors.ITensorMPS reexports + TestITensorMPSExportedNames.ITENSORMPS_EXPORTED_NAMES + ], + ) + end + @testset "Aliases" begin + @test ITensorMPS.Experimental.dmrg === ITensorTDVP.dmrg + @test ITensorMPS.dmrg === ITensors.ITensorMPS.dmrg + end + @testset "Not exported" begin + @test ITensorMPS.sortmergeterms === ITensors.ITensorMPS.sortmergeterms + # Should we fix this in ITensors.jl by adding: + # ```julia + # using .ITensorMPS: sortmergeterms + # ``` + # ? + @test_broken ITensorMPS.sortmergeterms === ITensors.sortmergeterms + for f in [ + :AbstractProjMPO, + :AbstractMPS, + :ProjMPS, + :makeL!, + :makeR!, + :set_terms, + :sortmergeterms, + :terms, + ] + @test getfield(ITensorMPS, f) === getfield(ITensors.ITensorMPS, f) + end + end +end +end diff --git a/test/utils/TestITensorMPSExportedNames.jl b/test/utils/TestITensorMPSExportedNames.jl index 52deb77..7993079 100644 --- a/test/utils/TestITensorMPSExportedNames.jl +++ b/test/utils/TestITensorMPSExportedNames.jl @@ -51,6 +51,7 @@ const ITENSORMPS_EXPORTED_NAMES = [ :coefficient, :common_siteind, :common_siteinds, + :contract, :convert_leaf_eltype, :correlation_matrix, :cutoff,