diff --git a/src/MPSKit.jl b/src/MPSKit.jl index c1b47dc40..0137080a4 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -19,7 +19,7 @@ export QP, LeftGaugedQP, RightGaugedQP export AbstractMPO export MPO, FiniteMPO, InfiniteMPO export JordanMPOTensor, JordanMPOTensorMap -export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian +export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian, WindowMPOHamiltonian export MultilineMPO export UntimedOperator, TimedOperator, MultipliedOperator, LazySum @@ -122,6 +122,7 @@ include("operators/abstractmpo.jl") include("operators/mpo.jl") include("operators/jordanmpotensor.jl") include("operators/mpohamiltonian.jl") # the mpohamiltonian objects +include("operators/windowhamiltonian.jl") include("operators/ortho.jl") include("operators/multilinempo.jl") include("operators/projection.jl") diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index b5db28f12..aff3e375f 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -128,6 +128,13 @@ function expectation_value( return dot(ψ, H, ψ, envs) / dot(ψ, ψ) end +function expectation_value( + ψ::WindowMPS, H::WindowMPOHamiltonian, + envs::AbstractMPSEnvironments = environments(ψ, H) + ) + return dot(ψ, H, ψ, envs) / dot(ψ, ψ) +end + function expectation_value( ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, envs::AbstractMPSEnvironments = environments(ψ, H) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index f4936c002..7f8831726 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -35,7 +35,7 @@ end """ propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init=copy(ψ₀)) -Calculate the propagator ``\\frac{1}{E₀ + z - H}|ψ₀⟩`` using the dynamical DMRG +Calculate the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG algorithm. """ function propagator end @@ -47,19 +47,17 @@ An alternative approach to the dynamical DMRG algorithm, without quadratic terms less controlled approximation. This algorithm minimizes the following cost function ```math -⟨ψ|(H - E)|ψ⟩ - ⟨ψ|ψ₀⟩ - ⟨ψ₀|ψ⟩ -``` -which is equivalent to the original approach if -```math -|ψ₀⟩ = (H - E)|ψ⟩ +⟨ψ|(z - H)|ψ⟩ - ⟨ψ|ψ₀⟩ - ⟨ψ₀|ψ⟩ ``` +Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩. + See also [`Jeckelmann`](@ref) for the original approach. """ struct NaiveInvert <: DDMRG_Flavour end function propagator( - A::AbstractFiniteMPS, z::Number, H::FiniteMPOHamiltonian, + A::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG{NaiveInvert}; init = copy(A) ) h_envs = environments(init, H) # environments for h @@ -105,11 +103,20 @@ end """ $(TYPEDEF) -The original flavour of dynamical DMRG, which minimizes the following (quadratic) cost function: +The original flavour of dynamical DMRG, which minimizes functional (14) from Jeckelmann2002. +```math +|| <ψ| (Re(z) - H)^2 + Im(z)^2 |ψ⟩ +Im(z) (<ψ₀|ψ⟩+<ψ|ψ₀⟩ || +``` + +This would achieve a minima at ```math -|| (H - E) |ψ₀⟩ - |ψ⟩ || +-Im(z) |ψ₀⟩ = ((Re(z) - H)^2 + Im(z)^2)|ψ⟩ ``` +Together with equation (11) from that same paper we can determine the full propagator (z-H)^-1 |ψ₀>. + +Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩. + See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative. ## References @@ -119,7 +126,7 @@ See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative. struct Jeckelmann <: DDMRG_Flavour end function propagator( - A::AbstractFiniteMPS, z, H::FiniteMPOHamiltonian, + A::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG{Jeckelmann}; init = copy(A) ) ω = real(z) @@ -176,7 +183,7 @@ function propagator( end function squaredenvs( - state::AbstractFiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H) + state::AbstractFiniteMPS, H, envs = environments(state, H) ) H² = conj(H) * H L = length(state) diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index 6372f65ab..9194e5cb5 100644 --- a/src/environments/finite_envs.jl +++ b/src/environments/finite_envs.jl @@ -53,14 +53,14 @@ function environments( return environments(below, O, above, leftstart, rightstart) end function environments( - below::WindowMPS, O::Union{InfiniteMPOHamiltonian, InfiniteMPO}, above = nothing; - lenvs = environments(below.left_gs, O), - renvs = environments(below.right_gs, O) + below::WindowMPS, O::WindowMPOHamiltonian, above = nothing; + lenvs = environments(below.left_gs, O.left_ham), + renvs = environments(below.right_gs, O.right_ham) ) leftstart = copy(lenvs.GLs[1]) rightstart = copy(renvs.GRs[end]) - return environments(below, O, above, leftstart, rightstart) + return environments(below, O.finite_ham, above, leftstart, rightstart) end function environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index daf7eb6a5..f2008f755 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -706,7 +706,7 @@ end function Base.similar(H::MPOHamiltonian, ::Type{T}) where {T <: Number} return MPOHamiltonian(similar.(parent(H), T)) end - +Base.circshift(H::InfiniteMPOHamiltonian, shift::Integer) = InfiniteMPOHamiltonian(circshift(parent(copy(H)), shift)) # Linear Algebra # -------------- function Base.:+( diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl new file mode 100644 index 000000000..071e0eca3 --- /dev/null +++ b/src/operators/windowhamiltonian.jl @@ -0,0 +1,40 @@ +""" +A WindowMPS is a finite MPS embedded between an infinite mps to the left, and an inifite mps to the right. +The associated hamiltonian has also an infinite part to the left, a finite hamiltonian in the middle, and an infinite part to the right. + +Acts simalar as just a finite hamiltonian, but we 'remember' the boundary hamiltonians. +""" + +# todo - what is the required interface for abstractmpo? +# support densempo windows? +struct WindowMPOHamiltonian{O} <: AbstractMPO{O} + left_ham::InfiniteMPOHamiltonian{O} + finite_ham::FiniteMPOHamiltonian{O} + right_ham::InfiniteMPOHamiltonian{O} +end + +#utility constructor +function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange) + + # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly + left_edge = (interval.start - 1) % length(ham) + left_ham = circshift(ham, -left_edge) + right_edge = (interval.stop + 1) % length(ham) + right_ham = circshift(ham, -right_edge + 1) + + finite_ham = FiniteMPOHamiltonian([ham[i] for i in interval]) + return WindowMPOHamiltonian(left_ham, finite_ham, right_ham) +end + +Base.parent(h::WindowMPOHamiltonian) = h.finite_ham +Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham)) + +# some basic linalg +for fun in (:(Base.:+), :(Base.:-), :(Base.:*)) + @eval $fun(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham, b.left_ham), $fun(a.finite_ham, b.finite_ham), $fun(a.right_ham, b.right_ham)) +end + +TensorKit.dot( + bra::WindowMPS, H::WindowMPOHamiltonian, ket::WindowMPS = bra, + envs = environments(bra, H, ket) +) = dot(bra.window, H.finite_ham, ket.window, envs) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 5c0bab9d0..7c5d48327 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -106,18 +106,27 @@ function WindowMPS(N::Int, V::VectorSpace, args...; kwargs...) return WindowMPS(fill(V, N), args...; kwargs...) end -function WindowMPS(ψ::InfiniteMPS{A, B}, L::Int) where {A, B} + +WindowMPS(ψ::InfiniteMPS, L::Int) = WindowMPS(ψ, 1:L) + +function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A, B} + + # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly + left_edge = (interval.start - 1) % length(ψ) + right_edge = (interval.stop + 1) % length(ψ) + + L = length(interval) CLs = Vector{Union{Missing, B}}(missing, L + 1) ALs = Vector{Union{Missing, A}}(missing, L) ARs = Vector{Union{Missing, A}}(missing, L) ACs = Vector{Union{Missing, A}}(missing, L) - ALs .= ψ.AL[1:L] - ARs .= ψ.AR[1:L] - ACs .= ψ.AC[1:L] - CLs .= ψ.C[0:L] + ALs .= ψ.AL[interval] + ARs .= ψ.AR[interval] + ACs .= ψ.AC[interval] + CLs .= ψ.C[(interval.start - 1):interval.stop] - return WindowMPS(ψ, FiniteMPS(ALs, ARs, ACs, CLs), ψ) + return WindowMPS(circshift(ψ, -left_edge), FiniteMPS(ALs, ARs, ACs, CLs), circshift(ψ, -right_edge + 1)) end #=========================================================================================== diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl index 2b7a821b7..146c43af3 100644 --- a/test/algorithms/dynamical_dmrg.jl +++ b/test/algorithms/dynamical_dmrg.jl @@ -12,7 +12,7 @@ using TensorKit: ℙ verbosity_conv = 1 -@testset "Dynamical DMRG" verbose = true begin +@testset "Dynamical DMRG FiniteMPS" verbose = true begin L = 10 H = force_planar(-transverse_field_ising(; L, g = -4)) gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^10), H; verbosity = verbosity_conv) @@ -32,3 +32,32 @@ verbosity_conv = 1 @test data ≈ predicted atol = 1.0e-8 end end + + +@testset "Dynamical DMRG WindowMPS" verbose = true begin + N = 20 + + H = transverse_field_ising(g = -4) + Ω = InfiniteMPS(ComplexSpace(2), ComplexSpace(20)) + + (Ω, _) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv)) + XΩ = WindowMPS(Ω, N) + H_w = WindowMPOHamiltonian(H, 1:N) + + + gs_en = expectation_value(XΩ, H_w) + vals = range(gs_en - 1.0, gs_en + 1.0, length = 5) + eta = 0.3im + predicted = [1 / (v + eta - gs_en) for v in vals] + + + @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert()) + alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8) + data = map(vals) do v + result, = propagator(XΩ, v + eta, H_w, alg) + return result + end + @test data ≈ predicted atol = 1.0e-8 + end +end + diff --git a/test/states/windowmps.jl b/test/states/windowmps.jl index 3cbddf09d..ab37bbab8 100644 --- a/test/states/windowmps.jl +++ b/test/states/windowmps.jl @@ -63,14 +63,15 @@ using TensorKit: ℙ e1 = expectation_value(window, (2, 3) => O) - window, envs, _ = find_groundstate(window, ham, DMRG(; verbosity = 0)) + w_ham = WindowMPOHamiltonian(ham, 1:length(window)) + window, envs, _ = find_groundstate(window, w_ham, DMRG(; verbosity = 0)) e2 = expectation_value(window, (2, 3) => O) @test real(e2) ≤ real(e1) - window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme = truncrank(20)), envs) - window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs) + window, envs = timestep(window, w_ham, 0.1, 0.0, TDVP2(; trscheme = truncrank(20)), envs) + window, envs = timestep(window, w_ham, 0.1, 0.0, TDVP(), envs) e3 = expectation_value(window, (2, 3) => O)