Add WindowMPOHamiltonian enabling propagators and dynamical DMRG on WindowMPS

src/MPSKit.jl

@@ -19,7 +19,7 @@ export QP, LeftGaugedQP, RightGaugedQP
 export AbstractMPO
 export MPO, FiniteMPO, InfiniteMPO
 export JordanMPOTensor
-export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian
+export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian, WindowMPOHamiltonian
 export MultilineMPO
 export UntimedOperator, TimedOperator, MultipliedOperator, LazySum
 
@@ -123,6 +123,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,8 +35,11 @@ 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 action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
 algorithm.
+
+Returns a tuple `(g, ψ)` where `g` is the approximation of the propagator matrix element
+``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and `ψ` is the MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``.
 """
 function propagator end
 
@@ -47,19 +50,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 ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
+
 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, init) # environments for h
@@ -105,10 +106,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.
+Writing ``ω = \\mathrm{Re}(z)`` and ``η = \\mathrm{Im}(z)``, this is
 ```math
-|| (H - E) |ψ₀⟩ - |ψ⟩ ||
+W(ψ) = ⟨ψ|(ω - H)^2 + η^2|ψ⟩ + η(⟨ψ₀|ψ⟩ + ⟨ψ|ψ₀⟩)
 ```
+which attains its minimum at
+```math
+((ω - H)^2 + η^2)|ψ⟩ = -η|ψ₀⟩
+```
+
+Together with equation (11) from that same paper we can determine the full propagator
+``\\frac{1}{z - H}|ψ₀⟩``.
+
+Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
 
 See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 
@@ -119,7 +130,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 +187,7 @@ function propagator(
 end
 
 function squaredenvs(
-        state::AbstractFiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H, state)
+        state::AbstractFiniteMPS, H, envs = environments(state, H, state)
     )
     H² = conj(H) * H
     L = length(state)

src/environments/finite_envs.jl

@@ -56,6 +56,17 @@ function environments(
     )
     return initialize_environments(below, operator, above, leftstart, rightstart)
 end
+function environments(
+        below::WindowMPS, O::WindowMPOHamiltonian, above = nothing;
+        # NOTE: the defaults are deliberately inlined rather than shared via `lenvs`/`renvs`
+        # kwargs: when explicit boundaries are passed (e.g. the squared environments in
+        # `squaredenvs`), we must not trigger the infinite fixed-point solve, which need not
+        # converge for e.g. `conj(H) * H`.
+        leftstart = copy(environments(below.left_gs, O.left_ham, below.left_gs).GLs[1]),
+        rightstart = copy(environments(below.right_gs, O.right_ham, below.right_gs).GRs[end])
+    )
+    return initialize_environments(below, O.finite_ham, above, leftstart, rightstart)
+end
 
 environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} =
     environments(below, nothing, above)

src/operators/mpohamiltonian.jl

@@ -707,7 +707,7 @@ end
 function Base.similar(H::MPOHamiltonian, ::Type{TorA}) where {TorA <: Union{Number, DenseVector}}
     return MPOHamiltonian(similar.(parent(H), TorA))
 end
-
+Base.circshift(H::InfiniteMPOHamiltonian, shift::Integer) = InfiniteMPOHamiltonian(circshift(parent(copy(H)), shift))
 # Linear Algebra
 # --------------
 function Base.:+(

src/operators/windowhamiltonian.jl

@@ -0,0 +1,86 @@
+# todo - what is the required interface for abstractmpo?
+# support densempo windows?
+"""
+$(TYPEDEF)
+
+The Hamiltonian counterpart of a [`WindowMPS`](@ref): a finite region embedded between an
+infinite environment to the left and to the right.
+It consists of an infinite Hamiltonian to the left, a finite Hamiltonian in the middle, and
+an infinite Hamiltonian to the right.
+
+Acts similar to just a finite Hamiltonian, but we "remember" the boundary Hamiltonians.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+## Constructors
+
+    WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
+
+Construct a `WindowMPOHamiltonian` by carving a finite `interval` out of an infinite
+Hamiltonian `ham`.
+The finite window consists of the sites in `interval`, while the left and right environments
+are copies of `ham` whose unit cells are circshifted so that they line up with the window
+boundaries.
+"""
+struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
+    "Hamiltonian acting on the infinite environment to the left of the window"
+    left_ham::InfiniteMPOHamiltonian{O}
+    "Hamiltonian acting on the finite window"
+    finite_ham::FiniteMPOHamiltonian{O}
+    "Hamiltonian acting on the infinite environment to the right of the window"
+    right_ham::InfiniteMPOHamiltonian{O}
+end
+
+function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
+    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))
+
+function Base.:+(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
+    # the finite window carries full Jordan virtual spaces at its boundaries, so it has to be
+    # block-diagonalized at every site (like an InfiniteMPOHamiltonian) rather than through
+    # the regular FiniteMPOHamiltonian addition
+    Ha, Hb = a.finite_ham, b.finite_ham
+    N = check_length(Ha, Hb)
+    finite_ham = similar(parent(Ha))
+    Vtriv = leftunitspace(first(physicalspace(Ha)))
+    for i in 1:N
+        A = cat(Ha[i].A, Hb[i].A; dims = (1, 4))
+        B = cat(Ha[i].B, Hb[i].B; dims = 1)
+        C = cat(Ha[i].C, Hb[i].C; dims = 3)
+        D = Ha[i].D + Hb[i].D
+
+        Vleft = ⊞(Vtriv, left_virtualspace(A), Vtriv)
+        Vright = ⊞(Vtriv, right_virtualspace(A), Vtriv)
+        V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright
+
+        finite_ham[i] = JordanMPOTensor(V, A, B, C, D)
+    end
+    return WindowMPOHamiltonian(
+        a.left_ham + b.left_ham, FiniteMPOHamiltonian(finite_ham), a.right_ham + b.right_ham
+    )
+end
+function Base.:*(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
+    return WindowMPOHamiltonian(
+        a.left_ham * b.left_ham, a.finite_ham * b.finite_ham, a.right_ham * b.right_ham
+    )
+end
+function VectorInterface.scale(H::WindowMPOHamiltonian, α::Number)
+    return WindowMPOHamiltonian(
+        scale(H.left_ham, α), scale(H.finite_ham, α), scale(H.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,34 @@ 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(ψ::InfiniteMPS, interval::UnitRange)
+
+Construct a [`WindowMPS`](@ref) from an infinite MPS `ψ` by promoting the sites in `interval`
+(or `1:L`) to the finite window while keeping `ψ` as the left and right infinite environments.
+The environment unit cells are circshifted so that they line up with the window boundaries.
+"""
+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/operators/windowhamiltonian.jl

