diff --git a/test/test_contract_mpo.jl b/test/test_contract_mpo.jl index ab33cc6..549b3e6 100644 --- a/test/test_contract_mpo.jl +++ b/test/test_contract_mpo.jl @@ -13,7 +13,7 @@ using ITensors: truncate! using ITensorTDVP: ITensorTDVP using Test: @test, @testset -@testset "Contract MPO (eltype=$elt)" for elt in ( +@testset "Contract MPO (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( Float32, Float64, Complex{Float32}, Complex{Float64} ), conserve_qns in [false, true] @@ -52,11 +52,13 @@ using Test: @test, @testset truncate!(psi_guess; maxdim=2) Hpsi = apply(H, psi; alg="fit", nsweeps=4, init_mps=psi_guess) @test ITensors.scalartype(Hpsi) == elt - @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1e-5 + @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = √eps(real(elt)) # Test with nsite=1 - Hpsi_guess = apply(H, psi; alg="naive", cutoff=1E-4) + Hpsi_guess = apply(H, psi; alg="naive", cutoff=1e-4) Hpsi = apply(H, psi; alg="fit", init_mps=Hpsi_guess, nsite=1, nsweeps=2) @test ITensors.scalartype(Hpsi) == elt - @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1e-4 + scale(::Type{Float32}) = 10^2 + scale(::Type{Float64}) = 10^6 + @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = √eps(real(elt)) * scale(real(elt)) end end diff --git a/test/test_dmrg.jl b/test/test_dmrg.jl index c31e175..159b73b 100644 --- a/test/test_dmrg.jl +++ b/test/test_dmrg.jl @@ -2,7 +2,7 @@ using ITensors: ITensors, MPO, OpSum, inner, randomMPS, siteinds using ITensorTDVP: ITensorTDVP using Test: @test, @testset -@testset "DMRG (eltype=$elt, nsite=$nsite)" for elt in ( +@testset "DMRG (eltype=$elt, nsite=$nsite, conserve_qns=$conserve_qns)" for elt in ( Float32, Float64, Complex{Float32}, Complex{Float64} ), nsite in [1, 2], diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl index 8bd2cc2..8e6d7f4 100644 --- a/test/test_dmrg_x.jl +++ b/test/test_dmrg_x.jl @@ -1,8 +1,8 @@ -using ITensors -using ITensorTDVP -using Random -using Test - +@eval module $(gensym()) +using ITensors: ITensors, MPO, MPS, OpSum, ProjMPO, inner, siteinds +using ITensorTDVP: dmrg_x +using Random: Random +using Test: @test, @testset @testset "DMRG-X (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( Float32, Float64, Complex{Float32}, Complex{Float64} ), @@ -35,24 +35,18 @@ using Test nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 ) ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...) - @test ITensors.scalartype(ϕ) == elt @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = √eps(real(elt)) @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = √eps(real(elt)) - ϕ̃ = dmrg_x(ProjMPO(H), ϕ; nsite=1, dmrg_x_kwargs...) - @test ITensors.scalartype(ϕ̃) == elt @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1 - scale(::Type{Float32}) = 10^2 scale(::Type{Float64}) = 10^5 - @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = √(eps(real(elt))) * scale(real(elt)) # Sometimes broken, sometimes not # @test abs(loginner(ϕ̃, ϕ) / n) ≈ 0.0 atol = 1e-6 end - -nothing +end diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index 32f3efe..42d0902 100644 --- a/test/test_linsolve.jl +++ b/test/test_linsolve.jl @@ -1,50 +1,41 @@ -using ITensors -using ITensorTDVP -using Test -using Random - -@testset "linsolve with conserve_qns=$conserve_qns and eltype=$eltype" for conserve_qns in - (false, true), - eltype in (Float32, Float64, ComplexF32, ComplexF64) +@eval module $(gensym()) +using ITensors: ITensors, MPO, OpSum, apply, randomMPS, siteinds +using ITensorTDVP: ITensorTDVP, dmrg +using KrylovKit: linsolve +using LinearAlgebra: norm +using Test: @test, @testset +using Random: Random +@testset "linsolve (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] N = 6 s = siteinds("S=1/2", N; conserve_qns) - 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 - H = ITensors.convert_leaf_eltype(eltype, MPO(os, s)) - + H = ITensors.convert_leaf_eltype(elt, MPO(os, s)) state = [isodd(n) ? "Up" : "Dn" for n in 1:N] - Random.seed!(1234) - x_c = randomMPS(eltype, s, state; linkdims=2) - e, x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0) - - @test ITensors.scalartype(x_c) == eltype - + x_c = randomMPS(elt, s, state; linkdims=2) + x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0) + @test ITensors.scalartype(x_c) == elt # Compute `b = H * x_c` b = apply(H, x_c; cutoff=1e-8) - - @test ITensors.scalartype(b) == eltype - + @test ITensors.scalartype(b) == elt # Starting guess - x0 = x_c + eltype(0.05) * randomMPS(eltype, s, state; linkdims=2) - - @test ITensors.scalartype(x0) == eltype - + x0 = x_c + elt(0.05) * randomMPS(elt, s, state; linkdims=2) + @test ITensors.scalartype(x0) == elt nsweeps = 10 cutoff = 1e-5 maxdim = 20 solver_kwargs = (; tol=1e-4, maxiter=20, krylovdim=30, ishermitian=true) x = linsolve(H, b, x0; nsweeps, cutoff, maxdim, solver_kwargs) - - @test ITensors.scalartype(x) == eltype - + @test ITensors.scalartype(x) == elt @test norm(x - x_c) < 1e-2 end - -nothing +end diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index e5dd9ed..38bcb32 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -1,10 +1,11 @@ -using ITensors -using ITensorTDVP -using KrylovKit -using Observers -using Random -using Test - +@eval module $(gensym()) +using ITensors: + ITensors, AbstractObserver, ITensor, MPO, MPS, OpSum, randomMPS, inner, op, siteinds +using ITensorTDVP: tdvp +using KrylovKit: exponentiate +using LinearAlgebra: norm +using Observers: observer +using Test: @test, @testset @testset "Basic TDVP (eltype=$elt)" for elt in ( Float32, Float64, Complex{Float32}, Complex{Float64} ) @@ -20,31 +21,22 @@ using Test H = MPO(elt, os, s) ψ0 = randomMPS(elt, s; linkdims=10) time_step = elt(0.1) * im - # Time evolve forward: ψ1 = tdvp(H, -time_step, ψ0; nsweeps=1, cutoff, nsite=1) - @test ITensors.scalartype(ψ1) == complex(elt) @test ψ1 ≈ tdvp(-time_step, H, ψ0; nsweeps=1, cutoff, nsite=1) @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1) - #Different backend solvers, default solver_backend = "exponentiate" @test ψ1 ≈ tdvp(H, ψ0, -time_step; nsweeps=1, cutoff, nsite=1, solver_backend="applyexp") - @test norm(ψ1) ≈ 1 rtol = eps(real(elt)) * 10^3 - ## Should lose fidelity: #@test abs(inner(ψ0,ψ1)) < 0.9 - # Average energy should be conserved: @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) - # Time evolve backwards: ψ2 = tdvp(H, time_step, ψ1; nsweeps=1, cutoff) - @test ITensors.scalartype(ψ2) == complex(elt) @test norm(ψ2) ≈ 1 rtol = eps(real(elt)) * 10^4 - # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end @@ -64,52 +56,36 @@ end for j in 1:(N - 1) os2 += "Sz", j, "Sz", j + 1 end - H1 = MPO(os1, s) H2 = MPO(os2, s) Hs = [H1, H2] - ψ0 = randomMPS(s; linkdims=10) - ψ1 = tdvp(Hs, -0.1im, ψ0; nsweeps=1, cutoff, nsite=1) - @test ψ1 ≈ tdvp(-0.1im, Hs, ψ0; nsweeps=1, cutoff, nsite=1) @test ψ1 ≈ tdvp(Hs, ψ0, -0.1im; nsweeps=1, cutoff, nsite=1) - @test norm(ψ1) ≈ 1.0 - ## Should lose fidelity: #@test abs(inner(ψ0,ψ1)) < 0.9 - # Average energy should be conserved: @test real(sum(H -> inner(ψ1', H, ψ1), Hs)) ≈ sum(H -> inner(ψ0', H, ψ0), Hs) - # Time evolve backwards: ψ2 = tdvp(Hs, +0.1im, ψ1; nsweeps=1, cutoff) - @test norm(ψ2) ≈ 1.0 - # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - @testset "Custom solver in TDVP" begin N = 10 cutoff = 1e-12 - s = siteinds("S=1/2", 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 - H = MPO(os, s) - ψ0 = randomMPS(s; linkdims=10) - function solver(PH, t, psi0; kwargs...) solver_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true @@ -117,34 +93,24 @@ end psi, info = exponentiate(PH, t, psi0; solver_kwargs...) return psi, info end - ψ1 = tdvp(solver, H, -0.1im, ψ0; cutoff, nsite=1) - @test norm(ψ1) ≈ 1.0 - ## Should lose fidelity: #@test abs(inner(ψ0,ψ1)) < 0.9 - # Average energy should be conserved: @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) - # Time evolve backwards: ψ2 = tdvp(H, +0.1im, ψ1; cutoff) - @test norm(ψ2) ≈ 1.0 - # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - @testset "Accuracy Test" begin N = 4 tau = 0.1 ttotal = 1.0 cutoff = 1e-12 - s = siteinds("S=1/2", N; conserve_qns=false) - os = OpSum() for j in 1:(N - 1) os += 0.5, "S+", j, "S-", j + 1 @@ -153,20 +119,15 @@ end end H = MPO(os, s) HM = prod(H) - Ut = exp(-im * tau * HM) - psi = MPS(s, n -> isodd(n) ? "Up" : "Dn") psi2 = deepcopy(psi) psix = prod(psi) - Sz_tdvp = Float64[] Sz_tdvp2 = Float64[] Sz_exact = Float64[] - c = div(N, 2) Szc = op("Sz", s[c]) - Nsteps = Int(ttotal / tau) for step in 1:Nsteps psix = noprime(Ut * psix) @@ -183,7 +144,6 @@ end solver_krylovdim=25, ) push!(Sz_tdvp, real(expect(psi, "Sz"; sites=c:c)[1])) - psi2 = tdvp( H, -im * tau, @@ -195,15 +155,12 @@ end solver_krylovdim=25, ) push!(Sz_tdvp2, real(expect(psi2, "Sz"; sites=c:c)[1])) - push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix))) F = abs(scalar(dag(psix) * prod(psi))) end - @test norm(Sz_tdvp - Sz_exact) < 1e-5 @test norm(Sz_tdvp2 - Sz_exact) < 1e-5 end - @testset "TEBD Comparison" begin N = 10 cutoff = 1e-12 @@ -232,17 +189,12 @@ end psi = MPS(s, n -> isodd(n) ? "Up" : "Dn") phi = deepcopy(psi) c = div(N, 2) - - # # Evolve using TEBD - # - Nsteps = convert(Int, ceil(abs(ttotal / tau))) Sz1 = zeros(Nsteps) En1 = zeros(Nsteps) Sz2 = zeros(Nsteps) En2 = zeros(Nsteps) - for step in 1:Nsteps psi = apply(gates, psi; cutoff) nsite = (step <= 3 ? 2 : 1) @@ -254,10 +206,7 @@ end En1[step] = real(inner(psi', H, psi)) En2[step] = real(inner(phi', H, phi)) end - - # # Evolve using TDVP - # struct TDVPObserver <: AbstractObserver end Sz2 = zeros(Nsteps) En2 = zeros(Nsteps) @@ -280,24 +229,19 @@ end @test norm(Sz1 - Sz2) < 1e-3 @test norm(En1 - En2) < 1e-3 end - @testset "Imaginary Time Evolution" for reverse_step in [true, false] N = 10 cutoff = 1e-12 tau = 1.0 ttotal = 50.0 - s = siteinds("S=1/2", 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 - H = MPO(os, s) - psi = randomMPS(s; linkdims=2) psi2 = deepcopy(psi) trange = 0.0:tau:ttotal @@ -310,23 +254,18 @@ end H, -tau, psi2; cutoff, nsite, reverse_step, normalize=true, solver_krylovdim=15 ) end - @test psi ≈ psi2 rtol = 1e-6 - en1 = inner(psi', H, psi) en2 = inner(psi2', H, psi2) @test en1 < -4.25 @test en1 ≈ en2 end - @testset "Observers" begin N = 10 cutoff = 1e-12 tau = 0.1 ttotal = 1.0 - s = siteinds("S=1/2", N; conserve_qns=true) - os = OpSum() for j in 1:(N - 1) os += 0.5, "S+", j, "S-", j + 1 @@ -334,14 +273,9 @@ end os += "Sz", j, "Sz", j + 1 end H = MPO(os, s) - c = div(N, 2) - - # # Using the ITensors observer system - # struct TDVPObserver <: AbstractObserver end - Nsteps = convert(Int, ceil(abs(ttotal / tau))) Sz1 = zeros(Nsteps) En1 = zeros(Nsteps) @@ -351,7 +285,6 @@ end En1[sweep] = real(inner(psi', H, psi)) end end - psi1 = MPS(s, n -> isodd(n) ? "Up" : "Dn") tdvp( H, @@ -362,37 +295,26 @@ end normalize=false, (observer!)=TDVPObserver(), ) - - # # Using Observers.jl - # - function measure_sz(; psi, bond, half_sweep) if bond == 1 && half_sweep == 2 return expect(psi, "Sz"; sites=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 - function identity_info(; info) return info end - obs = observer("Sz" => measure_sz, "En" => measure_en, "info" => identity_info) - step_measure_sz(; psi) = expect(psi, "Sz"; sites=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, @@ -404,14 +326,11 @@ end (observer!)=obs, (step_observer!)=step_obs, ) - Sz2 = filter(!isnothing, obs.Sz) En2 = filter(!isnothing, obs.En) infos = obs.info - Sz2_step = step_obs.Sz En2_step = step_obs.En - @test Sz1 ≈ Sz2 @test En1 ≈ En2 @test Sz1 ≈ Sz2_step @@ -419,5 +338,4 @@ end @test all(x -> x.converged == 1, infos) @test length(values(infos)) == 180 end - -nothing +end