QuantumKitHub · Yue-Zhengyuan · Apr 6, 2026 · Apr 6, 2026
diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl
@@ -19,11 +19,9 @@ Construct the tensor
     └-----------------------------------┘
 ```
 """
-function _tensor_Ra(
-        benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1}
-    ) where {T <: Number, S <: ElementarySpace}
+function _tensor_Ra(benv::BondEnv, b::MPSTensor)
     return @autoopt @tensor Ra[DX1 Db1; DX0 Db0] := (
-        benv[DX1 DY1; DX0 DY0] * b[Db0 DY0; db] * conj(b[Db1 DY1; db])
+        benv[DX1 DY1; DX0 DY0] * b[Db0 db; DY0] * conj(b[Db1 db; DY1])
     )
 end
 
@@ -44,10 +42,10 @@ Construct the tensor
 ```
 """
 function _tensor_Sa(
-        benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2}
+        benv::BondEnv, b::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2}
     ) where {T <: Number, S <: ElementarySpace}
-    return @autoopt @tensor Sa[DX1 Db1; da] := (
-        benv[DX1 DY1; DX0 DY0] * conj(b[Db1 DY1; db]) * a2b2[DX0 DY0; da db]
+    return @autoopt @tensor Sa[DX1 da; Db1] := (
+        benv[DX1 DY1; DX0 DY0] * conj(b[Db1 db; DY1]) * a2b2[DX0 DY0; da db]
     )
 end
 
@@ -67,11 +65,9 @@ Construct the tensor
     └-----------------------------------┘
 ```
 """
-function _tensor_Rb(
-        benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1}
-    ) where {T <: Number, S <: ElementarySpace}
+function _tensor_Rb(benv::BondEnv, a::MPSTensor)
     return @autoopt @tensor Rb[Da1 DY1; Da0 DY0] := (
-        benv[DX1 DY1; DX0 DY0] * a[DX0 Da0; da] * conj(a[DX1 Da1; da])
+        benv[DX1 DY1; DX0 DY0] * a[DX0 da; Da0] * conj(a[DX1 da; Da1])
     )
 end
 
@@ -92,10 +88,10 @@ Construct the tensor
 ```
 """
 function _tensor_Sb(
-        benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2}
+        benv::BondEnv, a::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2}
     ) where {T <: Number, S <: ElementarySpace}
-    return @autoopt @tensor Sb[Da1 DY1; db] := (
-        benv[DX1 DY1; DX0 DY0] * conj(a[DX1 Da1; da]) * a2b2[DX0 DY0; da db]
+    return @autoopt @tensor Sb[Da1 db; DY1] := (
+        benv[DX1 DY1; DX0 DY0] * conj(a[DX1 da; Da1]) * a2b2[DX0 DY0; da db]
     )
 end
 
@@ -116,30 +112,10 @@ Calculate the inner product <a1,b1|a2,b2>
 ```
 """
 function inner_prod(
-        benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2}
+        benv::BondEnv, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2}
     ) where {T <: Number, S <: ElementarySpace}
     return @autoopt @tensor benv[DX1 DY1; DX0 DY0] *
-        conj(a1b1[DX1 DY1; da db]) *
-        a2b2[DX0 DY0; da db]
-end
-
-"""
-$(SIGNATURES)
-
-Calculate the fidelity between two evolution steps
-```
-        |⟨a1,b1|a2,b2⟩|^2
-    --------------------------
-    ⟨a1,b1|a1,b1⟩⟨a2,b2|a2,b2⟩
-```
-"""
-function fidelity(
-        benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2}
-    ) where {T <: Number, S <: ElementarySpace}
-    b12 = inner_prod(benv, a1b1, a2b2)
-    b11 = inner_prod(benv, a1b1, a1b1)
-    b22 = inner_prod(benv, a2b2, a2b2)
-    return abs2(b12) / abs(b11 * b22)
+        conj(a1b1[DX1 DY1; da db]) * a2b2[DX0 DY0; da db]
 end
 
 """
@@ -153,56 +129,53 @@ Contract the axis between `a` and `b` tensors
 ```
 """
 function _combine_ab(
-        a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 1, 2}
+        a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2}
     ) where {T <: Number, S <: ElementarySpace}
     return @tensor ab[DX DY; da db] := a[DX da; D] * b[D; db DY]
 end
-function _combine_ab(
-        a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 2, 1}
-    ) where {T <: Number, S <: ElementarySpace}
-    return @tensor ab[DX DY; da db] := a[DX D; da] * b[D DY; db]
+function _combine_ab(a::MPSTensor, b::MPSTensor)
+    return @tensor ab[DX DY; da db] := a[DX da; D] * b[D db; DY]
 end
 
 """
 $(SIGNATURES)
 
 Calculate the cost function
 ```
