diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index e3935c756..ffee4ebb7 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -471,8 +471,9 @@ function swap!(
     O₁, O₂ = mpo[i], mpo[i + 1]
 
     if inv
+        # `τ'` must not be the first tensor here (QuantumKitHub/TensorKit.jl#622)
         @plansor O₂₁[-1 -2 -3; -4 -5 -6] :=
-            τ'[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ[2 3; -1 -4]
+            O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ[2 3; -1 -4] * τ'[-3 -6; 4 5]
     else
         @plansor O₂₁[-1 -2 -3; -4 -5 -6] :=
             τ[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ'[2 3; -1 -4]
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index c3f1c4b7d..c3f40e129 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -516,7 +516,9 @@ function _instantiate_operators(lattice, local_operators)
     # A comprehension rather than `map`: all of the work happens in here, and type inference
     # limits `map` calls that are nested within another `map` (JuliaLang/julia#64862).
     return [
-        instantiate_operator(lattice, _decompose_once!(representatives, term))
+        Base.inferencebarrier(instantiate_operator)(
+            lattice, _decompose_once!(representatives, term)
+        )
             for term in collect(local_operators)
     ]
 end
@@ -538,19 +540,23 @@ function _decompose_once!(representatives, O::AbstractTensorMap)
     # cached so that it is not recomputed for every term
     norm_squared_new = real(inner(O, O))
     for (O_rep, norm_squared_rep, mpo) in representatives
-        λ = _proportionality(
+        λ = Base.inferencebarrier(_proportionality)(
             O, O_rep; norm_squared_new, norm_squared_old = norm_squared_rep
         )
         isnothing(λ) && continue
         isone(λ) && return mpo
-        Os = parent(mpo)
-        return FiniteMPO([i == lastindex(Os) ? λ * Os[i] : Os[i] for i in eachindex(Os)])
+        return Base.inferencebarrier(_rescale_last)(mpo, λ)
     end
     mpo = FiniteMPO(O)
     push!(representatives, (O, norm_squared_new, mpo))
     return mpo
 end
 
+function _rescale_last(mpo, λ)
+    Os = parent(mpo)
+    return FiniteMPO([i == lastindex(Os) ? λ * Os[i] : Os[i] for i in eachindex(Os)])
+end
+
 """
     _assign_channels!(nonzero_keys, nonzero_opps, local_mpos)
 
diff --git a/test/mpo/finite.jl b/test/mpo/finite.jl
index b16d8081e..627e5c055 100644
--- a/test/mpo/finite.jl
+++ b/test/mpo/finite.jl
@@ -118,3 +118,20 @@ end
         @test convert(TensorMap, mpo) ≈ O
     end
 end
+
+@testset "swap" begin
+    Vs = (
+        ℂ^2, Vect[FermionParity](0 => 1, 1 => 1),
+        Vect[FibonacciAnyon](:I => 1, :τ => 1),
+    )
+    for V in Vs
+        O = randn(ComplexF64, V ⊗ V ← V ⊗ V)
+        mpo = FiniteMPO(O)
+        # compare the dense operators: `isapprox` on the MPOs is unreliable at its default
+        # tolerance (QuantumKitHub/MPSKit.jl#555)
+        mpo₁ = MPSKit.swap(MPSKit.swap(mpo, 1), 1; inv = true)
+        @test convert(TensorMap, mpo₁) ≈ O
+        mpo₂ = MPSKit.swap(MPSKit.swap(mpo, 1; inv = true), 1)
+        @test convert(TensorMap, mpo₂) ≈ O
+    end
+end