Three code fixes carved out of the docs pass (#449)

docs/src/changelog.md

@@ -51,13 +51,18 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   Accordingly, `marek_gap` and `correlation_length` now return a `TensorKit.SectorDict` of
   per-sector results by default; pass `sector = ...` to obtain a single sector's result as before.
 - All `trscheme` keyword arguments are renamed to `trunc` ([#482](https://github.com/QuantumKitHub/MPSKit.jl/pull/482)).
+- `correlator` now throws an `ArgumentError` when the sites are not ordered as `i < j`.
+  Previously such a call only logged an `@error` and then continued into a contraction that is
+  not the requested correlator. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
 
 ### Deprecated
 
 ### Removed
 
 ### Fixed
 
+- `isfinite(::WindowMPOHamiltonian)` was undefined. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
+- `excitations(::InfiniteMPO, ::QuasiparticleAnsatz, ::InfiniteQP, lenvs, renvs)` referenced `H_eff`  before assigning. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
 - `Base.:+`/`-` on `FiniteMPS` returned a wrong state for near-parallel operands carried by
   different tensor networks, e.g. `norm(E₀ * gs - H * gs)` coming out as `2 * norm(gs) * E₀`
   instead of ~0. The lazy gauge sweep in `CView.getindex` re-derived `AL`/`C` entries that were

src/algorithms/correlators.jl

@@ -4,6 +4,7 @@
 
 Compute the 2-point correlator <ψ|O1[i]O2[j]|ψ> for inserting `O1` at `i` and `O2` at `j`.
 Also accepts ranges for `j`.
+The sites must be ordered as `i < j`; other orderings throw an `ArgumentError`.
 """
 function correlator end
 
@@ -14,7 +15,7 @@ end
 function correlator(
         state::AbstractMPS, O₁::MPOTensor, O₂::MPOTensor, i::Int, js::AbstractRange{Int}
     )
-    first(js) > i || @error "i should be smaller than j ($i, $(first(js)))"
+    first(js) > i || throw(ArgumentError("i should be smaller than j ($i, $(first(js)))"))
     S₁ = _firstspace(O₁)
     isunitspace(S₁) || throw(ArgumentError("O₁ should start with a trivial leg."))
     S₂ = _lastspace(O₂)

src/algorithms/excitation/quasiparticleexcitation.jl

@@ -214,7 +214,7 @@ function excitations(
         num = 1
     )
     E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs)
-    H_eff = EffectiveExcitationHamiltonian(H_eff, lenvs, renvs, E)
+    H_eff = EffectiveExcitationHamiltonian(H, lenvs, renvs, E)
 
     Es, ϕs, convhist = eigsolve(ϕ₀, num, :LM, alg.alg) do ϕ
         return H_eff(ϕ, alg.alg_environments)

src/operators/windowhamiltonian.jl

@@ -32,6 +32,7 @@ struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
     "Hamiltonian acting on the infinite environment to the right of the window"
     right_ham::InfiniteMPOHamiltonian{O}
 end
+Base.isfinite(::Type{<:WindowMPOHamiltonian}) = true
 
 function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
     left_edge = (interval.start - 1) % length(ham)

test/algorithms/correlators.jl

@@ -41,4 +41,12 @@ end
     @test isapprox(last(G), last(G2), atol = 1.0e-2)
     @test isapprox(G[1], expectation_value(ψ, (1, 2) => S_z_S_z()), atol = 1.0e-2)
     @test isapprox(G[2], expectation_value(ψ, (1, 3) => S_z_S_z()), atol = 1.0e-2)
+
+    # the sites have to be ordered `i < j`: every `correlator` method funnels its
+    # arguments through the range version, so all three signatures have to throw
+    @test_throws ArgumentError correlator(ψ, Z_mpo, Z_mpo, 3, 2)
+    @test_throws ArgumentError correlator(ψ, Z_mpo, Z_mpo, 3, 1:2)
+    @test_throws ArgumentError correlator(ψ, S_z_S_z(), 3, 2)
+    # equal sites are not allowed either
+    @test_throws ArgumentError correlator(ψ, Z_mpo, Z_mpo, 2, 2)
 end

test/algorithms/excitations.jl

@@ -52,6 +52,35 @@ verbosity_conv = 1
         @test abs(energies[1]) > abs(energies[2]) # has a minimum at pi/2
     end
 
+    @testset "infinite (mpo) - quasiparticle input" begin
+        # regression test: handing an `InfiniteQP` to an `InfiniteMPO`, rather than a
+        # momentum, used to build the effective Hamiltonian out of the (still undefined)
+        # `H_eff` instead of `H`, so every call through this entry point threw an
+        # `UndefVarError`. The momentum entry point is unaffected, as it converts to a
+        # `MultilineMPO` first. The classical Ising transfer matrix is used here (rather than
+        # the six-vertex model above) because its quasiparticle eigenproblem is well enough
+        # conditioned to compare the two entry points at tight tolerance.
+        H = classical_ising(; β = 0.4)
+        ψ = InfiniteMPS([ℂ^2], [ℂ^12])
+        ψ, envs, = leading_boundary(ψ, H, VUMPS(; maxiter = 400, verbosity = 0, tol = 1.0e-12))
+
+        alg = QuasiparticleAnsatz(; tol = 1.0e-10)
+        for momentum in (0.0, Float64(pi / 2))
+            ϕ₀ = MPSKit.LeftGaugedQP(rand, ψ, ψ; momentum)
+            Es, ϕs = excitations(H, alg, ϕ₀, envs, envs; num = 1)
+
+            # same spectrum as the momentum entry point, which routes through MultilineMPO
+            Es_momentum, = excitations(H, alg, momentum, ψ, envs; num = 1)
+            @test Es[1] ≈ Es_momentum[1] atol = 1.0e-8
+
+            # the operator really is `H`: the result is an eigenpair of the corresponding
+            # effective excitation Hamiltonian
+            E = MPSKit.effective_excitation_renormalization_energy(H, ϕs[1], envs, envs)
+            H_eff = MPSKit.EffectiveExcitationHamiltonian(H, envs, envs, E)
+            @test norm(H_eff(ϕs[1], alg.alg_environments) - Es[1] * ϕs[1]) < 1.0e-8
+        end
+    end
+
     @testset "finite" begin
         verbosity = verbosity_conv
         H_inf = force_planar(transverse_field_ising())

test/operators/windowhamiltonian.jl

@@ -27,6 +27,10 @@ using TensorKit: ℙ
             ψ = WindowMPS(gs, interval)
             Hw = WindowMPOHamiltonian(H, interval)
             @test length(Hw) == length(interval)
+            # the window itself is finite, even though its environments are not
+            @test isfinite(Hw)
+            @test isfinite(typeof(Hw))
+            @test isfinite(Hw) == isfinite(typeof(Hw))
             return expectation_value(ψ, Hw)
         end
         @test energies[1] ≈ energies[2] atol = 1.0e-8