-    f(a,b)  = ‖ |a1,b1⟩ - |a2,b2⟩ ‖^2
-    = ⟨a1,b1|a1,b1⟩ - 2 Re⟨a1,b1|a2,b2⟩ + ⟨a2,b2|a2,b2⟩
+    f(a,b) = ‖ |ψ1⟩ - |ψ2⟩ ‖^2
+    = ⟨ψ1|benv|ψ1⟩ - 2 Re⟨ψ1|benv|ψ2⟩ + ⟨ψ2|benv|ψ2⟩
+```
+and the fidelity
+```
+        |⟨ψ1|benv|ψ2⟩|²
+    ------------------------
+    ⟨ψ1|benv|ψ1⟩⟨ψ2|benv|ψ2⟩
 ```
 """
-function cost_function_als(
-        benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2}
-    ) where {T <: Number, S <: ElementarySpace}
-    t1 = inner_prod(benv, a1b1, a1b1)
-    t2 = inner_prod(benv, a2b2, a2b2)
-    t3 = inner_prod(benv, a1b1, a2b2)
-    return real(t1) + real(t2) - 2 * real(t3)
+function cost_function_als(benv, ψ1, ψ2)
+    b12 = inner_prod(benv, ψ1, ψ2)
+    b11 = inner_prod(benv, ψ1, ψ1)
+    b22 = inner_prod(benv, ψ2, ψ2)
+    cost = real(b11) + real(b22) - 2 * real(b12)
+    fid = abs2(b12) / abs(b11 * b22)
+    return cost, fid
 end
 
 """
 $(SIGNATURES)
 
-Solve the equations `Rx x = Sx` (x = a, b) with initial guess `x0`
-```
-    ┌---------------------------┐
-    |   ┌----┐                  |
-    └---|    |--- 1 -- x -- 2 --┘
-        |    |         ↓
-        | Rx |        -3
-        |    |
-    ┌---|    |--- -1       -2 --┐
-    |   └----┘                  |
-    └---------------------------┘
-```
-"""
-function _solve_ab(
-        Rx::AbstractTensorMap{T, S, 2, 2},
-        Sx::AbstractTensorMap{T, S, 2, 1},
-        x0::AbstractTensorMap{T, S, 2, 1},
-    ) where {T <: Number, S <: ElementarySpace}
-    f(x) = (@tensor Sx2[-1 -2; -3] := Rx[-1 -2; 1 2] * x[1 2; -3])
-    x1, info = linsolve(f, Sx, x0, 0, 1)
+Solve the equations `Rx x = Sx` with initial guess `x0`.
+"""
+function _solve_als(
+        Rx::AbstractTensorMap{T, S, N, N},
+        Sx::GenericMPSTensor{S, N},
+        x0::GenericMPSTensor{S, N}; kwargs...
+    ) where {T, S, N}
+    @assert N >= 2
+    pR = (codomainind(Rx), domainind(Rx))
+    pX = ((1, (3:(N + 1))...), (2,))
+    pRX = ((1, N + 1, (2:(N - 1))...), (N,))
+    f(x) = tensorcontract(Rx, pR, false, x, pX, false, pRX)
+    x1, info = linsolve(f, Sx, x0, 0, 1; kwargs...)
     return x1, info
 end
diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl
@@ -1,14 +1,16 @@
 """
 Construct the environment (norm) tensor
 ```
-    C1---T1---------T1---C2   r-1
-    |    ‖          ‖    |
-    T4===XX==     ==YY===T2    r
-    |    ‖          ‖    |
-    C4---T3---------T3---C3   r+1
-    c-1  c         c+1  c+2
+    -1  C1---E1---------E1---C2
+        |    ‖          ‖    |
+    0   E4===XX==     ==YY===E2
+        |    ‖          ‖    |
+    1   C4---E3---------E3---C3
+        -1   0          1    2
 ```
-where `XX = X' X` and `YY = Y' Y` (stacked together).
+with `X` at position `[row, col]`.
+`X, Y` are unitary tensors produced when finding the reduced site tensors 
+with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together).
 
 Axis order: `[DX1 DY1; DX0 DY0]`, as in
 ```
