diff --git a/examples/treetensornetworks/solvers/01_tdvp.jl b/examples/treetensornetworks/solvers/01_tdvp.jl deleted file mode 100644 index af8943ca..00000000 --- a/examples/treetensornetworks/solvers/01_tdvp.jl +++ /dev/null @@ -1,42 +0,0 @@ -using ITensors -using ITensorNetworks - -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) -ψ = randomMPS(s, "↑"; linkdims=10) - -@show inner(ψ', H, ψ) / inner(ψ, ψ) - -ϕ = tdvp( - H, - -1.0, - ψ; - nsweeps=20, - reverse_step=false, - normalize=true, - maxdim=30, - cutoff=1e-10, - 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) - -ϕ3 = ITensorNetworks.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1) - -@show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3) diff --git a/examples/treetensornetworks/solvers/02_dmrg-x.jl b/examples/treetensornetworks/solvers/02_dmrg-x.jl deleted file mode 100644 index 0a17ad04..00000000 --- a/examples/treetensornetworks/solvers/02_dmrg-x.jl +++ /dev/null @@ -1,44 +0,0 @@ -using ITensors -using ITensorNetworks -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 - end - end - return os -end - -n = 10 -s = siteinds("S=1/2", n) - -using Random -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) - -dmrg_x_kwargs = ( - nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1 -) - -ϕ = 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 diff --git a/examples/treetensornetworks/solvers/03_models.jl b/examples/treetensornetworks/solvers/03_models.jl deleted file mode 100644 index 03a4bc97..00000000 --- a/examples/treetensornetworks/solvers/03_models.jl +++ /dev/null @@ -1,20 +0,0 @@ -using ITensors - -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/treetensornetworks/solvers/03_solvers.jl b/examples/treetensornetworks/solvers/03_solvers.jl deleted file mode 100644 index d8f00580..00000000 --- a/examples/treetensornetworks/solvers/03_solvers.jl +++ /dev/null @@ -1,37 +0,0 @@ -using DifferentialEquations -using ITensors -using ITensorNetworks -using KrylovKit: exponentiate - -function ode_solver( - H::TimeDependentSum, - time_step, - ψ₀; - current_time=0.0, - outputlevel=0, - solver_alg=Tsit5(), - kwargs..., -) - if outputlevel ≥ 3 - println(" In ODE solver, current_time = $current_time, time_step = $time_step") - end - - time_span = (current_time, current_time + time_step) - u₀, ITensor_from_vec = to_vec(ψ₀) - f(ψ::ITensor, p, t) = H(t)(ψ) - f(u::Vector, p, t) = to_vec(f(ITensor_from_vec(u), p, t))[1] - prob = ODEProblem(f, u₀, time_span) - sol = solve(prob, solver_alg; kwargs...) - uₜ = sol.u[end] - return ITensor_from_vec(uₜ), nothing -end - -function krylov_solver( - H::TimeDependentSum, time_step, ψ₀; current_time=0.0, outputlevel=0, kwargs... -) - if outputlevel ≥ 3 - println(" In Krylov solver, current_time = $current_time, time_step = $time_step") - end - ψₜ, info = exponentiate(H(current_time), time_step, ψ₀; kwargs...) - return ψₜ, info -end diff --git a/examples/treetensornetworks/solvers/03_tdvp_time_dependent.jl b/examples/treetensornetworks/solvers/03_tdvp_time_dependent.jl deleted file mode 100644 index 18b571a8..00000000 --- a/examples/treetensornetworks/solvers/03_tdvp_time_dependent.jl +++ /dev/null @@ -1,153 +0,0 @@ -using DifferentialEquations -using ITensors -using ITensorNetworks -using KrylovKit -using LinearAlgebra -using Random - -Random.seed!(1234) - -# Define the time-independent model -include("03_models.jl") - -# Define the solvers needed for TDVP -include("03_solvers.jl") - -# 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 -# solver. -outputlevel = 1 - -# 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 - -# 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 ω⃗] - -# 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 ℋ⃗₀] - -# Initial state, ψ₀ = ψ(0) -# Initialize as complex since that is what DifferentialEquations.jl -# expects. -ψ₀ = complex.(randomMPS(s, j -> isodd(j) ? "↑" : "↓"; linkdims=start_linkdim)) - -@show norm(ψ₀) - -println() -println("#"^100) -println("Running TDVP with ODE solver") -println("#"^100) -println() - -function ode_solver(H⃗₀, time_step, ψ₀; kwargs...) - return ode_solver( - -im * TimeDependentSum(f⃗, H⃗₀), - time_step, - ψ₀; - solver_alg=ode_alg, - ode_kwargs..., - kwargs..., - ) -end - -ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel) - -println() -println("Finished running TDVP with ODE solver") -println() - -println() -println("#"^100) -println("Running TDVP with Krylov solver") -println("#"^100) -println() - -function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...) - return krylov_solver( - -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs... - ) -end - -ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel) - -println() -println("Finished running TDVP with Krylov solver") -println() - -println() -println("#"^100) -println("Running full state evolution with ODE solver") -println("#"^100) -println() - -@disable_warn_order begin - ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel) -end - -println() -println("Finished full state evolution with ODE solver") -println() - -@show norm(ψₜ_ode) -@show norm(ψₜ_krylov) -@show norm(ψₜ_full) - -@show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full)) -@show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full)) diff --git a/examples/treetensornetworks/solvers/04_tdvp_observers.jl b/examples/treetensornetworks/solvers/04_tdvp_observers.jl deleted file mode 100644 index 9cb02c8e..00000000 --- a/examples/treetensornetworks/solvers/04_tdvp_observers.jl +++ /dev/null @@ -1,82 +0,0 @@ -using ITensors -using ITensorNetworks -using Observers - -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 -cutoff = 1e-12 -tau = 0.1 -ttotal = 1.0 - -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 - end - return nothing -end - -function current_time(; current_time, bond, half_sweep) - if bond == 1 && half_sweep == 2 - return current_time - end - return nothing -end - -function measure_sz(; psi, bond, half_sweep) - if bond == 1 && half_sweep == 2 - return expect(psi, "Sz"; vertices=[N ÷ 2]) - end - return nothing -end - -function return_state(; psi, bond, half_sweep) - if bond == 1 && half_sweep == 2 - return psi - end - return nothing -end - -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, -) - -res = results(obs) -steps = res["steps"] -times = res["times"] -psis = res["psis"] -Sz = res["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() -end diff --git a/examples/treetensornetworks/solvers/05_tdvp_nonuniform_timesteps.jl b/examples/treetensornetworks/solvers/05_tdvp_nonuniform_timesteps.jl deleted file mode 100644 index eeae88c0..00000000 --- a/examples/treetensornetworks/solvers/05_tdvp_nonuniform_timesteps.jl +++ /dev/null @@ -1,47 +0,0 @@ -using ITensors -using ITensorNetworks - -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 - 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 -) - -res = results(obs) -times = res["times"] -psis = res["psis"] - -println("\nResults") -println("=======") -print("step = ", 0) -print(", time = ", zero(ComplexF64)) -print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; vertices=[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"; vertices=[N ÷ 2]); digits=3)) - println() -end diff --git a/examples/treetensornetworks/solvers/05_utils.jl b/examples/treetensornetworks/solvers/05_utils.jl deleted file mode 100644 index a8c5feeb..00000000 --- a/examples/treetensornetworks/solvers/05_utils.jl +++ /dev/null @@ -1,60 +0,0 @@ -using ITensors -using ITensorNetworks -using Observers -using Printf - -using ITensorNetworks: tdvp_solver, tdvp_step, process_sweeps, TDVPOrder - -function tdvp_nonuniform_timesteps( - solver, - PH, - psi::MPS; - time_steps, - reverse_step=true, - time_start=0.0, - order=2, - (step_observer!)=Observer(), - kwargs..., -) - nsweeps = length(time_steps) - maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, kwargs...) - tdvp_order = TDVPOrder(order, Base.Forward) - current_time = time_start - for sw in 1:nsweeps - sw_time = @elapsed begin - psi, PH, info = tdvp_step( - tdvp_order, - solver, - PH, - time_steps[sw], - psi; - kwargs..., - current_time, - reverse_step, - sweep=sw, - maxdim=maxdim[sw], - mindim=mindim[sw], - cutoff=cutoff[sw], - noise=noise[sw], - ) - end - current_time += time_steps[sw] - - update!(step_observer!; psi, sweep=sw, outputlevel, current_time) - - if outputlevel ≥ 1 - print("After sweep ", sw, ":") - print(" maxlinkdim=", maxlinkdim(psi)) - @printf(" maxerr=%.2E", info.maxtruncerr) - print(" current_time=", round(current_time; digits=3)) - print(" time=", round(sw_time; digits=3)) - println() - flush(stdout) - end - end - return psi -end - -function tdvp_nonuniform_timesteps(H, psi::MPS; kwargs...) - return tdvp_nonuniform_timesteps(tdvp_solver(; kwargs...), H, psi; kwargs...) -end diff --git a/test/test_treetensornetworks/test_solvers/Project.toml b/test/test_treetensornetworks/test_solvers/Project.toml new file mode 100644 index 00000000..1f448702 --- /dev/null +++ b/test/test_treetensornetworks/test_solvers/Project.toml @@ -0,0 +1,9 @@ +[deps] +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" +Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 9943caa2..1bc39a19 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -1,12 +1,24 @@ -using ITensors -using ITensorNetworks -using ITensorNetworks: exponentiate_updater -using KrylovKit: exponentiate -using Observers -using Random -using Test - -#ToDo: Add tests for different signatures and functionality of extending the params +@eval module $(gensym()) +using Graphs: dst, edges, src +using ITensors: ITensor, contract, dag, inner, noprime, normalize, prime, scalar +using ITensorNetworks: + ITensorNetworks, + OpSum, + TTN, + apply, + expect, + mpo, + mps, + op, + random_mps, + random_ttn, + siteinds, + tdvp +using LinearAlgebra: norm +using NamedGraphs: named_binary_tree, named_comb_tree +using Observers: observer +using Test: @testset, @test + @testset "MPS TDVP" begin @testset "Basic TDVP" begin N = 10 @@ -193,7 +205,7 @@ using Test cutoff, normalize=false, updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), - updater=exponentiate_updater, + updater=ITensorNetworks.exponentiate_updater, ) # TODO: What should `expect` output? Right now # it outputs a dictionary. @@ -253,7 +265,6 @@ using Test for step in 1:Nsteps state = apply(gates, state; cutoff) - #normalize!(state) nsites = (step <= 3 ? 2 : 1) phi = tdvp( @@ -279,7 +290,7 @@ using Test phi = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) - obs = Observer( + obs = observer( "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c], "En" => (; state) -> real(inner(state', H, state)), ) @@ -361,12 +372,12 @@ using Test measure_sz(; state) = expect("Sz", state; vertices=[c])[c] measure_en(; state) = real(inner(state', H, state)) - sweep_obs = Observer("Sz" => measure_sz, "En" => measure_en) + sweep_obs = observer("Sz" => measure_sz, "En" => measure_en) get_info(; info) = info step_measure_sz(; state) = expect("Sz", state; vertices=[c])[c] step_measure_en(; state) = real(inner(state', H, state)) - region_obs = Observer( + region_obs = observer( "Sz" => step_measure_sz, "En" => step_measure_en, "info" => get_info ) @@ -411,7 +422,7 @@ end H = TTN(os, s) - ψ0 = normalize!(random_ttn(s)) + ψ0 = normalize(random_ttn(s)) # Time evolve forward: ψ1 = tdvp(H, -0.1im, ψ0; root_vertex, nsweeps=1, cutoff, nsites=2) @@ -453,7 +464,7 @@ end H2 = TTN(os2, s) Hs = [H1, H2] - ψ0 = normalize!(random_ttn(s; link_space=10)) + ψ0 = normalize(random_ttn(s; link_space=10)) ψ1 = tdvp(Hs, -0.1im, ψ0; nsweeps=1, cutoff, nsites=1) @@ -564,7 +575,6 @@ end for step in 1:Nsteps state = apply(gates, state; cutoff, maxdim) - #normalize!(state) nsites = (step <= 3 ? 2 : 1) phi = tdvp( @@ -589,7 +599,7 @@ end # phi = TTN(s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") - obs = Observer( + obs = observer( "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c], "En" => (; state) -> real(inner(state', H, state)), ) @@ -622,7 +632,7 @@ end os = ITensorNetworks.heisenberg(c) H = TTN(os, s) - state = normalize!(random_ttn(s; link_space=2)) + state = normalize(random_ttn(s; link_space=2)) trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) @@ -641,100 +651,5 @@ end @test inner(state', H, state) < -2.47 end - - # TODO: verify quantum number suport in ITensorNetworks - - # @testset "Observers" begin - # cutoff = 1e-12 - # tau = 0.1 - # ttotal = 1.0 - - # tooth_lengths = fill(2, 3) - # c = named_comb_tree(tooth_lengths) - # s = siteinds("S=1/2", c; conserve_qns=true) - - # os = ITensorNetworks.heisenberg(c) - # H = TTN(os, s) - - # c = (2, 2) - - # # - # # Using the ITensors observer system - # # - # struct TDVPObserver <: AbstractObserver end - - # Nsteps = convert(Int, ceil(abs(ttotal / tau))) - # Sz1 = zeros(Nsteps) - # En1 = zeros(Nsteps) - # function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...) - # if bond == 1 && half_sweep == 2 - # Sz1[sweep] = expect("Sz", psi; vertices=[c])[c] - # En1[sweep] = real(inner(psi', H, psi)) - # end - # end - - # psi1 = productMPS(s, n -> isodd(n) ? "Up" : "Dn") - # tdvp( - # H, - # -im * ttotal, - # psi1; - # time_step=-im * tau, - # cutoff, - # normalize=false, - # (observer!)=TDVPObserver(), - # root_vertex=N, - # ) - - # # - # # Using Observers.jl - # # - - # function measure_sz(; psi, bond, half_sweep) - # if bond == 1 && half_sweep == 2 - # return expect("Sz", psi; vertices=[c])[c] - # end - # return nothing - # end - - # function measure_en(; psi, bond, half_sweep) - # if bond == 1 && half_sweep == 2 - # return real(inner(psi', H, psi)) - # end - # return nothing - # end - - # obs = Observer("Sz" => measure_sz, "En" => measure_en) - - # step_measure_sz(; psi) = expect("Sz", psi; vertices=[c])[c] - - # step_measure_en(; psi) = real(inner(psi', H, psi)) - - # step_obs = Observer("Sz" => step_measure_sz, "En" => step_measure_en) - - # psi2 = MPS(s, n -> isodd(n) ? "Up" : "Dn") - # tdvp( - # H, - # -im * ttotal, - # psi2; - # time_step=-im * tau, - # cutoff, - # normalize=false, - # (observer!)=obs, - # (step_observer!)=step_obs, - # root_vertex=N, - # ) - - # Sz2 = results(obs)["Sz"] - # En2 = results(obs)["En"] - - # Sz2_step = results(step_obs)["Sz"] - # En2_step = results(step_obs)["En"] - - # @test Sz1 ≈ Sz2 - # @test En1 ≈ En2 - # @test Sz1 ≈ Sz2_step - # @test En1 ≈ En2_step - # end end - -nothing +end