@@ -0,0 +1,63 @@
+println("
+---------------------------------
+|   WindowMPOHamiltonian tests   |
+---------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+
+@testset "WindowMPOHamiltonian" begin
+    # a uniform state, represented with a two-site unit cell, so that carving out a window
+    # at different offsets exercises the (non-trivial) circshift alignment of the infinite
+    # environments while the physics stays translationally invariant.
+    H1 = force_planar(transverse_field_ising(; g = 1.5))
+    gs1, = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^12]), H1, VUMPS(; verbosity = 0))
+    gs = repeat(gs1, 2)
+    H = repeat(H1, 2)
+
+    L = 6
+
+    @testset "interval offset (circshift)" begin
+        # the window energy must be independent of the offset for a uniform state
+        energies = map((1:L, 2:(L + 1), 3:(L + 2))) do interval
+            ψ = WindowMPS(gs, interval)
+            Hw = WindowMPOHamiltonian(H, interval)
+            @test length(Hw) == length(interval)
+            return expectation_value(ψ, Hw)
+        end
+        @test energies[1] ≈ energies[2] atol = 1.0e-8
+        @test energies[1] ≈ energies[3] atol = 1.0e-8
+    end
+
+    @testset "linear algebra" begin
+        ψ = WindowMPS(gs, 1:L)
+        Hw = WindowMPOHamiltonian(H, 1:L)
+        e = expectation_value(ψ, Hw)
+
+        @test expectation_value(ψ, Hw + Hw) ≈ 2 * e atol = 1.0e-8
+        @test expectation_value(ψ, Hw - Hw) ≈ 0 atol = 1.0e-8
+        @test expectation_value(ψ, 3 * Hw) ≈ 3 * e atol = 1.0e-8
+        @test expectation_value(ψ, Hw * 3) ≈ 3 * e atol = 1.0e-8
+        @test expectation_value(ψ, -Hw) ≈ -e atol = 1.0e-8
+    end
+
+    @testset "distinct summands" begin
+        # add two genuinely different Hamiltonians: the ZZ Ising string has a single bulk
+        # level while the XY hopping needs two, so the boundary/bulk virtual spaces of the
+        # two summands differ and have to be block-diagonalized consistently.
+        ψ = WindowMPS(gs1, 1:L)
+        Ha = WindowMPOHamiltonian(H1, 1:L)
+        Hb = WindowMPOHamiltonian(force_planar(XY_model(; g = 0.7)), 1:L)
+        @test right_virtualspace(Ha.finite_ham, 3) != right_virtualspace(Hb.finite_ham, 3)
+
+        ea = expectation_value(ψ, Ha)
+        eb = expectation_value(ψ, Hb)
+        @test expectation_value(ψ, Ha + Hb) ≈ ea + eb atol = 1.0e-8
+        @test expectation_value(ψ, Hb + Ha) ≈ ea + eb atol = 1.0e-8
+        @test expectation_value(ψ, Ha - Hb) ≈ ea - eb 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)