@@ -21,7 +23,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in
     └---------------------┘
 ```
 """
-function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv)
+function bondenv_ctm(
+        row::Int, col::Int, X::T, Y::T, env::CTMRGEnv
+    ) where {T}
     Nr, Nc = size(env.corners)[[2, 3]]
     cm1 = _prev(col, Nc)
     cp1 = _next(col, Nc)
@@ -32,29 +36,29 @@ function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv
     c2 = env.corners[2, rm1, cp2]
     c3 = env.corners[3, rp1, cp2]
     c4 = env.corners[4, rp1, cm1]
-    t1X, t1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1]
-    t2 = env.edges[2, row, cp2]
-    t3X, t3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1]
-    t4 = env.edges[4, row, cm1]
+    e1X, e1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1]
+    e2 = env.edges[2, row, cp2]
+    e3X, e3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1]
+    e4 = env.edges[4, row, cm1]
     #= index labels
 
-    C1--χ4--T1X---------χ6---------T1Y--χ8---C2     r-1
+    C1--χ4--E1X---------χ6---------E1Y--χ8---C2     r-1
     |        ‖                      ‖        |
     χ2      DNX                    DNY      χ10
     |        ‖                      ‖        |
-    T4==DWX==XX===DX==       ==DY===YY==DEY==T2     r
+    E4==DWX==XX===DX==       ==DY===YY==DEY==E2     r
     |        ‖                      ‖        |
     χ1      DSX                    DSY       χ9
     |        ‖                      ‖        |
-    C4--χ3--T3X---------χ5---------T3Y--χ7---C3     r+1
+    C4--χ3--E3X---------χ5---------E3Y--χ7---C3     r+1
     c-1      c                      c+1     c+2
     =#
+    X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y)
     @autoopt @tensor benv[DX1 DY1; DX0 DY0] :=
-        c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] *
-        X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] *
-        c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] *
-        Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8]
-
+        c4[χ3 χ1] * e4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * e3X[χ5 DSX0 DSX1 χ3] *
+        twistdual(X, 1)[pX DNX0 DX0 DSX0 DWX0] * conj(X[pX DNX1 DX1 DSX1 DWX1]) * e1X[χ4 DNX0 DNX1 χ6] *
+        c3[χ9 χ7] * e2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * e3Y[χ7 DSY0 DSY1 χ5] *
+        twistdual(Y, 1)[pY DNY0 DEY0 DSY0 DY0] * conj(Y[pY DNY1 DEY1 DSY1 DY1]) * e1Y[χ6 DNY0 DNY1 χ8]
     normalize!(benv, Inf)
     return benv
 end
diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl
@@ -1,3 +1,13 @@
 const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace}
+const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace}
+const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace}
+# Orthogonal tensors obtained PEPSTensor/PEPOTensor
+# with one physical leg factored out by `_qr_bond`
 const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace}
 const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace}
+
+"Convert tensor `t` connected by the bond to be truncated to a `PEPSTensor`."
+_prepare_site_tensor(t::PEPSTensor) = t
+_prepare_site_tensor(t::PEPOTensor) = first(fuse_physicalspaces(t))
+_prepare_site_tensor(t::PEPSOrth) = permute(insertleftunit(t, 1), ((1,), (2, 3, 4, 5)))
+_prepare_site_tensor(t::PEPOOrth) = permute(t, ((1,), (2, 3, 4, 5)))
diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl
@@ -1,5 +1,5 @@
 """
-Replace bond environment `benv` by its positive approximant `Z† Z`
+Find positive approximant `Z† Z` of a norm tensor `benv`
 (returns the "half environment" `Z`)
 ```
     ┌-----------------┐     ┌---------------┐
@@ -11,9 +11,9 @@ Replace bond environment `benv` by its positive approximant `Z† Z`
     └-----------------┘     └---------------┘
 ```
 """
-function positive_approx(benv::BondEnv)
+function positive_approx(benv::AbstractTensorMap{T, S, N, N}) where {T, S, N}
     # eigen-decomposition: benv = U D U'
-    D, U = eigh_full((benv + benv') / 2)
+    D, U = eigh_full!(project_hermitian(benv))
     # determine if `env` is (mostly) positive or negative
     sgn = sign(sum(D.data))
     # When optimizing the truncation of a bond,