diff --git a/Project.toml b/Project.toml index a44bfa7..f16f2e6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorTDVP" uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342" authors = ["Matthew Fishman and contributors"] -version = "0.2.6" +version = "0.3.0" [deps] ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" diff --git a/README.md b/README.md index 2f23235..098670b 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,21 @@ pkg> add ITensorTDVP ## News +### ITensorTDVP.jl v0.3 Release Notes + +#### Breaking changes + +- `ITensorTDVP.dmrg` and `ITensorTDVP.dmrg_x` now output a tuple containing the eigenvalue and eigenvector, while before they just output the eigenvector. You should update code like this: +```julia +psi = dmrg_x(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6) +psi = ITensorTDVP.dmrg(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6) +``` +to: +```julia +energy, psi = dmrg_x(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6) +energy, psi = ITensorTDVP.dmrg(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6) +``` + ### ITensorTDVP.jl v0.2 Release Notes #### Breaking changes diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl index c386c0a..26e6dc2 100644 --- a/examples/01_tdvp.jl +++ b/examples/01_tdvp.jl @@ -36,11 +36,11 @@ function main() e2, ϕ2 = ITensors.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), e2 - ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1) + e3, ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1) - @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3) + @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3), e3 return nothing end diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl index 6d9e670..77fc311 100644 --- a/examples/02_dmrg-x.jl +++ b/examples/02_dmrg-x.jl @@ -36,11 +36,11 @@ function main() nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1 ) - ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...) + e, ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...) @show inner(ψ', H, ψ) / inner(ψ, ψ) @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2 - @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ) + @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ), e @show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2 return nothing end diff --git a/examples/Project.toml b/examples/Project.toml index d262170..2e1c761 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,5 +1,6 @@ [deps] ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/src/dmrg.jl b/src/dmrg.jl index f736f88..ef6f412 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -21,8 +21,7 @@ function dmrg_solver( maxiter=solver_maxiter, verbosity=solver_verbosity, ) - psi = vecs[1] - return psi, info + return vecs[1], (; info, eigval=vals[1]) end return solver end @@ -36,9 +35,13 @@ function dmrg( solver_krylovdim=default_solver_krylovdim(eigsolve), solver_maxiter=default_solver_maxiter(eigsolve), solver_verbosity=default_solver_verbosity(), + (observer!)=default_observer!(), kwargs..., ) reverse_step = false + info_ref! = Ref{Any}() + info_observer! = values_observer(; info=info_ref!) + observer! = compose_observers(observer!, info_observer!) psi = alternating_update( dmrg_solver( eigsolve; @@ -52,7 +55,8 @@ function dmrg( H, psi0; reverse_step, + observer!, kwargs..., ) - return psi + return info_ref![].eigval, psi end diff --git a/src/dmrg_x.jl b/src/dmrg_x.jl index 7a82deb..d640237 100644 --- a/src/dmrg_x.jl +++ b/src/dmrg_x.jl @@ -5,12 +5,19 @@ function dmrg_x_solver(PH, t, psi0; current_time, outputlevel) H = contract(PH, ITensor(true)) D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) + u′ = uniqueind(D, U) max_overlap, max_ind = findmax(abs, array(psi0 * dag(U))) U_max = U * dag(onehot(eltype(U), u => max_ind)) - return U_max, nothing + D_max = D[u′ => max_ind, u => max_ind] + return U_max, (; eigval=D_max) end -function dmrg_x(PH, psi0::MPS; reverse_step=false, kwargs...) - psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, kwargs...) - return psi +function dmrg_x( + PH, psi0::MPS; reverse_step=false, (observer!)=default_observer!(), kwargs... +) + info_ref = Ref{Any}() + info_observer! = values_observer(; info=info_ref) + observer! = compose_observers(observer!, info_observer!) + psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, observer!, kwargs...) + return info_ref[].eigval, psi end diff --git a/src/update_observer.jl b/src/update_observer.jl index d4bb6b4..beed55d 100644 --- a/src/update_observer.jl +++ b/src/update_observer.jl @@ -6,3 +6,25 @@ using ITensors.ITensorMPS: ITensorMPS function update_observer!(observer::ITensorMPS.AbstractObserver; kwargs...) return ITensorMPS.measure!(observer; kwargs...) end + +struct ValuesObserver{Values<:NamedTuple} + values::Values +end +function update_observer!(observer::ValuesObserver; kwargs...) + for key in keys(observer.values) + observer.values[key][] = kwargs[key] + end + return observer +end +values_observer(; kwargs...) = ValuesObserver(NamedTuple(kwargs)) + +struct ComposedObservers{Observers<:Tuple} + observers::Observers +end +compose_observers(observers...) = ComposedObservers(observers) +function update_observer!(observer::ComposedObservers; kwargs...) + for observerᵢ in observer.observers + update_observer!(observerᵢ; kwargs...) + end + return observer +end diff --git a/test/test_dmrg.jl b/test/test_dmrg.jl index 159b73b..8a60dac 100644 --- a/test/test_dmrg.jl +++ b/test/test_dmrg.jl @@ -21,12 +21,14 @@ using Test: @test, @testset psi = randomMPS(elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=20) nsweeps = 10 maxdim = [10, 20, 40, 100] - psi = ITensorTDVP.dmrg( + e, psi = ITensorTDVP.dmrg( H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1 ) + @test inner(psi', H, psi) ≈ e e2, psi2 = ITensors.dmrg(H, psi; nsweeps, maxdim, cutoff, outputlevel=0) @test ITensors.scalartype(psi2) == elt @test e2 isa real(elt) + @test e ≈ e2 rtol = √(eps(real(elt))) * 10 @test inner(psi', H, psi) ≈ inner(psi2', H, psi2) rtol = √(eps(real(elt))) * 10 end end diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl index 8e6d7f4..34ce26b 100644 --- a/test/test_dmrg_x.jl +++ b/test/test_dmrg_x.jl @@ -34,13 +34,15 @@ using Test: @test, @testset dmrg_x_kwargs = ( nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 ) - ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...) + e, ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...) @test ITensors.scalartype(ϕ) == elt + @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e @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...) + ẽ, ϕ̃ = dmrg_x(ProjMPO(H), ϕ; nsite=1, dmrg_x_kwargs...) @test ITensors.scalartype(ϕ̃) == elt + @test inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) ≈ ẽ @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1 scale(::Type{Float32}) = 10^2 scale(::Type{Float64}) = 10^5 diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index 42d0902..086f9f5 100644 --- a/test/test_linsolve.jl +++ b/test/test_linsolve.jl @@ -22,7 +22,7 @@ using Random: Random state = [isodd(n) ? "Up" : "Dn" for n in 1:N] Random.seed!(1234) x_c = randomMPS(elt, s, state; linkdims=2) - x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0) + e, 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)