diff --git a/src/MPSKit.jl b/src/MPSKit.jl index b29174216..694dd8024 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 -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") diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index c147ba391..d60213af2 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 de04e8078..93bfd6637 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/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) diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index fee1555bb..cb11cd53f 100644 --- a/src/environments/finite_envs.jl +++ b/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) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 37c52b9fb..f89d46a9f 100644 --- a/src/operators/mpohamiltonian.jl +++ b/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.:+( diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl new file mode 100644 index 000000000..f4ace7c45 --- /dev/null +++ b/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) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 5c0bab9d0..fb2396da6 100644 --- a/src/states/windowmps.jl +++ b/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 #=========================================================================================== 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/operators/windowhamiltonian.jl b/test/operators/windowhamiltonian.jl new file mode 100644 index 000000000..46d504e2a --- /dev/null +++ b/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 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)