diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 464097003..6d8f5ff89 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -123,6 +123,15 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   return value regardless, so this adds no extra cost. ([#509](https://github.com/QuantumKitHub/MPSKit.jl/pull/509))
 - `make_time_mpo` with `TaylorCluster` on a Hamiltonian whose virtual bond dimension varies along
   the chain are now correctly handled. ([#511](https://github.com/QuantumKitHub/MPSKit.jl/pull/511))
+- The converting constructor of `JordanMPO_AC_Hamiltonian` assigned the "finished" block `E` to
+  the "ending" field `B` whenever `B` was absent, raising a `convert` `MethodError` from deep
+  inside `AC_hamiltonian`. This is reached when an `MPOHamiltonian` whose scalartype differs from
+  the state's has an on-site term on a site where no interaction ends, such as a long-range one.
+  ([#493](https://github.com/QuantumKitHub/MPSKit.jl/pull/493))
+- Raised the `BlockTensorKit` compat lower bound to 0.3.19, which fixes a silent data-loss bug: with
+  `TensorKit` 0.17.2 and `BlockTensorKit` <= 0.3.18, a permuted contraction into a sparse block
+  tensor (e.g. an environment sweep near a `FiniteMPS` chain boundary when the operator and state
+  have different `scalartype`s) could silently drop data instead of erroring.
 
 ### Performance
 
diff --git a/test/groundstate/groundstate.jl b/test/groundstate/groundstate.jl
index 816c04809..b31b1512a 100644
--- a/test/groundstate/groundstate.jl
+++ b/test/groundstate/groundstate.jl
@@ -20,132 +20,152 @@ verbosity_conv = 1
     D = 6
     L = 10
 
-    H = force_planar(transverse_field_ising(; g, L))
+    models = [
+        "nearest-neighbour" => force_planar(transverse_field_ising(; g, L)),
+        "long-range, real scalartype" => force_planar(long_range_ising(Float64; g, L)),
+    ]
 
-    @testset "DMRG" begin
-        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^D)
-        v₀ = variance(ψ₀, H)
+    @testset "$name" for (name, H) in models
+        @testset "DMRG" begin
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^D)
+            v₀ = variance(ψ₀, H)
 
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2)
-        )
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2)
+            )
 
-        ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10), envs
-        )
-        v = variance(ψ, H)
+            ψ, envs, δ = find_groundstate(
+                ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10), envs
+            )
+            v = variance(ψ, H)
 
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀
-        @test v < 1.0e-2
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀
+            @test v < 1.0e-2
+        end
 
-        # the algorithm object carries no scratch space of its own - the sweep's allocator is
-        # obtained per solve - so re-using one across solves has to reproduce the answer
-        alg = DMRG(; verbosity = verbosity_conv, maxiter = 10)
-        ψ1, = find_groundstate(ψ₀, H, alg)
-        ψ2, = find_groundstate(ψ₀, H, alg)
-        @test expectation_value(ψ1, H) ≈ expectation_value(ψ2, H) atol = 1.0e-10
-    end
+        @testset "DMRG2" begin
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^D)
+            v₀ = variance(ψ₀, H)
+            trunc = truncrank(floor(Int, D * 1.5))
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, DMRG2(; verbosity = verbosity_full, maxiter = 2, trunc)
+            )
 
-    @testset "DMRG2" begin
-        ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D)
-        v₀ = variance(ψ₀, H)
-        trunc = truncrank(floor(Int, D * 1.5))
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG2(; verbosity = verbosity_full, maxiter = 2, trunc)
-        )
+            ψ, envs, δ = find_groundstate(
+                ψ, H, DMRG2(; verbosity = verbosity_conv, maxiter = 10, trunc), envs
+            )
+            v = variance(ψ, H)
 
-        ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG2(; verbosity = verbosity_conv, maxiter = 10, trunc), envs
-        )
-        v = variance(ψ, H)
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀
+            @test v < 1.0e-2
+        end
 
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀
-        @test v < 1.0e-2
-    end
+        @testset "CBEDMRG" begin
+            # start from a small bond so the bond expansion is exercised
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+            v₀ = variance(ψ₀, H)
+            expand = OptimalExpand(; trunc = truncrank(D ÷ 2))
+            trunc = truncrank(D)
 
-    @testset "CBEDMRG" begin
-        # start from a small bond so the bond expansion is exercised
-        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
-        v₀ = variance(ψ₀, H)
-        expand = OptimalExpand(; trunc = truncrank(D ÷ 2))
-        trunc = truncrank(D)
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc)
+            )
 
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc)
-        )
+            ψ, envs, δ = find_groundstate(
+                ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_expand = expand, trunc), envs
+            )
+            v = variance(ψ, H)
+
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀
+            @test v < 1.0e-2
+            # the bond should have grown to the truncation target
+            @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+        end
 
