From c50d042e4544a2151e588b94186fd23545ced06b Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 24 Jun 2026 20:36:26 -0400
Subject: [PATCH 1/9] try and update the JordanMPOTensor structure

---
 src/operators/jordanmpotensor.jl | 43 ++++++++++----------------------
 1 file changed, 13 insertions(+), 30 deletions(-)

diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index af95282cc..20b03f54f 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -12,38 +12,21 @@ A tensor map that represents the upper triangular block form of a matrix product
 ```
 """
 struct JordanMPOTensor{
-        E, S,
-        TA <: AbstractTensorMap{E, S, 2, 2},
-        TB <: AbstractTensorMap{E, S, 2, 1},
-        TC <: AbstractTensorMap{E, S, 1, 2},
-        TD <: AbstractTensorMap{E, S, 1, 1},
-    } <: AbstractBlockTensorMap{E, S, 2, 2}
-    V::TensorMapSumSpace{S, 2, 2}
-    A::SparseBlockTensorMap{TA, E, S, 2, 2, 4}
-    B::SparseBlockTensorMap{TB, E, S, 2, 1, 3}
-    C::SparseBlockTensorMap{TC, E, S, 1, 2, 3}
-    D::SparseBlockTensorMap{TD, E, S, 1, 1, 2}
+        T <: Number, S, A <: DenseVector{T},
+    } <: AbstractBlockTensorMap{T, S, 2, 2}
+    tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4}
+    scalars::Dict{CartesianIndex{2}, T}
+
     # uninitialized constructor
-    function JordanMPOTensor{E, S, TA, TB, TC, TD}(
+    function JordanMPOTensor{T, S, A}(
             ::UndefInitializer, V::TensorMapSumSpace{S, 2, 2}
-        ) where {E, S, TA, TB, TC, TD}
-        allVs = eachspace(V)
-
-        # Note that this is a bit of a hack using end to get the last index:
-        # it should be 1 or end depending on this being an "edge" tensor or a "bulk" tensor
-        VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)])
-        A = SparseBlockTensorMap{TA}(undef, VA)
-
-        VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4)
-        B = SparseBlockTensorMap{TB}(undef, VB)
-
-        VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1)
-        C = SparseBlockTensorMap{TC}(undef, VC)
-
-        VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1)
-        D = SparseBlockTensorMap{TD}(undef, VD)
-
-        return new{E, S, TA, TB, TC, TD}(V, A, B, C, D)
+        ) where {T, S, A}
+        tensors = SparseBlockTensorMap{tensormaptype(S, 2, 2, A)}(undef, V)
+        scalars = Dict{CartesianIndex{2}, T}()
+        rows, cols = size(tensors, 1), size(tensors, 4)
+        cols > 1 && (scalars[CartesianIndex(1, 1)] = one(T))
+        rows > 1 && (scalars[CartesianIndex(rows, cols)] = one(T))
+        return new{T, S, A}(tensors, scalars)
     end
 
     # constructor from data

From d57b8e6949c8e19e001018d994f88c15899eaa12 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 24 Jun 2026 21:07:07 -0400
Subject: [PATCH 2/9] Rework JordanMPOTensor internals over tensors/scalars

Drive JordanMPOTensor off a single full-space `tensors` SparseBlockTensorMap
plus a `scalars` Dict of identity coefficients, instead of the four reduced-leg
A/B/C/D fields. The four-block decomposition is preserved as derived
`getproperty` accessors (with a 5-arg constructor that writes back), so existing
consumers keep working. Generic `getindex` materializes scalar entries as
`TensorMap`s; the `.A` accessor preserves `BraidingTensor`s for hot kernels.

Drop the now-redundant `JordanMPOTensorMap` alias and its export.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/MPSKit.jl                    |   2 +-
 src/operators/jordanmpotensor.jl | 288 +++++++++++++++----------------
 2 files changed, 141 insertions(+), 149 deletions(-)

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index c39977731..a4e1fe605 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -18,7 +18,7 @@ export QP, LeftGaugedQP, RightGaugedQP
 # operators:
 export AbstractMPO
 export MPO, FiniteMPO, InfiniteMPO
-export JordanMPOTensor, JordanMPOTensorMap
+export JordanMPOTensor
 export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian
 export MultilineMPO
 export UntimedOperator, TimedOperator, MultipliedOperator, LazySum
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index 20b03f54f..b47efa2af 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -1,5 +1,5 @@
 """
-    JordanMPOTensor{E,S,TA,TB,TC,TD} <: AbstractBlockTensorMap{E,S,2,2}
+    JordanMPOTensor{T,S,A} <: AbstractBlockTensorMap{T,S,2,2}
 
 A tensor map that represents the upper triangular block form of a matrix product operator (MPO).
 
@@ -10,6 +10,10 @@ A tensor map that represents the upper triangular block form of a matrix product
 0 & 0 & 1
 \\end{pmatrix}
 ```
