windowhamiltonian poc

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")

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)

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)

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}}

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.:+(

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)

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
 
 #===========================================================================================

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
+

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)