-        ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_expand = expand, trunc), envs
-        )
-        v = variance(ψ, H)
+        @testset "CBEDMRG (SketchedExpand)" begin
+            # randomized bond expansion at single-site cost. The sketch is redrawn every sweep, so an
+            # aggressive expansion (a large fraction of the bond) keeps the single-site Galerkin error
+            # noisy; a gentle per-sweep increment lets it converge like the deterministic expanders.
+            Random.seed!(1234)
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+            v₀ = variance(ψ₀, H)
+            expand = SketchedExpand(; trunc = truncrank(2), oversampling = 4)
+            trunc = truncrank(D)
 
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀
-        @test v < 1.0e-2
-        # the bond should have grown to the truncation target
-        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
-    end
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc)
+            )
 
-    @testset "CBEDMRG (SketchedExpand)" begin
-        # randomized bond expansion at single-site cost. The sketch is redrawn every sweep, so an
-        # aggressive expansion (a large fraction of the bond) keeps the single-site Galerkin error
-        # noisy; a gentle per-sweep increment lets it converge like the deterministic expanders.
-        Random.seed!(1234)
-        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
-        v₀ = variance(ψ₀, H)
-        expand = SketchedExpand(; trunc = truncrank(2), oversampling = 4)
-        trunc = truncrank(D)
+            ψ, envs, δ = find_groundstate(
+                ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trunc), envs
+            )
+            v = variance(ψ, H)
+
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀
+            @test v < 1.0e-2
+            # the bond should have grown to the truncation target
+            @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+        end
 
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc)
-        )
+        @testset "DMRG3S" begin
+            # start from a small bond so the post-expansion is exercised, mirroring CBEDMRG above
+            Random.seed!(1234)
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+            v₀ = variance(ψ₀, H)
+            alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7))  # TODO: match final constructor API
+            trunc = truncrank(D)
 
-        ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trunc), envs
-        )
-        v = variance(ψ, H)
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_gauge, trunc)
+            )
 
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀
-        @test v < 1.0e-2
-        # the bond should have grown to the truncation target
-        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
-    end
+            ψ, envs, δ = find_groundstate(
+                ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_gauge, trunc), envs
+            )
+            v = variance(ψ, H)
+
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀
+            @test v < 1.0e-2
+            # the bond should have grown to the truncation target
+            @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+        end
 
-    @testset "DMRG3S" begin
-        # start from a small bond so the post-expansion is exercised, mirroring CBEDMRG above
-        Random.seed!(1234)
-        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
-        v₀ = variance(ψ₀, H)
-        alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7))  # TODO: match final constructor API
-        trunc = truncrank(D)
+        @testset "GradientGrassmann" begin
+            ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^D)
+            v₀ = variance(ψ₀, H)
 
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_gauge, trunc)
-        )
+            # test logging
+            ψ, envs, δ = find_groundstate(
+                ψ₀, H, GradientGrassmann(; verbosity = verbosity_full, maxiter = 2)
+            )
 
-        ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_gauge, trunc), envs
-        )
-        v = variance(ψ, H)
+            # an explicit `tol` keeps the optimizer from overshooting past the point where the
+            # gradient is floating-point noise: pushed further, the CG line search can hit a
+            # 0/0 in its step-size formula and feed a NaN tangent into the Grassmann retraction
+            ψ, envs, δ = find_groundstate(
+                ψ, H, GradientGrassmann(; tol, verbosity = verbosity_conv, maxiter = 50), envs
+            )
+            v = variance(ψ, H)
 
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀
-        @test v < 1.0e-2
-        # the bond should have grown to the truncation target
-        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+            # test using low variance
+            @test sum(δ) ≈ 0 atol = 1.0e-3
+            @test v < v₀ && v < 1.0e-2
+        end
     end
 
     fast_tests || @testset "DMRG3S escapes local minimum (Hubig et al. 2015, Sec. VII A)" begin
@@ -173,25 +193,6 @@ verbosity_conv = 1
         @test E_escape < E_stuck - 1.0
         @test isapprox(E_escape, -8.6824724; atol = 1.0e-4)
     end
-
-    @testset "GradientGrassmann" begin
-        ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D)
-        v₀ = variance(ψ₀, H)
-
-        # test logging
-        ψ, envs, δ = find_groundstate(
-            ψ₀, H, GradientGrassmann(; verbosity = verbosity_full, maxiter = 2)
-        )
-
-        ψ, envs, δ = find_groundstate(
-            ψ, H, GradientGrassmann(; verbosity = verbosity_conv, maxiter = 50), envs
-        )
-        v = variance(ψ, H)
-
-        # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
-        @test v < v₀ && v < 1.0e-2
-    end
 end
 
 @testset "InfiniteMPS ground state" verbose = true begin