+
+The genuine operators are stored in `tensors` over the full virtual space, while the
+scalar multiples of the identity (in particular the diagonal `1`s) are kept separately in
+`scalars`, keyed by their `(row, col)` virtual position.
 """
 struct JordanMPOTensor{
         T <: Number, S, A <: DenseVector{T},
@@ -17,6 +21,14 @@ struct JordanMPOTensor{
     tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4}
     scalars::Dict{CartesianIndex{2}, T}
 
+    # constructor from fields
+    function JordanMPOTensor{T, S, A}(
+            tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4},
+            scalars::Dict{CartesianIndex{2}, T}
+        ) where {T, S, A}
+        return new{T, S, A}(tensors, scalars)
+    end
+
     # uninitialized constructor
     function JordanMPOTensor{T, S, A}(
             ::UndefInitializer, V::TensorMapSumSpace{S, 2, 2}
@@ -28,27 +40,8 @@ struct JordanMPOTensor{
         rows > 1 && (scalars[CartesianIndex(rows, cols)] = one(T))
         return new{T, S, A}(tensors, scalars)
     end
-
-    # constructor from data
-    function JordanMPOTensor{E, S, TA, TB, TC, TD}(
-            V::TensorMapSumSpace,
-            A::SparseBlockTensorMap{<:Any, <:Any, S, 2, 2},
-            B::SparseBlockTensorMap{<:Any, <:Any, S, 2, 1},
-            C::SparseBlockTensorMap{<:Any, <:Any, S, 1, 2},
-            D::SparseBlockTensorMap{<:Any, <:Any, S, 1, 1}
-        ) where {E, S, TA, TB, TC, TD}
-        return new{E, S, TA, TB, TC, TD}(V, A, B, C, D)
-    end
 end
 
-const JordanMPOTensorMap{T, S, A <: DenseVector{T}} = JordanMPOTensor{
-    T, S,
-    Union{TensorMap{T, S, 2, 2, A}, BraidingTensor{T, S, A}},
-    TensorMap{T, S, 2, 1, A},
-    TensorMap{T, S, 1, 2, A},
-    TensorMap{T, S, 1, 1, A},
-}
-
 function JordanMPOTensor{E, S}(::UndefInitializer, V::TensorMapSumSpace{S}) where {E, S}
     return jordanmpotensortype(S, E)(undef, V)
 end
@@ -74,7 +67,21 @@ function JordanMPOTensor(
     VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1)
     VD == space(D) || throw(SpaceMismatch(lazy"D-block has incompatible spaces:\n$VD\n$(space(D))"))
 
-    return JordanMPOTensor{E, S, TA, TB, TC, TD}(V, A, B, C, D)
+    W = jordanmpotensortype(S, storagetype(A))(undef, V)
+    cols = size(W, 4)
+    for (I, v) in nonzero_pairs(A)
+        W[I[1] + 1, 1, 1, I[4] + 1] = v
+    end
+    for (I, v) in nonzero_pairs(B)
+        W[I[1] + 1, 1, 1, cols] = insertrightunit(v, 3)
+    end
+    for (I, v) in nonzero_pairs(C)
+        W[1, 1, 1, I[3] + 1] = insertleftunit(v, 1)
+    end
+    if nonzero_length(D) > 0
+        W[1, 1, 1, cols] = insertrightunit(insertleftunit(only(D), 1), 3)
+    end
+    return W
 end
 function JordanMPOTensor(
         V::TensorMapSumSpace{S, 2, 2},
@@ -105,13 +112,8 @@ function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E,
 end
 
 function jordanmpotensortype(::Type{S}, ::Type{E}) where {S <: VectorSpace, E}
-    TA = tensormaptype(S, 2, 2, E)
-    T = scalartype(TA)
-    Tτ = BraidingTensor{T, S, storagetype(TA)}
-    TB = tensormaptype(S, 2, 1, E)
-    TC = tensormaptype(S, 1, 2, E)
-    TD = tensormaptype(S, 1, 1, E)
-    return JordanMPOTensor{T, S, Union{TA, Tτ}, TB, TC, TD}
+    TT = tensormaptype(S, 2, 2, E)
+    return JordanMPOTensor{scalartype(TT), S, storagetype(TT)}
 end
 function jordanmpotensortype(::Type{O}) where {O <: AbstractTensorMap}
     return jordanmpotensortype(spacetype(O), storagetype(O))
@@ -123,31 +125,77 @@ end
 
 # Properties
 # ----------
-TensorKit.space(W::JordanMPOTensor) = W.V
-Base.eltype(::Type{JordanMPOTensor{E, S, TA, TB, TC, TD}}) where {E, S, TA, TB, TC, TD} = TA
+TensorKit.space(W::JordanMPOTensor) = space(getfield(W, :tensors))
+function Base.eltype(::Type{<:JordanMPOTensor{T, S, A}}) where {T, S, A}
+    return tensormaptype(S, 2, 2, A)
+end
+
+function Base.getproperty(W::JordanMPOTensor, sym::Symbol)
+    sym === :A && return _jordan_A(W)
+    sym === :B && return _jordan_B(W)
+    sym === :C && return _jordan_C(W)
+    sym === :D && return _jordan_D(W)
+    return getfield(W, sym)
+end
+
+# reduced-leg blocks reconstructed from `tensors` and `scalars`
+function _jordan_A(W::JordanMPOTensor)
+    rows, cols = size(W, 1), size(W, 4)
+    Tτ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}
+    TA = tensormaptype(spacetype(W), 2, 2, storagetype(W))
+    tensors = getfield(W, :tensors)
+    VA = space(eachspace(tensors)[2:(end - 1), 1, 1, 2:(end - 1)])
+    A = SparseBlockTensorMap{Union{TA, Tτ}}(undef, VA)
+    for (I, v) in nonzero_pairs(tensors)
+        (1 < I[1] < rows && 1 < I[4] < cols) || continue
+        A[I[1] - 1, 1, 1, I[4] - 1] = v
+    end
+    for (K, c) in getfield(W, :scalars)
+        i, j = K[1], K[2]
+        (1 < i < rows && 1 < j < cols) || continue
+        τ = Tτ(eachspace(A)[i - 1, 1, 1, j - 1])
+        A[i - 1, 1, 1, j - 1] = isone(c) ? τ : scale!(TensorMap(τ), c)
+    end
+    return A
+end
+function _jordan_B(W::JordanMPOTensor)
+    rows, cols = size(W, 1), size(W, 4)
+    TB = tensormaptype(spacetype(W), 2, 1, storagetype(W))
+    VB = removeunit(space(eachspace(getfield(W, :tensors))[2:(end - 1), 1, 1, end]), 4)
+    B = SparseBlockTensorMap{TB}(undef, VB)
+    for I in nonzero_keys(W)
+        (1 < I[1] < rows && I[4] == cols) || continue
+        B[I[1] - 1, 1, 1] = removeunit(W[I], 4)
+    end
+    return B
+end
+function _jordan_C(W::JordanMPOTensor)
+    cols = size(W, 4)
+    TC = tensormaptype(spacetype(W), 1, 2, storagetype(W))
+    VC = removeunit(space(eachspace(getfield(W, :tensors))[1, 1, 1, 2:(end - 1)]), 1)
+    C = SparseBlockTensorMap{TC}(undef, VC)
+    for I in nonzero_keys(W)
+        (I[1] == 1 && 1 < I[4] < cols) || continue
+        C[1, 1, I[4] - 1] = removeunit(W[I], 1)
+    end
+    return C
+end
+function _jordan_D(W::JordanMPOTensor)
+    cols = size(W, 4)
+    TD = tensormaptype(spacetype(W), 1, 1, storagetype(W))
+    VD = removeunit(removeunit(space(eachspace(getfield(W, :tensors))[1, 1, 1, end:end]), 4), 1)
+    D = SparseBlockTensorMap{TD}(undef, VD)
+    K = CartesianIndex(1, 1, 1, cols)
+    if haskey(W, K)
+        D[1, 1] = removeunit(removeunit(W[K], 4), 1)
+    end
+    return D
+end
 
 function Base.haskey(W::JordanMPOTensor, I::CartesianIndex{4})
     Base.checkbounds(W, I.I...)
-    # only has braiding tensors if sizes are large enough
-    sz = size(W)
-    (
-        sz[1] > 1 && I == CartesianIndex(1, 1, 1, 1) ||
-            sz[4] > 1 && I == CartesianIndex(sz[1], 1, 1, sz[4])
-    ) && return true
-
-    row, col = I.I[1], I.I[4]
-
-    if row == 1 && col == sz[4]
-        return haskey(W.D, CartesianIndex(1, 1))
-    elseif row == 1
-        return haskey(W.C, CartesianIndex(1, 1, col - 1))
-    elseif col == sz[4]
-        return haskey(W.B, CartesianIndex(row - 1, 1, 1))
-    elseif 1 < row < sz[1] && 1 < col < sz[4]
-        return haskey(W.A, CartesianIndex(row - 1, 1, 1, col - 1))
-    else
-        return false
-    end
+    return haskey(getfield(W, :scalars), CartesianIndex(I[1], I[4])) ||
+        haskey(getfield(W, :tensors), I)
 end
 
 Base.parent(W::JordanMPOTensor) = parent(SparseBlockTensorMap(W))
@@ -157,51 +205,30 @@ BlockTensorKit.issparse(W::JordanMPOTensor) = true
 # Converters
 # ----------
 function BlockTensorKit.SparseBlockTensorMap(W::JordanMPOTensor)
-    τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[1])
     W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W), spacetype(W), 2, 2}}(
         undef_blocks, space(W)
     )
-    if size(W, 4) > 1
-        W′[1, 1, 1, 1] = τ
-    end
-    if size(W, 1) > 1
-        W′[end, 1, 1, end] = τ
-    end
-
-    Ia = CartesianIndex(1, 0, 0, 1)
-    for (I, v) in nonzero_pairs(W.A)
-        W′[I + Ia] = v
-    end
-    Ib = CartesianIndex(1, 0, 0)
-    for (I, v) in nonzero_pairs(W.B)
-        W′[I + Ib, size(W′, 4)] = insertrightunit(v, 3)
+    for (I, v) in nonzero_pairs(getfield(W, :tensors))
+        W′[I] = v
     end
-    Ic = CartesianIndex(0, 0, 1)
-    for (I, v) in nonzero_pairs(W.C)
-        W′[1, I + Ic] = insertleftunit(v, 1)
+    for (K, c) in getfield(W, :scalars)
+        τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(
+            eachspace(W)[K[1], 1, 1, K[2]]
+        )
+        W′[K[1], 1, 1, K[2]] = isone(c) ? τ : scale!(TensorMap(τ), c)
     end
-    if nonzero_length(W.D) > 0
-        W′[1, 1, 1, end] = insertrightunit(insertleftunit(only(W.D), 1), 3)
-    end
-
     return W′
 end
 
 for f in (:real, :complex)
     @eval function Base.$f(W::JordanMPOTensor)
-        E = $f(scalartype(W))
-        W′ = similar(W, E)
-        for (I, v) in nonzero_pairs(W.A)
-            W′.A[I] = $f(v)
-        end
-        for (I, v) in nonzero_pairs(W.B)
-            W′.B[I] = $f(v)
+        W′ = similar(W, $f(scalartype(W)))
+        for (I, v) in nonzero_pairs(getfield(W, :tensors))
+            getfield(W′, :tensors)[I] = $f(v)
         end
-        for (I, v) in nonzero_pairs(W.C)
-            W′.C[I] = $f(v)
-        end
-        for (I, v) in nonzero_pairs(W.D)
-            W′.D[I] = $f(v)
+        empty!(getfield(W′, :scalars))
+        for (K, c) in getfield(W, :scalars)
+            getfield(W′, :scalars)[K] = $f(c)
         end
         return W′
     end
@@ -213,19 +240,13 @@ end
 @inline Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4}) = W[I.I...]
 @propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4})
     @assert I[2] == I[3] == 1
-    i = I[1]
-    j = I[4]
-    if (size(W, 4) > 1 && i == 1 && j == 1) ||
-            (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4))
-        return BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[1])
-    elseif i == 1 && j == size(W, 4)
-        return insertrightunit(insertleftunit(only(W.D), 1), 3)
-    elseif i == 1
-        return insertleftunit(W.C[1, 1, j - 1], 1)
-    elseif j == size(W, 4)
-        return insertrightunit(W.B[i - 1, 1, 1], 3)
-    elseif 1 < i < size(W, 1) && 1 < j < size(W, 4)
-        return W.A[i - 1, 1, 1, j - 1]
+    i, j = I[1], I[4]
+    c = get(getfield(W, :scalars), CartesianIndex(i, j), nothing)
+    if c !== nothing
+        τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[i, 1, 1, j])
+        return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c)
+    elseif haskey(getfield(W, :tensors), CartesianIndex(i, 1, 1, j))
+        return getfield(W, :tensors)[i, 1, 1, j]
     else
         return zeros(storagetype(W), eachspace(W)[i, 1, 1, j])
     end
@@ -238,21 +259,14 @@ end
         W::JordanMPOTensor, v::MPOTensor, I::Vararg{Int, 4}
     )
     @assert I[2] == I[3] == 1
-    i = I[1]
-    j = I[4]
-    if i == 1 && j == size(W, 4)
-        W.D[1] = removeunit(removeunit(v, 4), 1)
-    elseif i == 1 && 1 < j < size(W, 4)
-        W.C[1, 1, j - 1] = removeunit(v, 1)
-    elseif j == size(W, 4) && 1 < i < size(W, 1)
-        W.B[i - 1, 1, 1] = removeunit(v, 4)
-    elseif 1 < i < size(W, 1) && 1 < j < size(W, 4)
-        W.A[i - 1, 1, 1, j - 1] = v
-    elseif (size(W, 4) > 1 && i == 1 && j == 1) ||
-            (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4))
-        v isa BraidingTensor || throw(ArgumentError("Cannot set BraidingTensor"))
+    i, j = I[1], I[4]
+    tensors, scalars = getfield(W, :tensors), getfield(W, :scalars)
+    if v isa BraidingTensor
+        haskey(tensors, CartesianIndex(i, 1, 1, j)) && delete!(tensors, CartesianIndex(i, 1, 1, j))
+        scalars[CartesianIndex(i, j)] = one(scalartype(W))
     else
-        throw(ArgumentError("Cannot set index ($i, 1, 1, $j)"))
+        delete!(scalars, CartesianIndex(i, j))
+        tensors[i, 1, 1, j] = v
     end
     return W
 end
@@ -263,31 +277,10 @@ end
 # Sparse functionality
 # --------------------
 function BlockTensorKit.nonzero_keys(W::JordanMPOTensor)
-    nrows = size(W, 1)
-    ncols = size(W, 4)
-    p = CartesianIndex{4}[]
-    ncols > 1 && push!(p, CartesianIndex(1, 1, 1, 1))
-    nrows > 1 && push!(p, CartesianIndex(nrows, 1, 1, ncols))
-
-    Ia = CartesianIndex(1, 0, 0, 1)
-    for I in nonzero_keys(W.A)
-        push!(p, I + Ia)
-    end
-
-    Ib = CartesianIndex(1, 0, 0)
-    for I in nonzero_keys(W.B)
-        push!(p, CartesianIndex((I + Ib).I..., ncols))
+    p = collect(CartesianIndex{4}, nonzero_keys(getfield(W, :tensors)))
+    for K in keys(getfield(W, :scalars))
+        push!(p, CartesianIndex(K[1], 1, 1, K[2]))
     end
-
-    Ic = CartesianIndex(0, 0, 1)
-    for I in nonzero_keys(W.C)
-        push!(p, CartesianIndex(1, (I + Ic).I...))
-    end
-
-    for I in nonzero_keys(W.D)
-        push!(p, CartesianIndex(1, 1, 1, ncols))
-    end
-
     return p
 end
 function BlockTensorKit.nonzero_values(W::JordanMPOTensor)
@@ -297,8 +290,7 @@ function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor)
     return Iterators.map(I -> I => W[I], nonzero_keys(W))
 end
 function BlockTensorKit.nonzero_length(W::JordanMPOTensor)
-    return nonzero_length(W.A) + nonzero_length(W.B) + nonzero_length(W.C) +
-        nonzero_length(W.D) + Int(size(W, 1) > 1) + Int(size(W, 4) > 1)
+    return nonzero_length(getfield(W, :tensors)) + length(getfield(W, :scalars))
 end
 
 # linalg
@@ -343,32 +335,32 @@ function add_physical_charge(O::JordanMPOTensor, charge::Sector)
         fuse(physicalspace(O), auxspace) ←
         fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)
     Odst = JordanMPOTensor{scalartype(O)}(undef, Vdst)
-    for (I, v) in nonzero_pairs(O)
-        Odst[I] = add_physical_charge(v, charge)
+    for (I, v) in nonzero_pairs(getfield(O, :tensors))
+        getfield(Odst, :tensors)[I] = add_physical_charge(v, charge)
     end
+    merge!(getfield(Odst, :scalars), getfield(O, :scalars))
     return Odst
 end
 
 # Utility
 # -------
 function Base.copy(W::JordanMPOTensor)
-    return JordanMPOTensor(W.V, copy(W.A), copy(W.B), copy(W.C), copy(W.D))
+    return typeof(W)(copy(getfield(W, :tensors)), copy(getfield(W, :scalars)))
 end
 function Base.copy!(Wdst::JordanMPOTensor, Wsrc::JordanMPOTensor)
     space(Wdst) == space(Wsrc) || throw(SpaceMismatch())
-    copy!(Wdst.A, Wsrc.A)
-    copy!(Wdst.B, Wsrc.B)
-    copy!(Wdst.C, Wsrc.C)
-    copy!(Wdst.D, Wsrc.D)
+    copy!(getfield(Wdst, :tensors), getfield(Wsrc, :tensors))
+    empty!(getfield(Wdst, :scalars))
+    merge!(getfield(Wdst, :scalars), getfield(Wsrc, :scalars))
     return Wdst
 end
 
 # Avoid falling back to `norm(W1 - W2)` which has to convert to SparseBlockTensorMap
 function Base.isapprox(W1::JordanMPOTensor, W2::JordanMPOTensor; kwargs...)
-    return isapprox(W1.A, W2.A; kwargs...) &&
-        isapprox(W1.B, W2.B; kwargs...) &&
-        isapprox(W1.C, W2.C; kwargs...) &&
-        isapprox(W1.D, W2.D; kwargs...)
+    isapprox(getfield(W1, :tensors), getfield(W2, :tensors); kwargs...) || return false
+    s1, s2 = getfield(W1, :scalars), getfield(W2, :scalars)
+    keys(s1) == keys(s2) || return false
+    return all(k -> isapprox(s1[k], s2[k]; kwargs...), keys(s1))
 end
 
 function Base.showarg(io::IO, W::JordanMPOTensor, toplevel::Bool)

From b76c582ab5b74388b8e153e7989e309bf2fd07fa Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 24 Jun 2026 21:22:43 -0400
Subject: [PATCH 3/9] Propagate JordanMPOTensor rework to consumers

- MPOHamiltonian getproperty `:A` delegates to the tensor-level accessor so
  identity levels stay `BraidingTensor`s (keeps `isidentitylevel` correct).
- Rewrite `scale!` to act directly on `tensors`/`scalars` (the `.A/.B/.C/.D`
  accessors share data via `removeunit`, so in-place scaling of a copy leaked
  into the source).
- Add a 2-field `JordanMPOTensor(tensors, scalars)` constructor and use it in
  the Adapt extension and the `force_planar` test helper.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 ext/MPSKitAdaptExt.jl            |  2 +-
 src/operators/jordanmpotensor.jl |  7 +++++++
 src/operators/mpohamiltonian.jl  | 29 +++++++++++++++++++----------
 test/setup/testsetup.jl          | 11 +----------
 4 files changed, 28 insertions(+), 21 deletions(-)

diff --git a/ext/MPSKitAdaptExt.jl b/ext/MPSKitAdaptExt.jl
index 78f4b41af..3648b0fc4 100644
--- a/ext/MPSKitAdaptExt.jl
+++ b/ext/MPSKitAdaptExt.jl
@@ -32,7 +32,7 @@ end
 # inline to improve type stability with closures
 @inline Adapt.adapt_structure(to, mpo::MPO) = MPO(map(adapt(to), mpo.O))
 @inline Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor) =
-    MPSKit.JordanMPOTensor(space(W), adapt(to, W.A), adapt(to, W.B), adapt(to, W.C), adapt(to, W.D))
+    MPSKit.JordanMPOTensor(adapt(to, W.tensors), copy(W.scalars))
 @inline Adapt.adapt_structure(to, mpo::MPOHamiltonian) =
     MPOHamiltonian(map(x -> adapt(to, x), mpo.W))
 @inline Adapt.adapt_structure(to, ml::MPSKit.Multiline) = MPSKit.Multiline(map(adapt(to), ml.data))
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index b47efa2af..3b9ea8e51 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -42,6 +42,13 @@ struct JordanMPOTensor{
     end
 end
 
+function JordanMPOTensor(
+        tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4},
+        scalars::Dict{CartesianIndex{2}, T}
+    ) where {T, S, A}
+    return JordanMPOTensor{T, S, A}(tensors, scalars)
+end
+
 function JordanMPOTensor{E, S}(::UndefInitializer, V::TensorMapSumSpace{S}) where {E, S}
     return jordanmpotensortype(S, E)(undef, V)
 end
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 401064493..665400069 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -633,7 +633,7 @@ Base.copy(H::MPOHamiltonian) = MPOHamiltonian(map(copy, parent(H)))
 
 function Base.getproperty(H::MPOHamiltonian, sym::Symbol)
     if sym === :A
-        return map(h -> h[2:(end - 1), 1, 1, 2:(end - 1)], parent(H))
+        return map(h -> h.A, parent(H))
     elseif sym === :B
         return map(h -> h[2:(end - 1), 1, 1, end], parent(H))
     elseif sym === :C
@@ -782,12 +782,18 @@ Base.:-(H::MPOHamiltonian, λs::AbstractVector{<:Number}) = H + (-λs)
 Base.:-(λs::AbstractVector{<:Number}, H::MPOHamiltonian) = λs + (-H)
 Base.:-(H1::MPOHamiltonian, H2::MPOHamiltonian) = H1 + (-H2)
 
+# scaling a Jordan MPO Hamiltonian scales every path exactly once, by scaling the
+# transitions out of the starting level (the top row, excluding the identity corner)
 function VectorInterface.scale!(
         H::MPOHamiltonian{O}, λ::Number
     ) where {O <: JordanMPOTensor}
-    for i in 1:length(H)
-        scale!(H[i].C, λ)
-        scale!(H[i].D, λ)
+    for W in parent(H)
+        for (I, v) in nonzero_pairs(W.tensors)
+            I[1] == 1 && scale!(v, λ)
+        end
+        for K in collect(keys(W.scalars))
+            (K[1] == 1 && K[2] != 1) && (W.scalars[K] *= λ)
+        end
     end
     return H
 end
@@ -795,12 +801,15 @@ function VectorInterface.scale!(
         Hdst::MPOHamiltonian{<:JordanMPOTensor},
         Hsrc::MPOHamiltonian{<:JordanMPOTensor}, λ::Number
     )
-    N = check_length(Hdst, Hsrc)
-    for i in 1:N
-        scale!(Hdst[i].C, Hsrc[i].C, λ)
-        scale!(Hdst[i].D, Hsrc[i].D, λ)
-        copy!(Hdst[i].A, Hsrc[i].A)
-        copy!(Hdst[i].B, Hsrc[i].B)
+    check_length(Hdst, Hsrc)
+    for (Wd, Ws) in zip(parent(Hdst), parent(Hsrc))
+        for (I, v) in nonzero_pairs(Ws.tensors)
+            Wd.tensors[I] = I[1] == 1 ? scale(v, λ) : copy(v)
+        end
+        empty!(Wd.scalars)
+        for (K, c) in Ws.scalars
+            Wd.scalars[K] = (K[1] == 1 && K[2] != 1) ? c * λ : c
+        end
     end
     return Hdst
 end
diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl
index 68c99f096..053d272f3 100644
--- a/test/setup/testsetup.jl
+++ b/test/setup/testsetup.jl
@@ -63,16 +63,7 @@ function force_planar(x::SparseBlockTensorMap)
     return SparseBlockTensorMap{valtype(data)}(data, force_planar(space(x)))
 end
 function force_planar(W::JordanMPOTensor)
-    V = force_planar(space(W))
-    TW = MPSKit.jordanmpotensortype(eltype(V[1]), scalartype(W))
-    dst = TW(undef, V)
-
-    for t in (:A, :B, :C, :D)
-        for (I, v) in nonzero_pairs(getproperty(W, t))
-            getproperty(dst, t)[I] = force_planar(v)
-        end
-    end
-    return dst
+    return JordanMPOTensor(force_planar(W.tensors), copy(W.scalars))
 end
 force_planar(mpo::MPOHamiltonian) = MPOHamiltonian(map(force_planar, parent(mpo)))
 force_planar(mpo::MPO) = MPO(map(force_planar, parent(mpo)))

From e58323b2a161b242b04246f99f470c37bc985472 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 24 Jun 2026 22:05:09 -0400
Subject: [PATCH 4/9] Fix dual boundary units and adapt scalartype; use
 property access

- 5-arg constructor: insert the trivial boundary legs with the duality dictated
  by V (via isdual of the corner space), so conjugated MPOs (_conj_mpo) embed
  correctly instead of hitting a SpaceMismatch.
- Adapt extension: convert the `scalars` coefficients to the adapted scalartype
  so precision changes (e.g. ComplexF64 -> ComplexF32) construct cleanly.
- Replace internal `getfield(W, :tensors/:scalars)` with `W.tensors`/`W.scalars`
  (the getproperty passthrough), keeping only the one passthrough getfield.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 ext/MPSKitAdaptExt.jl            | 10 ++--
 src/operators/jordanmpotensor.jl | 86 ++++++++++++++++----------------
 2 files changed, 51 insertions(+), 45 deletions(-)

diff --git a/ext/MPSKitAdaptExt.jl b/ext/MPSKitAdaptExt.jl
index 3648b0fc4..9d550c75d 100644
--- a/ext/MPSKitAdaptExt.jl
+++ b/ext/MPSKitAdaptExt.jl
@@ -1,6 +1,6 @@
 module MPSKitAdaptExt
 
-using TensorKit: space, spacetype
+using TensorKit: space, spacetype, scalartype
 using MPSKit
 using BlockTensorKit: nonzero_pairs
 using Adapt
@@ -31,8 +31,12 @@ end
 
 # inline to improve type stability with closures
 @inline Adapt.adapt_structure(to, mpo::MPO) = MPO(map(adapt(to), mpo.O))
-@inline Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor) =
-    MPSKit.JordanMPOTensor(adapt(to, W.tensors), copy(W.scalars))
+@inline function Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor)
+    tensors = adapt(to, W.tensors)
+    T = scalartype(tensors)
+    scalars = Dict{CartesianIndex{2}, T}(K => convert(T, c) for (K, c) in W.scalars)
+    return MPSKit.JordanMPOTensor(tensors, scalars)
+end
 @inline Adapt.adapt_structure(to, mpo::MPOHamiltonian) =
     MPOHamiltonian(map(x -> adapt(to, x), mpo.W))
 @inline Adapt.adapt_structure(to, ml::MPSKit.Multiline) = MPSKit.Multiline(map(adapt(to), ml.data))
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index 3b9ea8e51..e7095b10b 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -76,17 +76,23 @@ function JordanMPOTensor(
 
     W = jordanmpotensortype(S, storagetype(A))(undef, V)
     cols = size(W, 4)
+    # the trivial boundary legs inherit their duality from the virtual spaces of V
+    corner = allVs[1, 1, 1, end]
+    dual_left = isdual(corner[1])
+    dual_right = !isdual(corner[4])
     for (I, v) in nonzero_pairs(A)
         W[I[1] + 1, 1, 1, I[4] + 1] = v
     end
     for (I, v) in nonzero_pairs(B)
-        W[I[1] + 1, 1, 1, cols] = insertrightunit(v, 3)
+        W[I[1] + 1, 1, 1, cols] = insertrightunit(v, 3; dual = dual_right)
     end
     for (I, v) in nonzero_pairs(C)
-        W[1, 1, 1, I[3] + 1] = insertleftunit(v, 1)
+        W[1, 1, 1, I[3] + 1] = insertleftunit(v, 1; dual = dual_left)
     end
     if nonzero_length(D) > 0
-        W[1, 1, 1, cols] = insertrightunit(insertleftunit(only(D), 1), 3)
+        W[1, 1, 1, cols] = insertrightunit(
+            insertleftunit(only(D), 1; dual = dual_left), 3; dual = dual_right
+        )
     end
     return W
 end
@@ -132,7 +138,7 @@ end
 
 # Properties
 # ----------
-TensorKit.space(W::JordanMPOTensor) = space(getfield(W, :tensors))
+TensorKit.space(W::JordanMPOTensor) = space(W.tensors)
 function Base.eltype(::Type{<:JordanMPOTensor{T, S, A}}) where {T, S, A}
     return tensormaptype(S, 2, 2, A)
 end
@@ -150,14 +156,13 @@ function _jordan_A(W::JordanMPOTensor)
     rows, cols = size(W, 1), size(W, 4)
     Tτ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}
     TA = tensormaptype(spacetype(W), 2, 2, storagetype(W))
-    tensors = getfield(W, :tensors)
-    VA = space(eachspace(tensors)[2:(end - 1), 1, 1, 2:(end - 1)])
+    VA = space(eachspace(W.tensors)[2:(end - 1), 1, 1, 2:(end - 1)])
     A = SparseBlockTensorMap{Union{TA, Tτ}}(undef, VA)
-    for (I, v) in nonzero_pairs(tensors)
+    for (I, v) in nonzero_pairs(W.tensors)
         (1 < I[1] < rows && 1 < I[4] < cols) || continue
         A[I[1] - 1, 1, 1, I[4] - 1] = v
     end
-    for (K, c) in getfield(W, :scalars)
+    for (K, c) in W.scalars
         i, j = K[1], K[2]
         (1 < i < rows && 1 < j < cols) || continue
         τ = Tτ(eachspace(A)[i - 1, 1, 1, j - 1])
@@ -168,7 +173,7 @@ end
 function _jordan_B(W::JordanMPOTensor)
     rows, cols = size(W, 1), size(W, 4)
     TB = tensormaptype(spacetype(W), 2, 1, storagetype(W))
-    VB = removeunit(space(eachspace(getfield(W, :tensors))[2:(end - 1), 1, 1, end]), 4)
+    VB = removeunit(space(eachspace(W.tensors)[2:(end - 1), 1, 1, end]), 4)
     B = SparseBlockTensorMap{TB}(undef, VB)
     for I in nonzero_keys(W)
         (1 < I[1] < rows && I[4] == cols) || continue
@@ -179,7 +184,7 @@ end
 function _jordan_C(W::JordanMPOTensor)
     cols = size(W, 4)
     TC = tensormaptype(spacetype(W), 1, 2, storagetype(W))
-    VC = removeunit(space(eachspace(getfield(W, :tensors))[1, 1, 1, 2:(end - 1)]), 1)
+    VC = removeunit(space(eachspace(W.tensors)[1, 1, 1, 2:(end - 1)]), 1)
     C = SparseBlockTensorMap{TC}(undef, VC)
     for I in nonzero_keys(W)
         (I[1] == 1 && 1 < I[4] < cols) || continue
@@ -190,7 +195,7 @@ end
 function _jordan_D(W::JordanMPOTensor)
     cols = size(W, 4)
     TD = tensormaptype(spacetype(W), 1, 1, storagetype(W))
-    VD = removeunit(removeunit(space(eachspace(getfield(W, :tensors))[1, 1, 1, end:end]), 4), 1)
+    VD = removeunit(removeunit(space(eachspace(W.tensors)[1, 1, 1, end:end]), 4), 1)
     D = SparseBlockTensorMap{TD}(undef, VD)
     K = CartesianIndex(1, 1, 1, cols)
     if haskey(W, K)
@@ -201,8 +206,7 @@ end
 
 function Base.haskey(W::JordanMPOTensor, I::CartesianIndex{4})
     Base.checkbounds(W, I.I...)
-    return haskey(getfield(W, :scalars), CartesianIndex(I[1], I[4])) ||
-        haskey(getfield(W, :tensors), I)
+    return haskey(W.scalars, CartesianIndex(I[1], I[4])) || haskey(W.tensors, I)
 end
 
 Base.parent(W::JordanMPOTensor) = parent(SparseBlockTensorMap(W))
@@ -215,10 +219,10 @@ function BlockTensorKit.SparseBlockTensorMap(W::JordanMPOTensor)
     W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W), spacetype(W), 2, 2}}(
         undef_blocks, space(W)
     )
-    for (I, v) in nonzero_pairs(getfield(W, :tensors))
+    for (I, v) in nonzero_pairs(W.tensors)
         W′[I] = v
     end
-    for (K, c) in getfield(W, :scalars)
+    for (K, c) in W.scalars
         τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(
             eachspace(W)[K[1], 1, 1, K[2]]
         )
@@ -230,12 +234,12 @@ end
 for f in (:real, :complex)
     @eval function Base.$f(W::JordanMPOTensor)
         W′ = similar(W, $f(scalartype(W)))
-        for (I, v) in nonzero_pairs(getfield(W, :tensors))
-            getfield(W′, :tensors)[I] = $f(v)
+        for (I, v) in nonzero_pairs(W.tensors)
+            W′.tensors[I] = $f(v)
         end
-        empty!(getfield(W′, :scalars))
-        for (K, c) in getfield(W, :scalars)
-            getfield(W′, :scalars)[K] = $f(c)
+        empty!(W′.scalars)
+        for (K, c) in W.scalars
+            W′.scalars[K] = $f(c)
         end
         return W′
     end
@@ -248,12 +252,12 @@ end
 @propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4})
     @assert I[2] == I[3] == 1
     i, j = I[1], I[4]
-    c = get(getfield(W, :scalars), CartesianIndex(i, j), nothing)
+    c = get(W.scalars, CartesianIndex(i, j), nothing)
     if c !== nothing
         τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[i, 1, 1, j])
         return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c)
-    elseif haskey(getfield(W, :tensors), CartesianIndex(i, 1, 1, j))
-        return getfield(W, :tensors)[i, 1, 1, j]
+    elseif haskey(W.tensors, CartesianIndex(i, 1, 1, j))
+        return W.tensors[i, 1, 1, j]
     else
         return zeros(storagetype(W), eachspace(W)[i, 1, 1, j])
     end
@@ -267,13 +271,12 @@ end
     )
     @assert I[2] == I[3] == 1
     i, j = I[1], I[4]
-    tensors, scalars = getfield(W, :tensors), getfield(W, :scalars)
     if v isa BraidingTensor
-        haskey(tensors, CartesianIndex(i, 1, 1, j)) && delete!(tensors, CartesianIndex(i, 1, 1, j))
-        scalars[CartesianIndex(i, j)] = one(scalartype(W))
+        haskey(W.tensors, CartesianIndex(i, 1, 1, j)) && delete!(W.tensors, CartesianIndex(i, 1, 1, j))
+        W.scalars[CartesianIndex(i, j)] = one(scalartype(W))
     else
-        delete!(scalars, CartesianIndex(i, j))
-        tensors[i, 1, 1, j] = v
+        delete!(W.scalars, CartesianIndex(i, j))
+        W.tensors[i, 1, 1, j] = v
     end
     return W
 end
@@ -284,8 +287,8 @@ end
 # Sparse functionality
 # --------------------
 function BlockTensorKit.nonzero_keys(W::JordanMPOTensor)
-    p = collect(CartesianIndex{4}, nonzero_keys(getfield(W, :tensors)))
-    for K in keys(getfield(W, :scalars))
+    p = collect(CartesianIndex{4}, nonzero_keys(W.tensors))
+    for K in keys(W.scalars)
         push!(p, CartesianIndex(K[1], 1, 1, K[2]))
     end
     return p
@@ -297,7 +300,7 @@ function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor)
     return Iterators.map(I -> I => W[I], nonzero_keys(W))
 end
 function BlockTensorKit.nonzero_length(W::JordanMPOTensor)
-    return nonzero_length(getfield(W, :tensors)) + length(getfield(W, :scalars))
+    return nonzero_length(W.tensors) + length(W.scalars)
 end
 
 # linalg
@@ -342,32 +345,31 @@ function add_physical_charge(O::JordanMPOTensor, charge::Sector)
         fuse(physicalspace(O), auxspace) ←
         fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)
     Odst = JordanMPOTensor{scalartype(O)}(undef, Vdst)
-    for (I, v) in nonzero_pairs(getfield(O, :tensors))
-        getfield(Odst, :tensors)[I] = add_physical_charge(v, charge)
+    for (I, v) in nonzero_pairs(O.tensors)
+        Odst.tensors[I] = add_physical_charge(v, charge)
     end
-    merge!(getfield(Odst, :scalars), getfield(O, :scalars))
+    merge!(Odst.scalars, O.scalars)
     return Odst
 end
 
 # Utility
 # -------
 function Base.copy(W::JordanMPOTensor)
-    return typeof(W)(copy(getfield(W, :tensors)), copy(getfield(W, :scalars)))
+    return typeof(W)(copy(W.tensors), copy(W.scalars))
 end
 function Base.copy!(Wdst::JordanMPOTensor, Wsrc::JordanMPOTensor)
     space(Wdst) == space(Wsrc) || throw(SpaceMismatch())
-    copy!(getfield(Wdst, :tensors), getfield(Wsrc, :tensors))
-    empty!(getfield(Wdst, :scalars))
-    merge!(getfield(Wdst, :scalars), getfield(Wsrc, :scalars))
+    copy!(Wdst.tensors, Wsrc.tensors)
+    empty!(Wdst.scalars)
+    merge!(Wdst.scalars, Wsrc.scalars)
     return Wdst
 end
 
 # Avoid falling back to `norm(W1 - W2)` which has to convert to SparseBlockTensorMap
 function Base.isapprox(W1::JordanMPOTensor, W2::JordanMPOTensor; kwargs...)
-    isapprox(getfield(W1, :tensors), getfield(W2, :tensors); kwargs...) || return false
-    s1, s2 = getfield(W1, :scalars), getfield(W2, :scalars)
-    keys(s1) == keys(s2) || return false
-    return all(k -> isapprox(s1[k], s2[k]; kwargs...), keys(s1))
+    isapprox(W1.tensors, W2.tensors; kwargs...) || return false
+    keys(W1.scalars) == keys(W2.scalars) || return false
+    return all(k -> isapprox(W1.scalars[k], W2.scalars[k]; kwargs...), keys(W1.scalars))
 end
 
 function Base.showarg(io::IO, W::JordanMPOTensor, toplevel::Bool)

From 6ed544e2b7d0f9632ee67d868fea7180f089ee22 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 25 Jun 2026 10:48:17 -0400
Subject: [PATCH 5/9] Optimize JordanMPOTensor block access in consumers

The A/B/C/D blocks are computed properties that rebuild a fresh
SparseBlockTensorMap on every access. Avoid redundant reconstruction:

- jordanmpotensor: rewrite `_jordan_B`/`_jordan_C` and `nonzero_values`/
  `nonzero_pairs` to iterate `tensors`/`scalars` directly instead of
  `nonzero_keys(W)` + `W[I]`; add `_identity_tensor` helper.
- hamiltonian_derivatives: bind `W.A/.B/.C/.D` once in the AC and AC2
  derivative builders (W*.A was rebuilt ~8x per AC2 build).
- ortho: cache blocks in left/right_canonicalize! and drop
  `typeof(W.A)`-style reconstructions.
- mpohamiltonian: answer `isidentitylevel` from `scalars` instead of
  reconstructing every site's A block.

Also expand the JordanMPOTensor docstring.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 .../derivatives/hamiltonian_derivatives.jl    | 61 +++++++------
 src/operators/jordanmpotensor.jl              | 78 +++++++++++++---
 src/operators/mpohamiltonian.jl               |  8 +-
 src/operators/ortho.jl                        | 90 ++++++++++---------
 4 files changed, 154 insertions(+), 83 deletions(-)

diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl
index fb93d67b6..47ee2b25a 100644
--- a/src/algorithms/derivatives/hamiltonian_derivatives.jl
+++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl
@@ -48,8 +48,12 @@ function AC_hamiltonian(
 end
 
 function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor)
+    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
+    # them once and reuse the locals throughout
+    WA, WB, WC, WD = W.A, W.B, W.C, W.D
+
     # onsite
-    D = nonzero_length(W.D) > 0 ? only(W.D) : missing
+    D = nonzero_length(WD) > 0 ? only(WD) : missing
 
     # not started
     I = size(W, 4) == 1 ? missing : removeunit(GR[1], 2)
@@ -58,25 +62,25 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens
     E = size(W, 1) == 1 ? missing : removeunit(GL[end], 2)
 
     # starting
-    C = if nonzero_length(W.C) > 0
+    C = if nonzero_length(WC) > 0
         GR_2 = GR[2:(end - 1)]
-        @plansor starting[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2]
+        @plansor starting[-1 -2; -3 -4] ≔ WC[-1; -3 1] * GR_2[-4 1; -2]
         only(starting)
     else
         missing
     end
 
     # ending
-    B = if nonzero_length(W.B) > 0
+    B = if nonzero_length(WB) > 0
         GL_2 = GL[2:(end - 1)]
-        @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4]
+        @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * WB[1 -2; -4]
         only(ending)
     else
         missing
     end
 
     # continuing
-    A = MPO_AC_Hamiltonian(GL[2:(end - 1)], W.A, GR[2:(end - 1)])
+    A = MPO_AC_Hamiltonian(GL[2:(end - 1)], WA, GR[2:(end - 1)])
 
     # obtaining storagetype of environments since these should have already mixed
     # the types of the operator and state
@@ -87,7 +91,7 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens
     O3 = typeof(A)
 
     # specialization for nearest neighbours
-    nonzero_length(W.A) == 0 && (A = missing)
+    nonzero_length(WA) == 0 && (A = missing)
 
     return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A)
 end
@@ -196,6 +200,11 @@ function AC2_hamiltonian(
 end
 
 function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor)
+    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
+    # them once and reuse the locals throughout
+    A1, B1, C1, D1 = W1.A, W1.B, W1.C, W1.D
+    A2, B2, C2, D2 = W2.A, W2.B, W2.C, W2.D
+
     # not started
     II = size(W2, 4) == 1 ? missing : transpose(removeunit(GR[1], 2))
 
@@ -203,22 +212,22 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda
     EE = size(W1, 1) == 1 ? missing : removeunit(GL[end], 2)
 
     # starting right
-    IC = if nonzero_length(W2.C) > 0
-        @plansor IC_[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2]
+    IC = if nonzero_length(C2) > 0
+        @plansor IC_[-1 -2; -3 -4] ≔ C2[-1; -3 1] * GR[2:(end - 1)][-4 1; -2]
         only(IC_)
     else
         missing
     end
 
     # onsite left
-    DE = nonzero_length(W1.D) > 0 ? only(W1.D) : missing
+    DE = nonzero_length(D1) > 0 ? only(D1) : missing
 
     # onsite right
-    ID = nonzero_length(W2.D) > 0 ? only(W2.D) : missing
+    ID = nonzero_length(D2) > 0 ? only(D2) : missing
 
     # starting left - ending right
-    CB = if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0
-        @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4]
+    CB = if nonzero_length(C1) > 0 && nonzero_length(B2) > 0
+        @plansor CB_[-1 -2; -3 -4] ≔ C1[-1; -3 1] * B2[1 -2; -4]
         # have to convert to complex if hamiltonian is real but states are complex
         scalartype(GL) <: Complex ? complex(only(CB_)) : only(CB_)
     else
@@ -226,8 +235,8 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda
     end
 
     # starting left - continuing right
-    CA = if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0
-        @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] *
+    CA = if nonzero_length(C1) > 0 && nonzero_length(A2) > 0
+        @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ C1[-1; -4 2] * A2[2 -2; -5 1] *
             GR[2:(end - 1)][-6 1; -3]
         only(CA_)
     else
@@ -235,24 +244,24 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda
     end
 
     # continuing left - ending right
-    AB = if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0
-        @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * W1.A[2 -2; -5 1] *
-            W2.B[1 -3; -6]
+    AB = if nonzero_length(A1) > 0 && nonzero_length(B2) > 0
+        @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * A1[2 -2; -5 1] *
+            B2[1 -3; -6]
         only(AB_)
     else
         missing
     end
 
     # ending left
-    BE = if nonzero_length(W1.B) > 0
-        @plansor BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4]
+    BE = if nonzero_length(B1) > 0
+        @plansor BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * B1[2 -2; -4]
         only(BE_)
     else
         missing
     end
 
     # continuing - continuing
-    AA = MPO_AC2_Hamiltonian(GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)])
+    AA = MPO_AC2_Hamiltonian(GL[2:(end - 1)], A1, A2, GR[2:(end - 1)])
 
     S = spacetype(GL)
     M = storagetype(GL)
@@ -261,16 +270,16 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda
     O3 = tensormaptype(S, 3, 3, M)
     O4 = typeof(AA)
 
-    if nonzero_length(W1.A) == 0 && nonzero_length(W2.A) == 0
+    if nonzero_length(A1) == 0 && nonzero_length(A2) == 0
         AA = missing
     else
-        mask1 = falses(size(W1.A, 1), size(W1.A, 4))
-        for I in nonzero_keys(W1.A)
+        mask1 = falses(size(A1, 1), size(A1, 4))
+        for I in nonzero_keys(A1)
             mask1[I[1], I[4]] = true
         end
 
-        mask2 = falses(size(W2.A, 1), size(W2.A, 4))
-        for I in nonzero_keys(W2.A)
+        mask2 = falses(size(A2, 1), size(A2, 4))
+        for I in nonzero_keys(A2)
             mask2[I[1], I[4]] = true
         end
 
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index e7095b10b..253055a26 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -1,19 +1,43 @@
 """
     JordanMPOTensor{T,S,A} <: AbstractBlockTensorMap{T,S,2,2}
 