@@ -206,6 +207,7 @@ end
     H_ref_realT = force_planar(transverse_field_ising(Float64; g))
     @test variance(ψ, H_ref_realT) ≈ v₀ atol = 1.0e-10
 
+
     @testset "VUMPS (unit cell $unit_cell_size, $schedname)" for unit_cell_size in [1, 3],
             (schedname, scheduler) in SCHEDULERS
 
@@ -227,6 +229,20 @@ end
         @test v < 1.0e-2
     end
 
+    @testset "VUMPS (long-range, real scalartype)" begin
+        H = force_planar(long_range_ising_infinite(Float64; g, L = 3))
+        ψ₀ = InfiniteMPS(fill(ℙ^2, 3), fill(ℙ^D, 3))
+        v₀ = variance(ψ₀, H)
+
+        ψ′, envs, δ = find_groundstate(ψ₀, H, VUMPS(; tol, verbosity = verbosity_conv, maxiter = 20))
+        v = variance(ψ′, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test v < v₀
+        @test v < 1.0e-2
+    end
+
     @testset "IDMRG" for unit_cell_size in [1, 3]
         ψ = unit_cell_size == 1 ? InfiniteMPS(ℙ^2, ℙ^D) : repeat(ψ, unit_cell_size)
         H = repeat(H_ref, unit_cell_size)
@@ -245,6 +261,20 @@ end
         @test v < 1.0e-2
     end
 
+    @testset "IDMRG (long-range, real scalartype)" begin
+        H = force_planar(long_range_ising_infinite(Float64; g, L = 3))
+        ψ₀ = InfiniteMPS(fill(ℙ^2, 3), fill(ℙ^D, 3))
+        v₀ = variance(ψ₀, H)
+
+        ψ, envs, δ = find_groundstate(ψ₀, H, IDMRG(; tol, verbosity = verbosity_conv, maxiter = 20))
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test v < v₀
+        @test v < 1.0e-2
+    end
+
     @testset "IDMRG2" begin
         ψ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 2)
         H = repeat(H_ref, 2)
diff --git a/test/hamiltonian/derivatives.jl b/test/hamiltonian/derivatives.jl
new file mode 100644
index 000000000..246637444
--- /dev/null
+++ b/test/hamiltonian/derivatives.jl
@@ -0,0 +1,122 @@
+println("
+----------------------------
+|   Derivative operators    |
+----------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using MPSKit: C_hamiltonian, AC_hamiltonian, AC2_hamiltonian
+using MPSKit: _transpose_front, _transpose_tail
+using TensorKit
+using BlockTensorKit: nonzero_length
+using Random
+
+Random.seed!(1234)
+
+# `JordanMPO_AC(2)_Hamiltonian` has two construction paths: a fast one, taken when every
+# block already has the storagetype of the environments, and a converting one
+# this converting path is reached when `scalartype(H) != scalartype(ψ)` and an on-site D block is passed
+# nearest-neighbour models have an ending B block on every site (but the first)
+
+@testset "Jordan block structure of the long-range models" begin
+    H = long_range_ising(Float64; L = 4)
+    Hi = long_range_ising_infinite(Float64; L = 3)
+
+    # nothing ends before the far end of the (long-range) bond
+    @test all(i -> nonzero_length(H[i].B) == 0, 1:3)
+    @test nonzero_length(H[4].B) == 1
+    @test all(i -> nonzero_length(Hi[i].B) == 0, 1:2)
+    @test nonzero_length(Hi[3].B) == 1
+
+    # while the left virtual space is non-trivial everywhere but on the first site of
+    # the finite chain, so `E` is present exactly where `B` is not
+    @test size(H[1], 1) == 1
+    @test all(i -> size(H[i], 1) > 1, 2:4)
+    @test all(i -> size(Hi[i], 1) > 1, 1:3)
+
+    # every site carries an on-site term, which is what forces the converting constructor
+    @test all(i -> nonzero_length(H[i].D) == 1, 1:4)
+    @test all(i -> nonzero_length(Hi[i].D) == 1, 1:3)
+end
+
+@testset "MPOHamiltonian derivatives: real operator, complex state" verbose = true begin
+    D = 8
+    L_inf, L_fin = 3, 4
+    ψ_inf = InfiniteMPS(randn, ComplexF64, fill(ℂ^2, L_inf), fill(ℂ^D, L_inf))
+    ψ_fin = FiniteMPS(randn, ComplexF64, L_fin, ℂ^2, ℂ^D)
+    models = [
+        "FiniteMPS, long-range" => (long_range_ising(Float64; L = L_fin), ψ_fin),
+        "FiniteMPS, nearest-neighbour" => (transverse_field_ising(Float64; g = 4.0, L = L_fin), ψ_fin),
+        "InfiniteMPS, long-range" => (long_range_ising_infinite(Float64; L = L_inf), ψ_inf),
+        "InfiniteMPS, nearest-neighbour" => (repeat(transverse_field_ising(Float64; g = 4.0), 3), ψ_inf),
+    ]
+
+    @testset "$name" for (name, (H, ψ)) in models
+        @test scalartype(H) <: Real
+        @test scalartype(ψ) <: Complex
+
+        # `complex(H)` takes the fast path, `H` the converting one: both must agree
+        Hc = complex(H)
+        @test scalartype(Hc) <: Complex
+        envs = environments(ψ, H, ψ)
+        envs_c = environments(ψ, Hc, ψ)
+        @test expectation_value(ψ, H, envs) ≈ expectation_value(ψ, Hc, envs_c)
+
+        L = length(ψ)
+        bonds = isfinite(H) ? (1:(L - 1)) : (1:L)
+
+        @testset "AC" begin
+            for site in 1:L
+                AC = ψ.AC[site]
+                @test AC_hamiltonian(site, ψ, H, ψ, envs)(AC) ≈
+                    AC_hamiltonian(site, ψ, Hc, ψ, envs_c)(AC)
+            end
+        end
+
+        @testset "AC2" begin
+            for site in bonds
+                AC2 = _transpose_front(ψ.AC[site]) * _transpose_tail(ψ.AR[site + 1])
+                @test AC2_hamiltonian(site, ψ, H, ψ, envs)(AC2) ≈
+                    AC2_hamiltonian(site, ψ, Hc, ψ, envs_c)(AC2)
+            end
+        end
+
+        @testset "C" begin
+            for site in bonds
+                C = ψ.C[site]
+                @test C_hamiltonian(site, ψ, H, ψ, envs)(C) ≈
+                    C_hamiltonian(site, ψ, Hc, ψ, envs_c)(C)
+            end
+        end
+    end
+end
+
+@testset "FiniteMPS derivatives reproduce the expectation value" verbose = true begin
+    # in mixed gauge the derivative operators contract everything but the center site(s),
+    # so their expectation value is the full energy, independent of where we sit
+    D = 8
+    L = 10
+    models = [
+        "long-range" => long_range_ising(Float64; L = L),
+        "nearest-neighbour" => transverse_field_ising(Float64; g = 4.0, L = L),
+    ]
+
+    @testset "$name" for (name, H) in models
+        for T in (Float64, ComplexF64)
+            ψ = normalize!(FiniteMPS(randn, T, L, ℂ^2, ℂ^D))
+            envs = environments(ψ, H, ψ)
+            E = expectation_value(ψ, H, envs)
+
+            for site in 1:L
+                AC = ψ.AC[site]
+                @test dot(AC, AC_hamiltonian(site, ψ, H, ψ, envs)(AC)) ≈ E
+            end
+            for site in 1:(L - 1)
+                AC2 = _transpose_front(ψ.AC[site]) * _transpose_tail(ψ.AR[site + 1])
+                @test dot(AC2, AC2_hamiltonian(site, ψ, H, ψ, envs)(AC2)) ≈ E
+            end
+        end
+    end
+end
diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl
index 86b82bd2d..107675c05 100644
--- a/test/setup/testsetup.jl
+++ b/test/setup/testsetup.jl
@@ -21,6 +21,7 @@ export f_plus_f_min, f_min_f_plus, f_num, f_hopping
 export force_planar
 export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model, kitaev_model
+export long_range_ising, long_range_ising_infinite
 export classical_ising_tensors, classical_ising, sixvertex
 export bad_initial_state
 export SCHEDULERS, with_scheduler
@@ -204,6 +205,23 @@ function kitaev_model(
     end
 end
 
+function long_range_ising(
+        T::Type{<:Number} = Float64, sym::Type{<:Sector} = Trivial; g = 4.0, L = 4
+    )
+    X = S_x(T, sym; spin = 1 // 2) * 2
+    ZZ = S_z_S_z(T, sym; spin = 1 // 2) * 4
+    lattice = fill(space(X, 1), L)
+    return FiniteMPOHamiltonian(lattice, ((i,) => -g * X for i in 1:L)..., (1, L) => -ZZ)
+end
+function long_range_ising_infinite(
+        T::Type{<:Number} = Float64, sym::Type{<:Sector} = Trivial; g = 4.0, L = 3
+    )
+    X = S_x(T, sym; spin = 1 // 2) * 2
+    ZZ = S_z_S_z(T, sym; spin = 1 // 2) * 4
+    lattice = PeriodicArray(fill(space(X, 1), L))
+    return InfiniteMPOHamiltonian(lattice, ((i,) => -g * X for i in 1:L)..., (1, L) => -ZZ)
+end
+
 function ising_bond_tensor(β)
     J = 1.0
     K = β * J