-A tensor map that represents the upper triangular block form of a matrix product operator (MPO).
+A single tensor of a matrix product operator (MPO) in upper triangular (Jordan) block form,
+as used to represent the local tensors of an [`MPOHamiltonian`](@ref).
+The virtual (row, column) structure is
 
 ```math
 \\begin{pmatrix}
 1 & C & D \\\\
-0 & A & B \\\\
-0 & 0 & 1
+⋅ & A & B \\\\
+⋅ & ⋅ & 1
 \\end{pmatrix}
 ```
 
-The genuine operators are stored in `tensors` over the full virtual space, while the
-scalar multiples of the identity (in particular the diagonal `1`s) are kept separately in
-`scalars`, keyed by their `(row, col)` virtual position.
+where `A` is the bulk of interacting operators, `C`/`B` are the operators that start/finish
+an interaction, `D` is the on-site term, and the diagonal `1`s are identities.
+
+## Representation
+
+Rather than storing the dense block matrix, the genuine operators and the identities are kept separately:
+
+- `tensors::SparseBlockTensorMap` holds the non-identity operators over the *full* virtual
+  space (so `A`, `B`, `C` and `D` all live at their `(row, 1, 1, col)` position).
+- `scalars::Dict{CartesianIndex{2},T}` holds the scalar multiples of the identity, keyed by
+  their `(row, col)` virtual position; the diagonal corner `1`s are stored here as well.
+
+An index is never present in both `tensors` and `scalars`.
+This keeps the ubiquitous identity blocks free of dense storage and lets identities be materialized lazily only when needed.
+
+## Type parameters
+
+- `T <: Number`: the `scalartype` of the tensors.
+- `S`: the `spacetype` of the tensors.
+- `A <: DenseVector{T}`: the storage type of the underlying tensors.
+
+## Block accessors
+
+The reduced-leg `A`, `B`, `C` and `D` blocks are exposed as properties (`W.A`, `W.B`, `W.C`, `W.D`),
+reconstructed on demand from `tensors` and `scalars`.
 """
 struct JordanMPOTensor{
         T <: Number, S, A <: DenseVector{T},
@@ -151,6 +175,12 @@ function Base.getproperty(W::JordanMPOTensor, sym::Symbol)
     return getfield(W, sym)
 end
 
+# materialize an identity scalar at virtual position `(i, j)` as a dense `TensorMap`
+@inline function _identity_tensor(W::JordanMPOTensor, i::Int, j::Int, c = one(scalartype(W)))
+    τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[i, 1, 1, j])
+    return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c)
+end
+
 # reduced-leg blocks reconstructed from `tensors` and `scalars`
 function _jordan_A(W::JordanMPOTensor)
     rows, cols = size(W, 1), size(W, 4)
@@ -175,9 +205,13 @@ function _jordan_B(W::JordanMPOTensor)
     TB = tensormaptype(spacetype(W), 2, 1, storagetype(W))
     VB = removeunit(space(eachspace(W.tensors)[2:(end - 1), 1, 1, end]), 4)
     B = SparseBlockTensorMap{TB}(undef, VB)
-    for I in nonzero_keys(W)
+    for (I, v) in nonzero_pairs(W.tensors)
         (1 < I[1] < rows && I[4] == cols) || continue
-        B[I[1] - 1, 1, 1] = removeunit(W[I], 4)
+        B[I[1] - 1, 1, 1] = removeunit(v, 4)
+    end
+    for (K, c) in W.scalars
+        (1 < K[1] < rows && K[2] == cols) || continue
+        B[K[1] - 1, 1, 1] = removeunit(_identity_tensor(W, K[1], K[2], c), 4)
     end
     return B
 end
@@ -186,9 +220,13 @@ function _jordan_C(W::JordanMPOTensor)
     TC = tensormaptype(spacetype(W), 1, 2, storagetype(W))
     VC = removeunit(space(eachspace(W.tensors)[1, 1, 1, 2:(end - 1)]), 1)
     C = SparseBlockTensorMap{TC}(undef, VC)
-    for I in nonzero_keys(W)
+    for (I, v) in nonzero_pairs(W.tensors)
         (I[1] == 1 && 1 < I[4] < cols) || continue
-        C[1, 1, I[4] - 1] = removeunit(W[I], 1)
+        C[1, 1, I[4] - 1] = removeunit(v, 1)
+    end
+    for (K, c) in W.scalars
+        (K[1] == 1 && 1 < K[2] < cols) || continue
+        C[1, 1, K[2] - 1] = removeunit(_identity_tensor(W, K[1], K[2], c), 1)
     end
     return C
 end
@@ -254,8 +292,7 @@ end
     i, j = I[1], I[4]
     c = get(W.scalars, CartesianIndex(i, j), nothing)
     if c !== nothing
-        τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[i, 1, 1, j])
-        return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c)
+        return _identity_tensor(W, i, j, c)
     elseif haskey(W.tensors, CartesianIndex(i, 1, 1, j))
         return W.tensors[i, 1, 1, j]
     else
@@ -294,10 +331,23 @@ function BlockTensorKit.nonzero_keys(W::JordanMPOTensor)
     return p
 end
 function BlockTensorKit.nonzero_values(W::JordanMPOTensor)
-    return Iterators.map(I -> W[I], nonzero_keys(W))
+    return Iterators.flatten(
+        (
+            nonzero_values(W.tensors),
+            (_identity_tensor(W, K[1], K[2], c) for (K, c) in W.scalars),
+        )
+    )
 end
 function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor)
-    return Iterators.map(I -> I => W[I], nonzero_keys(W))
+    return Iterators.flatten(
+        (
+            nonzero_pairs(W.tensors),
+            (
+                CartesianIndex(K[1], 1, 1, K[2]) => _identity_tensor(W, K[1], K[2], c)
+                    for (K, c) in W.scalars
+            ),
+        )
+    )
 end
 function BlockTensorKit.nonzero_length(W::JordanMPOTensor)
     return nonzero_length(W.tensors) + length(W.scalars)
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 665400069..067bd2a78 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -649,9 +649,11 @@ function isidentitylevel(H::InfiniteMPOHamiltonian{<:JordanMPOTensor}, i::Int)
     if i == 1 || i == size(H[1], 1)
         return true
     else
-        return all(H.A) do A
-            return haskey(A, CartesianIndex(i - 1, 1, 1, i - 1)) &&
-                A[i - 1, 1, 1, i - 1] isa BraidingTensor
+        # a diagonal level is an identity level iff every site stores a unit identity
+        # scalar there; pure identities live in `scalars` (genuine/scaled operators do not)
+        return all(parent(H)) do W
+            c = get(W.scalars, CartesianIndex(i, i), nothing)
+            return c !== nothing && isone(c)
         end
     end
 end
diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl
index a91442213..ac0a01ffb 100644
--- a/src/operators/ortho.jl
+++ b/src/operators/ortho.jl
@@ -10,14 +10,18 @@ function left_canonicalize!(
     S = spacetype(W)
     d = sqrt(dim(P))
 
+    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
+    # them once and reuse the locals throughout
+    WA, WB, WC, WD = W.A, W.B, W.C, W.D
+
     # orthogonalize second column against first
     Wi = W[1, 1, 1, 1]
     WI = removeunit(Wi, 1)
-    @plansor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r]
+    @plansor t[l; r] := conj(WI[p; p' l]) * WC[p; p' r]
     # TODO: the following is currently broken due to a TensorKit bug
-    # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r]
+    # @plansor C′[p; p' r] := WC[p; p' r] - WI[p; p' l] * t[l; r]
     @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r]
-    add!(C′, W.C)
+    add!(C′, WC)
 
     # QR of second column
     if size(W, 1) == 1
@@ -26,24 +30,24 @@ function left_canonicalize!(
 
         if dim(R) == 0 # fully truncated
             V = _rightunit ⊞ _rightunit
-            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
-            Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}())
+            Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
+            Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}())
         else
             V = ⊞(_rightunit, space(R, 1), _rightunit)
             scale!(Q, d)
             scale!(R, inv(d))
-            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1))
+            Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1))
             Q2 = transpose(Q, ((2,), (1, 3)))
         end
-        H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D)
+        H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, WB, Q2, WD)
     else
-        tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,)))
+        tmp = transpose(cat(insertleftunit(C′, 1), WA; dims = 1), ((3, 1, 2), (4,)))
         Q, R = left_orth!(tmp; alg)
 
         if dim(R) == 0 # fully truncated
             V = _rightunit ⊞ _rightunit
-            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
-            Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}())
+            Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
+            Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}())
         else
             scale!(Q, d)
             scale!(R, inv(d))
@@ -52,33 +56,34 @@ function left_canonicalize!(
             Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
             V = ⊞(_rightunit, right_virtualspace(Q′), _rightunit)
         end
-        H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, W.B, Q2, W.D)
+        H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD)
     end
 
     # absorb into next site
     W′ = H[i + 1]
+    W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D
 
     if size(W′, 4) > 1 && dim(R) != 0
-        @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r]
+        @plansor A′[l p; p' r] := R[l; r'] * W′A[r' p; p' r]
     else
-        A′ = similar(W′.A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.A))
+        A′ = similar(W′A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′A))
     end
 
     if size(W′, 4) > 1
-        @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r]
-        C′ = add!(removeunit(C′, 1), W′.C)
+        @plansor C′[l p; p' r] := t[l; r'] * W′A[r' p; p' r]
+        C′ = add!(removeunit(C′, 1), W′C)
     else
-        C′ = W′.C # empty
+        C′ = W′C # empty
     end
 
     if dim(R) != 0
-        @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p']
+        @plansor B′[l p; p'] := R[l; r] * W′B[r p; p']
     else
-        B′ = similar(W′.B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.B))
+        B′ = similar(W′B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′B))
     end
 
-    @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p']
-    D′ = add!(removeunit(D′, 1), W′.D)
+    @plansor D′[l p; p'] := t[l; r] * W′B[r p; p']
+    D′ = add!(removeunit(D′, 1), W′D)
 
     H[i + 1] = JordanMPOTensor(
         right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′),
@@ -99,14 +104,18 @@ function right_canonicalize!(
     S = spacetype(W)
     d = sqrt(dim(P))
 
+    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
+    # them once and reuse the locals throughout
+    WA, WB, WC, WD = W.A, W.B, W.C, W.D
+
     # orthogonalize second row against last
     Wi = W[end, 1, 1, end]
     WI = removeunit(Wi, 4)
-    @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p']
+    @plansor t[l; r] := conj(WI[r p; p']) * WB[l p; p']
     # TODO: the following is currently broken due to a TensorKit bug
-    # @plansor B′[l p; p'] := W.B[l p; p'] - WI[r p; p'] * t[l; r]
+    # @plansor B′[l p; p'] := WB[l p; p'] - WI[r p; p'] * t[l; r]
     @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r]
-    add!(B′, W.B)
+    add!(B′, WB)
 
     # LQ of second row
     if size(W, 4) == 1
@@ -115,24 +124,24 @@ function right_canonicalize!(
 
         if dim(R) == 0
             V = _rightunit ⊞ _rightunit
-            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
-            Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P)
+            Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
+            Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P)
         else
             V = ⊞(_rightunit, space(Q, 1), _rightunit)
             scale!(Q, d)
             scale!(R, inv(d))
-            Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}())
+            Q1 = typeof(WA)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}())
             Q2 = transpose(Q, ((1, 3), (2,)))
         end
-        H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, W.C, W.D)
+        H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD)
     else
         B′′ = insertleftunit(B′, 4)
-        tmp = transpose(cat(B′′, W.A; dims = 4), ((1,), (3, 4, 2)))
+        tmp = transpose(cat(B′′, WA; dims = 4), ((1,), (3, 4, 2)))
         R, Q = right_orth!(tmp; alg)
         if dim(R) == 0
             V = _rightunit ⊞ _rightunit
-            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
-            Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P)
+            Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}())
+            Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P)
         else
             scale!(Q, d)
             scale!(R, inv(d))
@@ -141,33 +150,34 @@ function right_canonicalize!(
             Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
             V = ⊞(_rightunit, left_virtualspace(Q′), _rightunit)
         end
-        H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, W.C, W.D)
+        H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD)
     end
 
     # absorb into previous site
     W′ = H[i - 1]
+    W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D
 
     if size(W′, 1) > 1 && dim(R) != 0
-        @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r]
+        @plansor A′[l p; p' r] := W′A[l p; p' r'] * R[r'; r]
     else
-        A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ left_virtualspace(H[i].A))
+        A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(H[i].A))
     end
 
     if size(W′, 1) > 1
-        @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r]
-        B′ = add!(removeunit(B′, 4), W′.B)
+        @plansor B′[l p; p' r] := W′A[l p; p' r'] * t[r'; r]
+        B′ = add!(removeunit(B′, 4), W′B)
     else
-        B′ = W′.B
+        B′ = W′B
     end
 
     if dim(R) != 0
-        @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
+        @plansor C′[p; p' r] := W′C[p; p' r'] * R[r'; r]
     else
-        C′ = similar(W′.C, codomain(W′.C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A))
+        C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A))
     end
 
-    @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r]
-    D′ = add!(removeunit(D′, 3), W′.D)
+    @plansor D′[p; p' r] := W′C[p; p' r'] * t[r'; r]
+    D′ = add!(removeunit(D′, 3), W′D)
     H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′)
     return H
 end

From b2b12c95667011d7e9058baf9c92e4b032e94473 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 25 Jun 2026 11:25:27 -0400
Subject: [PATCH 6/9] Densify JordanMPOTensor to TensorMap and key scalars by
 CartesianIndex{4}

- SparseBlockTensorMap(::JordanMPOTensor) now produces a uniform concrete
  `TensorMap` block type with the identity scalars folded in, instead of an
  abstract eltype carrying BraidingTensors.
- Store `scalars` as `Dict{CartesianIndex{4},T}` keyed by the full
  `(row, 1, 1, col)` virtual position, matching `tensors`; update getindex,
  setindex!, haskey, the nonzero_* iterators, the block accessors, scale!,
  isidentitylevel and the Adapt extension accordingly.
- `_identity_tensor` now takes the CartesianIndex{4} position directly.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 ext/MPSKitAdaptExt.jl            |  2 +-
 src/operators/jordanmpotensor.jl | 61 ++++++++++++++------------------
 src/operators/mpohamiltonian.jl  |  6 ++--
 3 files changed, 31 insertions(+), 38 deletions(-)

diff --git a/ext/MPSKitAdaptExt.jl b/ext/MPSKitAdaptExt.jl
index 9d550c75d..912b46569 100644
--- a/ext/MPSKitAdaptExt.jl
+++ b/ext/MPSKitAdaptExt.jl
@@ -34,7 +34,7 @@ end
 @inline function Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor)
     tensors = adapt(to, W.tensors)
     T = scalartype(tensors)
-    scalars = Dict{CartesianIndex{2}, T}(K => convert(T, c) for (K, c) in W.scalars)
+    scalars = Dict{CartesianIndex{4}, T}(K => convert(T, c) for (K, c) in W.scalars)
     return MPSKit.JordanMPOTensor(tensors, scalars)
 end
 @inline Adapt.adapt_structure(to, mpo::MPOHamiltonian) =
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index 253055a26..c0e5cb9b8 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -22,8 +22,8 @@ Rather than storing the dense block matrix, the genuine operators and the identi
 
 - `tensors::SparseBlockTensorMap` holds the non-identity operators over the *full* virtual
   space (so `A`, `B`, `C` and `D` all live at their `(row, 1, 1, col)` position).
-- `scalars::Dict{CartesianIndex{2},T}` holds the scalar multiples of the identity, keyed by
-  their `(row, col)` virtual position; the diagonal corner `1`s are stored here as well.
+- `scalars::Dict{CartesianIndex{4},T}` holds the scalar multiples of the identity, keyed by
+  their `(row, 1, 1, col)` virtual position; the diagonal corner `1`s are stored here as well.
 
 An index is never present in both `tensors` and `scalars`.
 This keeps the ubiquitous identity blocks free of dense storage and lets identities be materialized lazily only when needed.
@@ -43,12 +43,12 @@ struct JordanMPOTensor{
         T <: Number, S, A <: DenseVector{T},
     } <: AbstractBlockTensorMap{T, S, 2, 2}
     tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4}
-    scalars::Dict{CartesianIndex{2}, T}
+    scalars::Dict{CartesianIndex{4}, T}
 
     # constructor from fields
     function JordanMPOTensor{T, S, A}(
             tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4},
-            scalars::Dict{CartesianIndex{2}, T}
+            scalars::Dict{CartesianIndex{4}, T}
         ) where {T, S, A}
         return new{T, S, A}(tensors, scalars)
     end
@@ -58,17 +58,17 @@ struct JordanMPOTensor{
             ::UndefInitializer, V::TensorMapSumSpace{S, 2, 2}
         ) where {T, S, A}
         tensors = SparseBlockTensorMap{tensormaptype(S, 2, 2, A)}(undef, V)
-        scalars = Dict{CartesianIndex{2}, T}()
+        scalars = Dict{CartesianIndex{4}, T}()
         rows, cols = size(tensors, 1), size(tensors, 4)
-        cols > 1 && (scalars[CartesianIndex(1, 1)] = one(T))
-        rows > 1 && (scalars[CartesianIndex(rows, cols)] = one(T))
+        cols > 1 && (scalars[CartesianIndex(1, 1, 1, 1)] = one(T))
+        rows > 1 && (scalars[CartesianIndex(rows, 1, 1, cols)] = one(T))
         return new{T, S, A}(tensors, scalars)
     end
 end
 
 function JordanMPOTensor(
         tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4},
-        scalars::Dict{CartesianIndex{2}, T}
+        scalars::Dict{CartesianIndex{4}, T}
     ) where {T, S, A}
     return JordanMPOTensor{T, S, A}(tensors, scalars)
 end
@@ -175,9 +175,9 @@ function Base.getproperty(W::JordanMPOTensor, sym::Symbol)
     return getfield(W, sym)
 end
 
-# materialize an identity scalar at virtual position `(i, j)` as a dense `TensorMap`
-@inline function _identity_tensor(W::JordanMPOTensor, i::Int, j::Int, c = one(scalartype(W)))
-    τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[i, 1, 1, j])
+# materialize an identity scalar at virtual position `K` as a dense `TensorMap`
+@inline function _identity_tensor(W::JordanMPOTensor, K::CartesianIndex{4}, c = one(scalartype(W)))
+    τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[K])
     return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c)
 end
 
@@ -193,7 +193,7 @@ function _jordan_A(W::JordanMPOTensor)
         A[I[1] - 1, 1, 1, I[4] - 1] = v
     end
     for (K, c) in W.scalars
-        i, j = K[1], K[2]
+        i, j = K[1], K[4]
         (1 < i < rows && 1 < j < cols) || continue
         τ = Tτ(eachspace(A)[i - 1, 1, 1, j - 1])
         A[i - 1, 1, 1, j - 1] = isone(c) ? τ : scale!(TensorMap(τ), c)
@@ -210,8 +210,8 @@ function _jordan_B(W::JordanMPOTensor)
         B[I[1] - 1, 1, 1] = removeunit(v, 4)
     end
     for (K, c) in W.scalars
-        (1 < K[1] < rows && K[2] == cols) || continue
-        B[K[1] - 1, 1, 1] = removeunit(_identity_tensor(W, K[1], K[2], c), 4)
+        (1 < K[1] < rows && K[4] == cols) || continue
+        B[K[1] - 1, 1, 1] = removeunit(_identity_tensor(W, K, c), 4)
     end
     return B
 end
@@ -225,8 +225,8 @@ function _jordan_C(W::JordanMPOTensor)
         C[1, 1, I[4] - 1] = removeunit(v, 1)
     end
     for (K, c) in W.scalars
-        (K[1] == 1 && 1 < K[2] < cols) || continue
-        C[1, 1, K[2] - 1] = removeunit(_identity_tensor(W, K[1], K[2], c), 1)
+        (K[1] == 1 && 1 < K[4] < cols) || continue
+        C[1, 1, K[4] - 1] = removeunit(_identity_tensor(W, K, c), 1)
     end
     return C
 end
@@ -244,7 +244,7 @@ end
 
 function Base.haskey(W::JordanMPOTensor, I::CartesianIndex{4})
     Base.checkbounds(W, I.I...)
-    return haskey(W.scalars, CartesianIndex(I[1], I[4])) || haskey(W.tensors, I)
+    return haskey(W.scalars, I) || haskey(W.tensors, I)
 end
 
 Base.parent(W::JordanMPOTensor) = parent(SparseBlockTensorMap(W))
@@ -254,17 +254,12 @@ BlockTensorKit.issparse(W::JordanMPOTensor) = true
 # Converters
 # ----------
 function BlockTensorKit.SparseBlockTensorMap(W::JordanMPOTensor)
-    W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W), spacetype(W), 2, 2}}(
-        undef_blocks, space(W)
-    )
+    W′ = SparseBlockTensorMap{eltype(W)}(undef_blocks, space(W))
     for (I, v) in nonzero_pairs(W.tensors)
         W′[I] = v
     end
     for (K, c) in W.scalars
-        τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(
-            eachspace(W)[K[1], 1, 1, K[2]]
-        )
-        W′[K[1], 1, 1, K[2]] = isone(c) ? τ : scale!(TensorMap(τ), c)
+        W′[K] = _identity_tensor(W, K, c)
     end
     return W′
 end
@@ -290,9 +285,10 @@ end
 @propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4})
     @assert I[2] == I[3] == 1
     i, j = I[1], I[4]
-    c = get(W.scalars, CartesianIndex(i, j), nothing)
+    K = CartesianIndex(i, 1, 1, j)
+    c = get(W.scalars, K, nothing)
     if c !== nothing
-        return _identity_tensor(W, i, j, c)
+        return _identity_tensor(W, K, c)
     elseif haskey(W.tensors, CartesianIndex(i, 1, 1, j))
         return W.tensors[i, 1, 1, j]
     else
@@ -310,9 +306,9 @@ end
     i, j = I[1], I[4]
     if v isa BraidingTensor
         haskey(W.tensors, CartesianIndex(i, 1, 1, j)) && delete!(W.tensors, CartesianIndex(i, 1, 1, j))
-        W.scalars[CartesianIndex(i, j)] = one(scalartype(W))
+        W.scalars[CartesianIndex(i, 1, 1, j)] = one(scalartype(W))
     else
-        delete!(W.scalars, CartesianIndex(i, j))
+        delete!(W.scalars, CartesianIndex(i, 1, 1, j))
         W.tensors[i, 1, 1, j] = v
     end
     return W
@@ -326,7 +322,7 @@ end
 function BlockTensorKit.nonzero_keys(W::JordanMPOTensor)
     p = collect(CartesianIndex{4}, nonzero_keys(W.tensors))
     for K in keys(W.scalars)
-        push!(p, CartesianIndex(K[1], 1, 1, K[2]))
+        push!(p, K)
     end
     return p
 end
@@ -334,7 +330,7 @@ function BlockTensorKit.nonzero_values(W::JordanMPOTensor)
     return Iterators.flatten(
         (
             nonzero_values(W.tensors),
-            (_identity_tensor(W, K[1], K[2], c) for (K, c) in W.scalars),
+            (_identity_tensor(W, K, c) for (K, c) in W.scalars),
         )
     )
 end
@@ -342,10 +338,7 @@ function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor)
     return Iterators.flatten(
         (
             nonzero_pairs(W.tensors),
-            (
-                CartesianIndex(K[1], 1, 1, K[2]) => _identity_tensor(W, K[1], K[2], c)
-                    for (K, c) in W.scalars
-            ),
+            (K => _identity_tensor(W, K, c) for (K, c) in W.scalars),
         )
     )
 end
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 067bd2a78..fa957acc2 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -652,7 +652,7 @@ function isidentitylevel(H::InfiniteMPOHamiltonian{<:JordanMPOTensor}, i::Int)
         # a diagonal level is an identity level iff every site stores a unit identity
         # scalar there; pure identities live in `scalars` (genuine/scaled operators do not)
         return all(parent(H)) do W
-            c = get(W.scalars, CartesianIndex(i, i), nothing)
+            c = get(W.scalars, CartesianIndex(i, 1, 1, i), nothing)
             return c !== nothing && isone(c)
         end
     end
@@ -794,7 +794,7 @@ function VectorInterface.scale!(
             I[1] == 1 && scale!(v, λ)
         end
         for K in collect(keys(W.scalars))
-            (K[1] == 1 && K[2] != 1) && (W.scalars[K] *= λ)
+            (K[1] == 1 && K[4] != 1) && (W.scalars[K] *= λ)
         end
     end
     return H
@@ -810,7 +810,7 @@ function VectorInterface.scale!(
         end
         empty!(Wd.scalars)
         for (K, c) in Ws.scalars
-            Wd.scalars[K] = (K[1] == 1 && K[2] != 1) ? c * λ : c
+            Wd.scalars[K] = (K[1] == 1 && K[4] != 1) ? c * λ : c
         end
     end
     return Hdst

From f39d1dfc67d9231def26cec8c02ec8cf4d53c862 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 25 Jun 2026 11:31:06 -0400
Subject: [PATCH 7/9] Make JordanMPOTensor(::SparseBlockTensorMap)
 folding-aware

The corner assertion now accepts entries that are approximately identity
braidings (not only literal `BraidingTensor`s), via the `_isbraiding` helper.
The corners are then folded back into scalar identities by the constructor, so
`JordanMPOTensor(SparseBlockTensorMap(W))` round-trips even after the converter
densifies identities into `TensorMap`s.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/operators/jordanmpotensor.jl | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index c0e5cb9b8..b65613e51 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -134,8 +134,18 @@ function JordanMPOTensor(
     )
 end
 
+# whether `t` is (approximately) an identity braiding tensor, so it can be folded back
+# into a scalar identity
+function _isbraiding(t::AbstractTensorMap)
+    t isa BraidingTensor && return true
+    τ = BraidingTensor{scalartype(t), spacetype(t), storagetype(t)}(space(t))
+    return isapprox(t, τ)
+end
+
 function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E, S}
-    @assert W[1, 1, 1, 1] isa BraidingTensor && W[end, 1, 1, end] isa BraidingTensor
+    # the diagonal corners must be (close to) identity braidings; the constructor below
+    # folds them back into scalar identities
+    @assert _isbraiding(W[1, 1, 1, 1]) && _isbraiding(W[end, 1, 1, end])
     # @assert all(I -> I[1] ≤ I[4], nonzero_keys(W))
 
     A = W[2:(end - 1), 1, 1, 2:(end - 1)]

From c43594115274b2fa6423673911d2c7c4b15cd44d Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 25 Jun 2026 12:11:01 -0400
Subject: [PATCH 8/9] Fold A-block identities and rework JordanMPOTensor
 indexing

- `_jordan_A` now returns a concrete `SparseBlockTensorMap{TensorMap}` with the
  identity scalars folded in, instead of `Union{TA,BraidingTensor}` carrying
  BraidingTensors. All four block accessors now share the same TensorMap eltype.
- getindex/setindex! are now written on `CartesianIndex{4}` as the primary form;
  the `Vararg{Int,4}` and linear-index methods delegate to it.
- setindex! accepts a `Number` or a `BraidingTensor`, both of which are stored as
  a scalar identity entry (the number as its coefficient, the braiding as one);
  any other tensor is stored in `tensors`.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/operators/jordanmpotensor.jl | 50 +++++++++++++++-----------------
 1 file changed, 23 insertions(+), 27 deletions(-)

diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index b65613e51..ab6c3e27b 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -194,19 +194,16 @@ end
 # reduced-leg blocks reconstructed from `tensors` and `scalars`
 function _jordan_A(W::JordanMPOTensor)
     rows, cols = size(W, 1), size(W, 4)
-    Tτ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}
     TA = tensormaptype(spacetype(W), 2, 2, storagetype(W))
     VA = space(eachspace(W.tensors)[2:(end - 1), 1, 1, 2:(end - 1)])
-    A = SparseBlockTensorMap{Union{TA, Tτ}}(undef, VA)
+    A = SparseBlockTensorMap{TA}(undef, VA)
     for (I, v) in nonzero_pairs(W.tensors)
         (1 < I[1] < rows && 1 < I[4] < cols) || continue
         A[I[1] - 1, 1, 1, I[4] - 1] = v
     end
     for (K, c) in W.scalars
-        i, j = K[1], K[4]
-        (1 < i < rows && 1 < j < cols) || continue
-        τ = Tτ(eachspace(A)[i - 1, 1, 1, j - 1])
-        A[i - 1, 1, 1, j - 1] = isone(c) ? τ : scale!(TensorMap(τ), c)
+        (1 < K[1] < rows && 1 < K[4] < cols) || continue
+        A[K[1] - 1, 1, 1, K[4] - 1] = _identity_tensor(W, K, c)
     end
     return A
 end
@@ -291,41 +288,40 @@ end
 # Indexing
 # --------
 
-@inline Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4}) = W[I.I...]
-@propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4})
+@inline Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4}) = W[CartesianIndex(I)]
+@propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4})
     @assert I[2] == I[3] == 1
-    i, j = I[1], I[4]
-    K = CartesianIndex(i, 1, 1, j)
-    c = get(W.scalars, K, nothing)
+    c = get(W.scalars, I, nothing)
     if c !== nothing
-        return _identity_tensor(W, K, c)
-    elseif haskey(W.tensors, CartesianIndex(i, 1, 1, j))
-        return W.tensors[i, 1, 1, j]
+        return _identity_tensor(W, I, c)
+    elseif haskey(W.tensors, I)
+        return W.tensors[I]
     else
-        return zeros(storagetype(W), eachspace(W)[i, 1, 1, j])
+        return zeros(storagetype(W), eachspace(W)[I])
     end
 end
 
-@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::CartesianIndex{4})
-    return setindex!(W, v, I.I...)
+@inline function Base.setindex!(W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::Vararg{Int, 4})
+    return setindex!(W, v, CartesianIndex(I))
 end
+@inline function Base.setindex!(W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::Int)
+    return setindex!(W, v, CartesianIndices(W)[I])
+end
+# a scalar or a `BraidingTensor` is stored as a scalar identity; any other tensor is stored
+# in `tensors`
 @propagate_inbounds function Base.setindex!(
-        W::JordanMPOTensor, v::MPOTensor, I::Vararg{Int, 4}
+        W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::CartesianIndex{4}
     )
     @assert I[2] == I[3] == 1
-    i, j = I[1], I[4]
-    if v isa BraidingTensor
-        haskey(W.tensors, CartesianIndex(i, 1, 1, j)) && delete!(W.tensors, CartesianIndex(i, 1, 1, j))
-        W.scalars[CartesianIndex(i, 1, 1, j)] = one(scalartype(W))
+    if v isa Number || v isa BraidingTensor
+        haskey(W.tensors, I) && delete!(W.tensors, I)
+        W.scalars[I] = v isa Number ? convert(scalartype(W), v) : one(scalartype(W))
     else
-        delete!(W.scalars, CartesianIndex(i, 1, 1, j))
-        W.tensors[i, 1, 1, j] = v
+        delete!(W.scalars, I)
+        W.tensors[I] = v
     end
     return W
 end
-@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::Int)
-    return setindex!(W, v, CartesianIndices(W)[I])
-end
 
 # Sparse functionality
 # --------------------

From dda06e0de41e76e2fc9757302a59a5d857330028 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 25 Jun 2026 14:46:08 -0400
Subject: [PATCH 9/9] some more clenaup

---
 src/operators/jordanmpotensor.jl | 38 ++++++++++++++++----------------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index ab6c3e27b..8dd20f6fb 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -134,8 +134,6 @@ function JordanMPOTensor(
     )
 end
 
-# whether `t` is (approximately) an identity braiding tensor, so it can be folded back
-# into a scalar identity
 function _isbraiding(t::AbstractTensorMap)
     t isa BraidingTensor && return true
     τ = BraidingTensor{scalartype(t), spacetype(t), storagetype(t)}(space(t))
@@ -143,19 +141,17 @@ function _isbraiding(t::AbstractTensorMap)
 end
 
 function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E, S}
-    # the diagonal corners must be (close to) identity braidings; the constructor below
-    # folds them back into scalar identities
+    # the diagonal corners must be (close to) identity braidings
     @assert _isbraiding(W[1, 1, 1, 1]) && _isbraiding(W[end, 1, 1, end])
     # @assert all(I -> I[1] ≤ I[4], nonzero_keys(W))
 
-    A = W[2:(end - 1), 1, 1, 2:(end - 1)]
-    B = W[2:(end - 1), 1, 1, end]
-    C = W[1, 1, 1, 2:(end - 1)]
-    D = W[1, 1, 1, end:end] # ensure still blocktensor to allow for sparse
-
-    return JordanMPOTensor(
-        space(W), A, removeunit(B, 4), removeunit(C, 1), removeunit(removeunit(D, 4), 1)
-    )
+    O = jordanmpotensortype(S, storagetype(W))(undef, space(W))
+    corners = (CartesianIndex(1, 1, 1, 1), CartesianIndex(size(W, 1), 1, 1, size(W, 4)))
+    for (I, v) in nonzero_pairs(W)
+        I in corners && continue # already set as identity scalars by the constructor
+        O[I] = v
+    end
+    return O
 end
 
 function jordanmpotensortype(::Type{S}, ::Type{E}) where {S <: VectorSpace, E}
@@ -380,11 +376,15 @@ end
 
 function _conj_mpo(W::JordanMPOTensor)
     V = left_virtualspace(W)' ⊗ physicalspace(W) ← physicalspace(W) ⊗ right_virtualspace(W)'
-    A = _conj_mpo(W.A)
-    @plansor B[-1 -2; -3] ≔ conj(W.B[-1 -3; -2])
-    @plansor C[-1; -2 -3] ≔ conj(W.C[-2; -1 -3])
-    D = copy(adjoint(W.D))
-    return JordanMPOTensor(V, A, B, C, D)
+    O = jordanmpotensortype(spacetype(W), storagetype(W))(undef, V)
+    for (I, v) in nonzero_pairs(W.tensors)
+        O.tensors[I] = _conj_mpo(v)
+    end
+    empty!(O.scalars)
+    for (K, c) in W.scalars
+        O.scalars[K] = conj(c)
+    end
+    return O
 end
 
 function add_physical_charge(O::JordanMPOTensor, charge::Sector)
@@ -427,6 +427,6 @@ function Base.showarg(io::IO, W::JordanMPOTensor, toplevel::Bool)
     return nothing
 end
 
-function TensorKit.type_repr(::Type{<:JordanMPOTensor{E, S}}) where {E, S}
-    return "JordanMPOTensor{$E, " * TensorKit.type_repr(S) * ", …}"
+function TensorKit.type_repr(::Type{JordanMPOTensor{E, S, A}}) where {E, S, A}
+    return "JordanMPOTensor{$E, " * TensorKit.type_repr(S) * ", $A}"
 end