diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index 0b8e27800..b51cf1d75 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -29,6 +29,7 @@ jobs:
           - states
           - operators
           - algorithms
+          - multifusion
           - other
         os:
           - ubuntu-latest
diff --git a/src/algorithms/ED.jl b/src/algorithms/ED.jl
index a3a5c4036..df145b064 100644
--- a/src/algorithms/ED.jl
+++ b/src/algorithms/ED.jl
@@ -1,6 +1,6 @@
 """
     exact_diagonalization(H::FiniteMPOHamiltonian;
-                          sector=first(sectors(oneunit(physicalspace(H, 1)))),
+                          sector=rightunit(H),
                           len::Int=length(H), num::Int=1, which::Symbol=:SR,
                           alg=Defaults.alg_eigsolve(; dynamic_tols=false))
                             -> vals, state_vecs, convhist
@@ -13,7 +13,7 @@ equivalent to dense eigenvectors.
 - `H::FiniteMPOHamiltonian`: the Hamiltonian to diagonalize.
 
 ### Keyword arguments
-- `sector=first(sectors(oneunit(physicalspace(H, 1))))`: the total charge of the
+- `sector=rightunit(H)`: the total charge of the
   eigenvectors, which is chosen trivial by default.
 - `len::Int=length(H)`: the length of the system.
 - `num::Int=1`: the number of eigenvectors to find.
@@ -32,7 +32,7 @@ equivalent to dense eigenvectors.
 """
 function exact_diagonalization(
         H::FiniteMPOHamiltonian;
-        sector = one(sectortype(H)), num::Int = 1, which::Symbol = :SR,
+        sector = rightunit(H), num::Int = 1, which::Symbol = :SR,
         alg = Defaults.alg_eigsolve(; dynamic_tols = false)
     )
     L = length(H)
@@ -44,7 +44,7 @@ function exact_diagonalization(
 
     # fuse from left to right
     ALs = Vector{Union{Missing, TA}}(missing, L)
-    left = oneunit(spacetype(H))
+    left = spacetype(H)(rightunit(sector) => 1)
     for i in 1:(middle_site - 1)
         P = physicalspace(H, i)
         ALs[i] = isomorphism(T, left ⊗ P ← fuse(left ⊗ P))
diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl
index acd3adebf..3f9e24e17 100644
--- a/src/algorithms/correlators.jl
+++ b/src/algorithms/correlators.jl
@@ -16,7 +16,7 @@ function correlator(
     )
     first(js) > i || @error "i should be smaller than j ($i, $(first(js)))"
     S₁ = _firstspace(O₁)
-    S₁ == oneunit(S₁) || throw(ArgumentError("O₁ should start with a trivial leg."))
+    isunitspace(S₁) || throw(ArgumentError("O₁ should start with a trivial leg."))
     S₂ = _lastspace(O₂)
     S₂ == S₁' || throw(ArgumentError("O₂ should end with a trivial leg."))
 
@@ -39,8 +39,9 @@ function correlator(
 end
 
 function correlator(
-        state::AbstractMPS, O₁₂::AbstractTensorMap{<:Any, S, 2, 2}, i::Int, j
+        state::AbstractMPS, O₁₂::AbstractTensorMap{<:Any, S, 2, 2}, i::Int,
+        j
     ) where {S}
-    O₁, O₂ = decompose_localmpo(add_util_leg(O₁₂))
+    O₁, O₂ = decompose_localmpo(add_util_mpoleg(O₁₂))
     return correlator(state, O₁, O₂, i, j)
 end
diff --git a/src/algorithms/excitation/chepigaansatz.jl b/src/algorithms/excitation/chepigaansatz.jl
index e92a618be..153b262f9 100644
--- a/src/algorithms/excitation/chepigaansatz.jl
+++ b/src/algorithms/excitation/chepigaansatz.jl
@@ -35,10 +35,10 @@ end
 
 function excitations(
         H, alg::ChepigaAnsatz, ψ::FiniteMPS, envs = environments(ψ, H);
-        sector = one(sectortype(ψ)), num::Int = 1, pos::Int = length(ψ) ÷ 2
+        sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) || throw(ArgumentError("invalid position $pos"))
-    sector == one(sector) || error("not yet implemented for charged excitations")
+    isunit(sector) || error("not yet implemented for charged excitations")
 
     # add random offset to kickstart Krylov process:
     AC = ψ.AC[pos]
@@ -100,10 +100,10 @@ end
 
 function excitations(
         H, alg::ChepigaAnsatz2, ψ::FiniteMPS, envs = environments(ψ, H);
-        sector = one(sectortype(ψ)), num::Int = 1, pos::Int = length(ψ) ÷ 2
+        sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) - 1 || throw(ArgumentError("invalid position $pos"))
-    sector == one(sector) || error("not yet implemented for charged excitations")
+    isunit(sector) || error("not yet implemented for charged excitations")
 
     # add random offset to kickstart Krylov process:
     @plansor AC2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3]
diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl
index 0eb400258..372a3383c 100644
--- a/src/algorithms/excitation/quasiparticleexcitation.jl
+++ b/src/algorithms/excitation/quasiparticleexcitation.jl
@@ -46,7 +46,6 @@ end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renvs; num::Int = 1)
     E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs)
     H_eff = EffectiveExcitationHamiltonian(H, lenvs, renvs, E)
-
     Es, ϕs, convhist = eigsolve(ϕ₀, num, :SR, alg.alg) do ϕ
         return H_eff(ϕ; alg.alg_environments...)
     end
@@ -86,14 +85,14 @@ Create and optimise infinite quasiparticle states.
 # Keywords
 - `num::Int`: number of excited states to compute
 - `solver`: algorithm for the linear solver of the quasiparticle environments
-- `sector=one(sectortype(left_ψ))`: charge of the quasiparticle state
+- `sector=leftunit(left_ψ)`: charge of the quasiparticle state
 - `parallel=true`: enable multi-threading over different momenta
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, momentum::Number, lmps::InfiniteMPS,
         lenvs = environments(lmps, H), rmps::InfiniteMPS = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), kwargs...
+        sector = leftunit(lmps), kwargs...
     )
     ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector, momentum)
     return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...)
@@ -103,15 +102,27 @@ function excitations(
         lenvs = environments(lmps, H), rmps = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
         verbosity = Defaults.verbosity, num = 1,
-        sector = one(sectortype(lmps)), parallel = true, kwargs...
-    )
-    Toutput = Core.Compiler.return_type(
-        excitations,
-        Tuple{
-            typeof(H), typeof(alg), eltype(momenta), typeof(lmps),
-            typeof(lenvs), typeof(rmps), typeof(renvs),
-        }
+        sector = leftunit(lmps), parallel = true, kwargs...
     )
+    # wrapper to evaluate sector as positional argument
+    Toutput = let
+        function wrapper(H, alg, p, lmps, lenvs, rmps, renvs, sector; kwargs...)
+            return excitations(
+                H, alg, p, lmps, lenvs, rmps, renvs;
+                sector, kwargs...
+            )
+        end
+        Core.Compiler.return_type(
+            wrapper,
+            Tuple{
+                typeof(H), typeof(alg),
+                eltype(momenta), typeof(lmps),
+                typeof(lenvs),
+                typeof(rmps), typeof(renvs), typeof(sector),
+            }
+        )
+    end
+
     results = similar(momenta, Toutput)
     scheduler = parallel ? :greedy : :serial
     tmap!(results, momenta; scheduler) do momentum
@@ -167,13 +178,13 @@ Create and optimise finite quasiparticle states.
 
 # Keywords
 - `num::Int`: number of excited states to compute
-- `sector=one(sectortype(left_ψ))`: charge of the quasiparticle state
+- `sector=leftunit(lmps)`: charge of the quasiparticle state
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, lmps::FiniteMPS,
         lenvs = environments(lmps, H), rmps::FiniteMPS = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), num = 1
+        sector = leftunit(lmps), num = 1
     )
     ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector)
     return excitations(H, alg, ϕ₀, lenvs, renvs; num)
@@ -262,7 +273,7 @@ function excitations(
         H::MultilineMPO, alg::QuasiparticleAnsatz, momentum::Real, lmps::MultilineMPS,
         lenvs = environments(lmps, H), rmps = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), kwargs...
+        sector = leftunit(lmps), kwargs...
     )
     ϕ₀ = LeftGaugedQP(randn, lmps, rmps; sector, momentum)
     return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...)
diff --git a/src/algorithms/fixedpoint.jl b/src/algorithms/fixedpoint.jl
index 81e7c69d9..188c65559 100644
--- a/src/algorithms/fixedpoint.jl
+++ b/src/algorithms/fixedpoint.jl
@@ -4,14 +4,14 @@
     fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec
     fixedpoint(A, x₀, which::Symbol, alg) -> val, vec
 
-Compute the fixedpoint of a linear operator `A` using the specified eigensolver `alg`. The
+Compute the fixed point of a linear operator `A` using the specified eigensolver `alg`. The
 fixedpoint is assumed to be unique.
 """
 function fixedpoint(A, x₀, which::Symbol, alg::Lanczos)
     vals, vecs, info = eigsolve(A, x₀, 1, which, alg)
 
     if info.converged == 0
-        @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
+        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     end
 
     return vals[1], vecs[1]
@@ -21,10 +21,10 @@ function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi)
     TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg)
 
     if info.converged == 0
-        @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
+        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     end
     if size(TT, 2) > 1 && TT[2, 1] != 0
-        @warnv 1 "non-unique fixedpoint detected"
+        @warnv 1 "non-unique fixed point detected"
     end
 
     return vals[1], vecs[1]
diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl
index 478fe339c..2088a6763 100644
--- a/src/algorithms/groundstate/gradient_grassmann.jl
+++ b/src/algorithms/groundstate/gradient_grassmann.jl
@@ -59,7 +59,7 @@ function find_groundstate(
         ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H)
     )::Tuple{S, P, Float64} where {S, P}
     !isa(ψ, FiniteMPS) || dim(ψ.C[end]) == 1 ||
-        @warn "This is not fully supported - split the mps up in a sum of mps's and optimize seperately"
+        @warn "This is not fully supported - split the mps up in a sum of mps's and optimize separately"
     normalize!(ψ)
 
     fg(x) = GrassmannMPS.fg(x, H, envs)
diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl
index dd3c3be9c..3d93a2e0b 100644
--- a/src/algorithms/timestep/taylorcluster.jl
+++ b/src/algorithms/timestep/taylorcluster.jl
@@ -186,14 +186,14 @@ function make_time_mpo(
     H′ = copy(parent(H))
 
     V_left = left_virtualspace(H[1])
-    V_left′ = ⊞(V_left, oneunit(V_left), oneunit(V_left))
+    V_left′ = ⊞(V_left, rightunitspace(V_left), rightunitspace(V_left))
     H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1]))
     for (I, v) in nonzero_pairs(H[1])
         H′[1][I] = v
     end
 
     V_right = right_virtualspace(H[end])
-    V_right′ = ⊞(oneunit(V_right), oneunit(V_right), V_right)
+    V_right′ = ⊞(rightunitspace(V_right), rightunitspace(V_right), V_right)
     H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′)
     for (I, v) in nonzero_pairs(H[end])
         H′[end][I[1], 1, 1, end] = v
diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl
index 7b2fe2169..6ba74c91c 100644
--- a/src/algorithms/timestep/wii.jl
+++ b/src/algorithms/timestep/wii.jl
@@ -94,14 +94,14 @@ function make_time_mpo(
     H′ = copy(parent(H))
 
     V_left = left_virtualspace(H[1])
-    V_left′ = ⊞(V_left, oneunit(V_left), oneunit(V_left))
+    V_left′ = ⊞(V_left, rightunitspace(V_left), rightunitspace(V_left))
     H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1]))
     for (I, v) in nonzero_pairs(H[1])
         H′[1][I] = v
     end
 
     V_right = right_virtualspace(H[end])
-    V_right′ = ⊞(oneunit(V_right), oneunit(V_right), V_right)
+    V_right′ = ⊞(rightunitspace(V_right), rightunitspace(V_right), V_right)
     H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′)
     for (I, v) in nonzero_pairs(H[end])
         H′[end][I[1], 1, 1, end] = v
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index b6f98de64..0689654b2 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -28,7 +28,7 @@ Return the density matrix of the infinite temperature state for a given Hamilton
 This is the identity matrix in the physical space, and the identity in the auxiliary space.
 """
 function infinite_temperature_density_matrix(H::MPOHamiltonian)
-    V = oneunit(spacetype(H))
+    V = first(left_virtualspace(H[1]))
     W = map(1:length(H)) do site
         return BraidingTensor{scalartype(H)}(physicalspace(H, site), V)
     end
@@ -78,7 +78,7 @@ end
 
 """
     transfer_spectrum(above::InfiniteMPS; below=above, tol=Defaults.tol, num_vals=20,
-                           sector=first(sectors(oneunit(left_virtualspace(above, 1)))))
+                           sector=leftunit(above))
 
 Calculate the partial spectrum of the left mixed transfer matrix corresponding to the
 overlap of a given `above` state and a `below` state. The `sector` keyword argument can be
@@ -89,7 +89,7 @@ domain of each eigenvector. The `tol` and `num_vals` keyword arguments are passe
 """
 function transfer_spectrum(
         above::InfiniteMPS; below = above, tol = Defaults.tol, num_vals = 20,
-        sector = first(sectors(oneunit(left_virtualspace(above, 1))))
+        sector = leftunit(above)
     )
     init = randomize!(
         similar(
@@ -275,13 +275,13 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O}, L = length(mpo)) wher
     V_wrap = left_virtualspace(mpo, 1)
     ST = storagetype(O)
 
-    util = isometry(storagetype(O), oneunit(V_wrap) ← one(V_wrap))
+    util = isometry(storagetype(O), rightunitspace(V_wrap) ← one(V_wrap))
     @plansor cup[-1; -2 -3] := id(ST, V_wrap)[-2; -3] * util[-1]
 
     local F_right
     for i in 1:L
-        V_left = i == 1 ? oneunit(V_wrap) : fuse(V_wrap ⊗ left_virtualspace(mpo, i))
-        V_right = i == L ? oneunit(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
+        V_left = i == 1 ? rightunitspace(V_wrap) : fuse(V_wrap ⊗ left_virtualspace(mpo, i))
+        V_right = i == L ? rightunitspace(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
         output[i] = similar(
             mpo[i], V_left * physicalspace(mpo, i) ← physicalspace(mpo, i) * V_right
         )
@@ -336,7 +336,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L = length(H))
     output = Vector{O}(undef, L)
     for site in 1:L
         V_left = if site == 1
-            oneunit(V_wrap)
+            rightunitspace(V_wrap)
         else
             vs = Vector{S}(undef, chi_)
             for (k, v) in indmap
@@ -345,7 +345,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L = length(H))
             SumSpace(vs)
         end
         V_right = if site == L
-            oneunit(V_wrap)
+            rightunitspace(V_wrap)
         else
             vs = Vector{S}(undef, chi_)
             for (k, v) in indmap
diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl
index 1ec89a5ee..14aaacf16 100644
--- a/src/environments/abstract_envs.jl
+++ b/src/environments/abstract_envs.jl
@@ -81,7 +81,7 @@ function environment_alg(
         tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim,
         verbosity = Defaults.VERBOSE_NONE
     )
-    max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1))
+    max_krylovdim = ceil(Int, dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1)))
     return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity)
 end
 function environment_alg(
diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl
index 29f40ac1f..52e6c9f4d 100644
--- a/src/operators/abstractmpo.jl
+++ b/src/operators/abstractmpo.jl
@@ -31,6 +31,9 @@ right_virtualspace(mpo::AbstractMPO) = map(Base.Fix1(right_virtualspace, mpo), e
 physicalspace(mpo::AbstractMPO, site::Int) = physicalspace(mpo[site])
 physicalspace(mpo::AbstractMPO) = map(Base.Fix1(physicalspace, mpo), eachsite(mpo))
 
+TensorKit.leftunit(mpo::AbstractMPO) = only(sectors(leftunitspace(physicalspace(mpo, 1))))
+TensorKit.rightunit(mpo::AbstractMPO) = only(sectors(rightunitspace(physicalspace(mpo, 1))))
+
 for ftype in (:spacetype, :sectortype, :storagetype)
     @eval TensorKit.$ftype(mpo::AbstractMPO) = $ftype(typeof(mpo))
     @eval TensorKit.$ftype(::Type{MPO}) where {MPO <: AbstractMPO} = $ftype(eltype(MPO))
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index 58b7e2ecf..e0c3b04de 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -29,7 +29,7 @@ function FiniteMPO(Os::AbstractVector{O}) where {O}
 end
 
 function FiniteMPO(O::AbstractTensorMap{T, S, N, N}) where {T, S, N}
-    return FiniteMPO(decompose_localmpo(add_util_leg(O)))
+    return FiniteMPO(decompose_localmpo(add_util_mpoleg(O)))
 end
 
 """
@@ -223,8 +223,9 @@ function Base.:*(mpo1::FiniteMPO{<:MPOTensor}, mpo2::FiniteMPO{<:MPOTensor})
     N = check_length(mpo1, mpo2)
     (S = spacetype(mpo1)) == spacetype(mpo2) || throw(SectorMismatch())
 
-    if (left_virtualspace(mpo1, 1) != oneunit(S) || left_virtualspace(mpo2, 1) != oneunit(S)) ||
-            (right_virtualspace(mpo1, N) != oneunit(S) || right_virtualspace(mpo2, N) != oneunit(S))
+
+    if (isunitspace(left_virtualspace(mpo1, 1)) || isunitspace(left_virtualspace(mpo2, 1))) ||
+            (isunitspace(right_virtualspace(mpo1, N)) || isunitspace(right_virtualspace(mpo2, N)))
         @warn "left/right virtual space is not trivial, fusion may not be unique"
         # this is a warning because technically any isomorphism that fuses the left/right
         # would work and for now I dont feel like figuring out if this is important
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index aa7c4004e..937797981 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -80,11 +80,11 @@ function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O <: JordanMPO
     # left end
     nlvls = size(W_mats[1], 1)
     @assert nlvls == 1 "left boundary should have a single level"
-    Vspaces[1] = SumSpace(oneunit(S))
+    Vspaces[1] = SumSpace(unitspace(S))
     # right end
     nlvls = size(W_mats[end], 2)
     @assert nlvls == 1 "right boundary should have a single level"
-    Vspaces[end] = SumSpace(oneunit(S))
+    Vspaces[end] = SumSpace(unitspace(S))
 
     # start filling spaces
     # note that we assume that the FSA does not contain "dead ends", as this would mess with the
@@ -102,7 +102,7 @@ function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O <: JordanMPO
             # start by assuming trivial spaces everywhere -- replace everything that we know
             # assume spacecheck errors will happen when filling the BlockTensors
             nlvls = size(W_mat, 2)
-            Vs_right = SumSpace(fill(oneunit(S), nlvls))
+            Vs_right = SumSpace(fill(unitspace(S), nlvls))
         end
 
         for I in eachindex(IndexCartesian(), W_mat)
@@ -185,7 +185,7 @@ function InfiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O <: MPOTens
     MissingS = Union{Missing, S}
     Vspaces = PeriodicArray([Vector{MissingS}(missing, nlvls) for _ in 1:L])
     for V in Vspaces
-        V[1] = V[end] = oneunit(S)
+        V[1] = V[end] = unitspace(S)
     end
 
     haschanged = true
@@ -235,7 +235,7 @@ function InfiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O <: MPOTens
         end
     end
 
-    foreach(Base.Fix2(replace!, missing => oneunit(S)), Vspaces)
+    foreach(Base.Fix2(replace!, missing => unitspace(S)), Vspaces)
     Vsumspaces = map(Vspaces) do V
         return SumSpace(collect(S, V))
     end
@@ -420,13 +420,17 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_opera
     E = scalartype(T)
     S = spacetype(T)
 
+    # avoid using one(S)
+    somempo = local_mpos[1].second[1]
+    _rightunit = space(somempo, 1) # should be rightunitspace, but MPOHamiltonians are always diagonal for now
+
     virtualsumspaces = Vector{SumSpace{S}}(undef, length(lattice) + 1)
-    virtualsumspaces[1] = SumSpace(fill(oneunit(S), 1))
-    virtualsumspaces[end] = SumSpace(fill(oneunit(S), 1))
+    virtualsumspaces[1] = SumSpace(fill(_rightunit, 1))
+    virtualsumspaces[end] = SumSpace(fill(_rightunit, 1))
 
     for i in 1:(length(lattice) - 1)
         n_channels = maximum(last, nonzero_keys[i]; init = 1) + 1
-        V = SumSpace(fill(oneunit(S), n_channels))
+        V = SumSpace(fill(_rightunit, n_channels))
         if n_channels > 2
             for ((key_L, key_R), O) in zip(nonzero_keys[i], nonzero_opps[i])
                 V[key_R == 0 ? end : key_R] = if O isa Number
@@ -504,9 +508,13 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, local_
     virtualspaces = PeriodicArray(
         [Vector{MissingS}(missing, operator_size) for _ in 1:length(nonzero_keys)]
     )
+    # avoid using one(S)
+    somempo = local_mpos[1].second[1]
+    _rightunit = space(somempo, 1) # should be a rightunitspace
+
     for V in virtualspaces
-        V[1] = oneunit(S)
-        V[end] = oneunit(S)
+        V[1] = _rightunit
+        V[end] = _rightunit
     end
 
     # start by filling in tensors -> space information available
@@ -552,7 +560,7 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, local_
         end
     end
 
-    foreach(Base.Fix2(replace!, missing => oneunit(S)), virtualspaces)
+    foreach(Base.Fix2(replace!, missing => _rightunit), virtualspaces)
     virtualsumspaces = map(virtualspaces) do V
         return SumSpace(collect(S, V))
     end
@@ -683,7 +691,8 @@ function Base.:+(
     ) where {O <: JordanMPOTensor}
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
-    Vtriv = oneunit(spacetype(H₁))
+    # TODO: check if spaces match
+    Vtriv = rightunitspace(first(physicalspace(H₁)))
 
     for i in 1:N
         A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
@@ -707,7 +716,8 @@ function Base.:+(
     ) where {O <: JordanMPOTensor}
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
-    Vtriv = oneunit(spacetype(H₁))
+    # TODO: check if spaces match
+    Vtriv = rightunitspace(first(physicalspace(H₁)))
     for i in 1:N
         A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
         B = cat(H₁[i].B, H₂[i].B; dims = 1)
diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl
index 34ce1b159..607cc9bfc 100644
--- a/src/operators/ortho.jl
+++ b/src/operators/ortho.jl
@@ -22,11 +22,11 @@ function left_canonicalize!(
         Q, R = left_orth!(tmp; alg)
 
         if dim(R) == 0 # fully truncated
-            V = oneunit(S) ⊞ oneunit(S)
+            V = unitspace(S) ⊞ unitspace(S)
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
             Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
         else
-            V = ⊞(oneunit(S), space(R, 1), oneunit(S))
+            V = ⊞(unitspace(S), space(R, 1), unitspace(S))
             scale!(Q, d)
             scale!(R, inv(d))
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1))
@@ -38,7 +38,7 @@ function left_canonicalize!(
         Q, R = left_orth!(tmp; alg)
 
         if dim(R) == 0 # fully truncated
-            V = oneunit(S) ⊞ oneunit(S)
+            V = unitspace(S) ⊞ unitspace(S)
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
             Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
         else
@@ -47,7 +47,7 @@ function left_canonicalize!(
             Q′ = transpose(Q, ((2, 3), (1, 4)))
             Q1 = Q′[2:end, 1, 1, 1]
             Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
-            V = ⊞(oneunit(S), right_virtualspace(Q′), oneunit(S))
+            V = ⊞(unitspace(S), right_virtualspace(Q′), unitspace(S))
         end
         H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D)
     end
@@ -108,11 +108,11 @@ function right_canonicalize!(
         R, Q = right_orth!(tmp; alg)
 
         if dim(R) == 0
-            V = oneunit(S) ⊞ oneunit(S)
+            V = unitspace(S) ⊞ unitspace(S)
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
             Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W))
         else
-            V = ⊞(oneunit(S), space(Q, 1), oneunit(S))
+            V = ⊞(unitspace(S), space(Q, 1), unitspace(S))
             scale!(Q, d)
             scale!(R, inv(d))
             Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
@@ -123,7 +123,7 @@ function right_canonicalize!(
         tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
         R, Q = right_orth!(tmp; alg)
         if dim(R) == 0
-            V = oneunit(S) ⊞ oneunit(S)
+            V = unitspace(S) ⊞ unitspace(S)
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
             Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W))
         else
@@ -132,7 +132,7 @@ function right_canonicalize!(
             Q′ = transpose(Q, ((1, 4), (2, 3)))
             Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
             Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
-            V = ⊞(oneunit(S), left_virtualspace(Q′), oneunit(S))
+            V = ⊞(unitspace(S), left_virtualspace(Q′), unitspace(S))
         end
         H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D)
     end
diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl
index 3fa45832d..4039dd15e 100644
--- a/src/states/abstractmps.jl
+++ b/src/states/abstractmps.jl
@@ -110,7 +110,7 @@ function isfullrank(V::TensorKit.TensorMapSpace; side = :both)
 end
 
 """
-    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defalts.alg_qr())
+    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defaults.alg_qr())
 
 Make the set of MPS tensors full rank by performing a series of orthogonalizations.
 """
@@ -246,3 +246,6 @@ physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ))
 Return an iterator over the sites of the MPS `state`.
 """
 eachsite(ψ::AbstractMPS) = eachindex(ψ)
+
+TensorKit.leftunit(ψ::AbstractMPS) = only(sectors(leftunitspace(left_virtualspace(ψ, 1))))
+TensorKit.rightunit(ψ::AbstractMPS) = only(sectors(rightunitspace(right_virtualspace(ψ, 1))))
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index b9e1e8734..b96886cef 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -27,10 +27,10 @@ By convention, we have that:
 ## Constructors
     FiniteMPS([f, eltype], physicalspaces::Vector{<:Union{S,CompositeSpace{S}}},
               maxvirtualspaces::Union{S,Vector{S}};
-              normalize=true, left=oneunit(S), right=oneunit(S)) where {S<:ElementarySpace}
+              normalize=true, left=unitspace(S), right=unitspace(S)) where {S<:ElementarySpace}
     FiniteMPS([f, eltype], N::Int, physicalspace::Union{S,CompositeSpace{S}},
               maxvirtualspaces::Union{S,Vector{S}};
-              normalize=true, left=oneunit(S), right=oneunit(S)) where {S<:ElementarySpace}
+              normalize=true, left=unitspace(S), right=unitspace(S)) where {S<:ElementarySpace}
     FiniteMPS(As::Vector{<:GenericMPSTensor}; normalize=false, overwrite=false)
 
 Construct an MPS via a specification of physical and virtual spaces, or from a list of
@@ -54,8 +54,8 @@ total charge can be constructed by passing a non-trivially charged vector space
 ### Keywords
 - `normalize=true`: normalize the constructed state
 - `overwrite=false`: overwrite the given input tensors
-- `left=oneunit(S)`: left-most virtual space
-- `right=oneunit(S)`: right-most virtual space
+- `left=unitspace(S)`: left-most virtual space
+- `right=unitspace(S)`: right-most virtual space
 """
 struct FiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractFiniteMPS
     ALs::Vector{Union{Missing, A}}
@@ -235,7 +235,7 @@ end
 
 function FiniteMPS(
         f, elt, Pspaces::Vector{<:Union{S, CompositeSpace{S}}}, maxVspaces::Vector{S};
-        normalize = true, left::S = oneunit(S), right::S = oneunit(S)
+        normalize = true, left::S = unitspace(S), right::S = unitspace(S)
     ) where {S <: ElementarySpace}
     N = length(Pspaces)
     length(maxVspaces) == N - 1 ||
@@ -292,7 +292,7 @@ end
 # construct from dense state
 # TODO: make planar?
 function FiniteMPS(ψ::AbstractTensor)
-    U = ones(scalartype(ψ), oneunit(spacetype(ψ)))
+    U = ones(scalartype(ψ), unitspace(spacetype(ψ)))
     A = _transpose_front(
         U * transpose(ψ * U', ((), reverse(ntuple(identity, numind(ψ) + 1))))
     )
@@ -366,10 +366,10 @@ function Base.convert(::Type{TensorMap}, ψ::FiniteMPS)
     end
 
     # remove utility legs
-    space(T, 1) == oneunit(spacetype(T)) || throw(ArgumentError("utility leg not trivial"))
-    space(T, numind(T)) == oneunit(spacetype(T))' ||
+    space(T, 1) == unitspace(spacetype(T)) || throw(ArgumentError("utility leg not trivial"))
+    space(T, numind(T)) == unitspace(spacetype(T))' ||
         throw(ArgumentError("utility leg not trivial"))
-    U = ones(scalartype(ψ), oneunit(spacetype(ψ)))
+    U = ones(scalartype(ψ), unitspace(spacetype(ψ)))
     UTU = transpose(
         U' * _transpose_tail(T * U), (reverse(ntuple(identity, numind(T) - 2)), ())
     )
@@ -410,12 +410,12 @@ end
 
 """
     max_virtualspaces(ψ::FiniteMPS)
-    max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=oneunit(S), right=oneunit(S))
+    max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=unitspace(S), right=unitspace(S))
 
 Compute the maximal virtual spaces of a given finite MPS or its physical spaces.
 """
 function max_virtualspaces(
-        Ps::Vector{<:Union{S, CompositeSpace{S}}}; left = oneunit(S), right = oneunit(S)
+        Ps::Vector{<:Union{S, CompositeSpace{S}}}; left = unitspace(S), right = unitspace(S)
     ) where {S <: ElementarySpace}
     Vs = similar(Ps, length(Ps) + 1)
     Vs[1] = left
diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl
index 6fd61f30e..2b10cd762 100644
--- a/src/states/quasiparticle_state.jl
+++ b/src/states/quasiparticle_state.jl
@@ -35,7 +35,7 @@ end
 #constructors
 function LeftGaugedQP(
         datfun, left_gs, right_gs = left_gs;
-        sector = one(sectortype(left_gs)), momentum = 0.0
+        sector = leftunit(left_gs), momentum = 0.0
     )
     # find the left null spaces for the TNS
     excitation_space = Vect[typeof(sector)](sector => 1)
@@ -56,7 +56,7 @@ function LeftGaugedQP(
 end
 function LeftGaugedQP(
         datfun, left_gs::MultilineMPS, right_gs::MultilineMPS = left_gs;
-        sector = one(sectortype(left_gs)), momentum = 0.0
+        sector = leftunit(left_gs), momentum = 0.0
     )
     # not sure why this is needed for type stability
     Tresult = leftgaugedqptype(eltype(parent(left_gs)), typeof(momentum))
@@ -69,7 +69,7 @@ end
 
 function RightGaugedQP(
         datfun, left_gs, right_gs = left_gs;
-        sector = one(sectortype(left_gs)), momentum = 0.0
+        sector = leftunit(left_gs), momentum = 0.0
     )
     # find the left null spaces for the TNS
     excitation_space = Vect[typeof(sector)](sector => 1)
@@ -228,7 +228,7 @@ auxiliarysector(state::QP) = only(sectors(auxiliaryspace(state)))
 eachsite(state::QP) = eachsite(state.left_gs)
 
 istopological(qp::QP) = qp.left_gs !== qp.right_gs
-istrivial(qp::QP) = !istopological(qp) && isone(auxiliarysector(qp))
+istrivial(qp::QP) = !istopological(qp) && isunit(auxiliarysector(qp))
 
 Base.copy(a::QP) = copy!(similar(a), a)
 Base.copyto!(a::QP, b::QP) = copy!(a, b)
@@ -308,7 +308,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S <: FiniteMPS}
     elt = scalartype(v)
 
     utl = auxiliaryspace(v)
-    ou = oneunit(utl)
+    ou = leftunitspace(utl)
     utsp = ou ⊕ ou
     upper = isometry(storagetype(site_type(v.left_gs)), utsp, ou)
     lower = left_null(upper)
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index 3e2ce109b..c33f5826b 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -119,7 +119,7 @@ function transferplot end
         sector_formatter = string
     )
     if sectors === nothing
-        sectors = [one(sectortype(h.args[1]))]
+        sectors = [unit(sectortype(h.args[1]))]
     end
 
     for sector in sectors
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index 49cb96648..6dabb136c 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -56,7 +56,20 @@ Add trivial one-dimensional utility spaces with trivial sector to the left and r
 given tensor map, i.e. as the first space of the codomain and the last space of the domain.
 """
 function add_util_leg(tensor::AbstractTensorMap{T, S, N1, N2}) where {T, S, N1, N2}
-    ou = oneunit(_firstspace(tensor))
+    ou = unitspace(_firstspace(tensor))
+
+    util_front = isomorphism(storagetype(tensor), ou * codomain(tensor), codomain(tensor))
+    util_back = isomorphism(storagetype(tensor), domain(tensor), domain(tensor) * ou)
+
+    return util_front * tensor * util_back
+end
+
+function add_util_mpoleg(tensor::AbstractTensorMap{T, S, N1, N2}) where {T, S, N1, N2}
+    # separate function for mpo because add_util_leg is also used for mps from tensors
+    # and the additional legs there can be different depending on the fusion tree
+
+    #TODO?: might need a left and right utility space if MPO has module virtual space (doesn't happen for hamiltonians)
+    ou = S(rightunit(first(blocksectors(tensor))) => 1)
 
     util_front = isomorphism(storagetype(tensor), ou * codomain(tensor), codomain(tensor))
     util_back = isomorphism(storagetype(tensor), domain(tensor), domain(tensor) * ou)
diff --git a/test/algorithms.jl b/test/algorithms.jl
index 207857951..3e488a886 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -727,7 +727,7 @@ module TestAlgorithms
 
         @test expectation_value(ψ, H) ≈
             expectation_value(ψ, 1 => -g * S_x()) + expectation_value(ψ, (1, 2) => -S_zz())
-        Z_mpo = MPSKit.add_util_leg(S_z())
+        Z_mpo = MPSKit.add_util_mpoleg(S_z())
         G = correlator(ψ, Z_mpo, Z_mpo, 1, 2:5)
         G2 = correlator(ψ, S_zz(), 1, 3:2:5)
         @test isapprox(G[2], G2[1], atol = 1.0e-2)
@@ -1060,7 +1060,7 @@ module TestAlgorithms
             FiniteMPS(physicalspace(H), maxVspaces[2:(end - 1)]), H; verbosity = 0
         )
         E₀ = expectation_value(gs, H)
-        @test E₀ ≈ first(vals_dense[one(U1Irrep)])
+        @test E₀ ≈ first(vals_dense[unit(U1Irrep)])
 
         for (sector, vals) in pairs(vals_dense)
             # ED tests
@@ -1077,7 +1077,7 @@ module TestAlgorithms
             Es, Bs = excitations(H, QuasiparticleAnsatz(; tol), gs; sector, num = 1)
             Es = Es .+ E₀
             # first excited state is second eigenvalue if sector is trivial
-            @test Es[1] ≈ vals[isone(sector) ? 2 : 1] atol = 1.0e-8
+            @test Es[1] ≈ vals[isunit(sector) ? 2 : 1] atol = 1.0e-8
         end
 
         # shifted charges tests
diff --git a/test/multifusion.jl b/test/multifusion.jl
new file mode 100644
index 000000000..79662e5ee
--- /dev/null
+++ b/test/multifusion.jl
@@ -0,0 +1,169 @@
+println("
+-----------------
+|   Multifusion |
+-----------------
+")
+
+module TestMultifusion
+
+    using ..TestSetup
+    using Test, TestExtras
+    using MPSKit
+    using TensorKit
+
+    I = IsingBimodule
+
+    M = I(1, 2, 0) # σ
+    Mop = I(2, 1, 0)
+    C0 = I(1, 1, 0) # unit of C
+    C1 = I(1, 1, 1)
+    D0 = I(2, 2, 0) # unit of D
+    D1 = I(2, 2, 1)
+    V = Vect[I](M => 1)
+    Vop = Vect[I](Mop => 1)
+    PD = Vect[I](D0 => 1, D1 => 1)
+    PC = Vect[I](C0 => 1, C1 => 1)
+
+    bad_fusions = [(PC, PD), (PD, PC), (V, V), (Vop, Vop), (V, PC), (Vop, PD), (V, PD), (Vop, PC), (V, Vop), (Vop, V)]
+
+    flippy(charge::IsingBimodule) = only(charge ⊗ D1)
+
+    function TFIM_multifusion(::Type{T} = ComplexF64; g = 1.0, L = Inf, twosite = false) where {T <: Number}
+        P = Vect[I](D0 => 1, D1 => 1)
+        X = zeros(T, P ← P)
+        for (s, f) in fusiontrees(X)
+            isunit(only(f.uncoupled)) ? X[s, f] .= g : X[s, f] .= -g
+        end
+        ZZ = zeros(T, P^2 ← P^2)
+        for (s, f) in fusiontrees(ZZ)
+            s.uncoupled == map(flippy, f.uncoupled) ? ZZ[s, f] .= 1 : nothing
+        end
+
+        if L == Inf
+            lattice = twosite ? PeriodicArray([P, P]) : PeriodicArray([P])
+            H₁ = InfiniteMPOHamiltonian(lattice, i => X for i in 1:length(lattice))
+            H₂ = InfiniteMPOHamiltonian(lattice, (i, i + 1) => ZZ for i in 1:length(lattice))
+        else
+            lattice = fill(P, L)
+            H₁ = FiniteMPOHamiltonian(lattice, i => X for i in 1:L)
+            H₂ = FiniteMPOHamiltonian(lattice, (i, i + 1) => ZZ for i in 1:(L - 1))
+        end
+        return H₁ + H₂
+    end
+
+    @testset "InfiniteMPS construction" begin
+        for (P, V) in bad_fusions
+            @test_throws ArgumentError InfiniteMPS([P], [V])
+        end
+    end
+
+    @testset "FiniteMPS construction" begin
+        for (P, V) in bad_fusions
+            @test_warn "no fusion channels available at site 2" FiniteMPS(rand(2:100), P, V)
+            @test_warn "no fusion channels available at site 2" FiniteMPS(rand(2:100), P, V; left = V, right = V)
+        end
+    end
+
+    @testset "Exact diagonalization" begin
+        H = TFIM_multifusion(; g = 0, L = 4)
+        E, ψ = exact_diagonalization(H; sector = D0) # test that it runs with kwarg
+        @test isapprox(E, [-3, -1, 1, 3]; atol = 1.0e-6)
+    end
+
+    @testset "Finite systems" begin
+        L = 6
+        g = 4
+        H = TFIM_multifusion(; g = g, L = L)
+        V = Vect[I](M => 8)
+        init = FiniteMPS(L, PD, V; left = Vect[I](M => 1), right = Vect[I](M => 1))
+        v₀ = variance(init, H)
+        ψ, envs, δ = find_groundstate(init, H, DMRG())
+        v = variance(ψ, H)
+        E = expectation_value(ψ, H, envs)
+
+        ψ2, envs2, δ2 = find_groundstate(init, H, DMRG2(; trscheme = trunctol(; atol = 1.0e-6)))
+        v2 = variance(ψ2, H)
+        E2 = expectation_value(ψ2, H, envs2)
+
+        @test δ ≈ 0 atol = 1.0e-3
+        @test δ2 ≈ 0 atol = 1.0e-3
+        @test v < v₀ && v2 < v₀
+
+        @test isapprox(E, E2; atol = 1.0e-6)
+
+        ED, _ = exact_diagonalization(H; sector = D0)
+        @test isapprox(E, first(ED); atol = 1.0e-6)
+
+        excE, qp = excitations(H, QuasiparticleAnsatz(), ψ2; sector = C1, num = 1)
+        @test 0 < variance(qp[1], H) < 1.0e-8
+
+        excE_DM, qp_DM = excitations(H, FiniteExcited(; gsalg = DMRG2(; trscheme = trunctol(; atol = 1.0e-6))), ψ2; num = 1)
+        @test isapprox(first(excE_DM), first(excE) + E2; atol = 1.0e-6)
+    end
+
+    @testset "Infinite systems" begin
+        # Multifusion: effectively studying the KW dual in SSB phase
+        g = 1 / 4
+        H = TFIM_multifusion(; g = g, L = Inf, twosite = true)
+        V = Vect[I](M => 48)
+        init = InfiniteMPS([PD, PD], [V, V])
+        v₀ = variance(init, H)
+        tol = 1.0e-10
+        ψ, envs, δ = find_groundstate(init, H, IDMRG(; tol = tol, maxiter = 400))
+        E = expectation_value(ψ, H, envs)
+        v = variance(ψ, H)
+
+        ψ2, envs2, δ2 = find_groundstate(init, H, IDMRG2(; tol = tol, trscheme = trunctol(; atol = 1.0e-6), maxiter = 400))
+        E2 = expectation_value(ψ2, H, envs2)
+        v2 = variance(ψ2, H)
+
+        ψ3, envs3, δ3 = find_groundstate(init, H, VUMPS(; tol = tol, maxiter = 400))
+        E3 = expectation_value(ψ3, H, envs3)
+        v3 = variance(ψ3, H)
+
+        @test isapprox(E, E2; atol = 1.0e-6)
+        @test isapprox(E, E3; atol = 1.0e-6)
+        for delta in [δ, δ2, δ3]
+            @test delta ≈ 0 atol = 1.0e-3
+        end
+        for var in [v, v2, v3]
+            @test var < v₀
+            @test var < 1.0e-8
+        end
+
+        @test first(transfer_spectrum(ψ2; sector = C0)) ≈ 1
+        @test !(abs(first(transfer_spectrum(ψ2; sector = C1))) ≈ 1) # testing injectivity
+
+        @test only(keys(entanglement_spectrum(ψ2))) == M
+
+        momentum = 0
+        excC1, qpC1 = excitations(H, QuasiparticleAnsatz(), momentum, ψ3; sector = C1)
+        @test isapprox(first(excC1), abs(2 * (g - 1)); atol = 1.0e-6) # charged excitation lower in energy
+        @test variance(qpC1[1], H) < 1.0e-8
+
+        # diagonal test (M = D): injective GS in symmetric phase
+        Hdual = TFIM_multifusion(; g = 1 / g, L = Inf, twosite = true)
+        Vdiag = Vect[I](D0 => 24, D1 => 24)
+        initdiag = InfiniteMPS([PD, PD], [Vdiag, Vdiag])
+        gsdiag, envsdiag = find_groundstate(initdiag, Hdual, VUMPS(; tol = tol, maxiter = 400))
+        Ediag = expectation_value(gsdiag, Hdual, envsdiag)
+
+        excD1, qpD1 = excitations(Hdual, QuasiparticleAnsatz(), momentum, gsdiag; sector = D1)
+        @test isapprox(first(excD1), abs(2 * (1 / g - 1)); atol = 1.0e-6) # charged excitation lower in energy
+        @test variance(qpD1[1], Hdual) < 1.0e-8
+
+        # comparison to Z2 Ising: injective in symmetric phase
+        HZ2 = repeat(transverse_field_ising(Z2Irrep; g = 1 / g, L = Inf), 2)
+        VZ2 = Z2Space(0 => 24, 1 => 24)
+        PZ2 = Z2Space(0 => 1, 1 => 1)
+        initZ2 = InfiniteMPS([PZ2, PZ2], [VZ2, VZ2])
+        gsZ2, envsZ2 = find_groundstate(initZ2, HZ2, VUMPS(; tol = tol, maxiter = 400))
+        EZ2 = expectation_value(gsZ2, HZ2, envsZ2)
+        @test isapprox(EZ2, Ediag; atol = 1.0e-6)
+
+        excZ2_1, qpZ2_1 = excitations(HZ2, QuasiparticleAnsatz(), momentum, gsZ2; sector = Z2Irrep(1))
+        @test isapprox(first(excZ2_1), first(excD1); atol = 1.0e-6)
+        @test variance(qpZ2_1[1], HZ2) < 1.0e-8
+    end
+
+end # module TestMultifusion
diff --git a/test/operators.jl b/test/operators.jl
index f5dce96b5..17961ef00 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -255,7 +255,7 @@ module TestOperators
                 FiniteMPOHamiltonian(grid, horizontal_operators) atol = 1.0e-4
 
             H5 = changebonds(H4 / 3 + 2H4 / 3, SvdCut(; trscheme = trunctol(; atol = 1.0e-12)))
-            psi = FiniteMPS(physicalspace(H5), V ⊕ oneunit(V))
+            psi = FiniteMPS(physicalspace(H5), V ⊕ rightunitspace(V))
             @test expectation_value(psi, H4) ≈ expectation_value(psi, H5)
         end
     end
@@ -408,7 +408,7 @@ module TestOperators
 
     @testset "DenseMPO" for ham in (transverse_field_ising(), heisenberg_XXX(; spin = 1))
         pspace = physicalspace(ham, 1)
-        ou = oneunit(pspace)
+        ou = rightunitspace(pspace)
 
         ψ = InfiniteMPS([pspace], [ou ⊕ pspace])
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 556b26dec..2d0e6adb4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,9 @@ include("setup.jl")
     if GROUP == "ALL" || GROUP == "ALGORITHMS"
         @time include("algorithms.jl")
     end
+    if GROUP == "ALL" || GROUP == "MULTIFUSION"
+        @time include("multifusion.jl")
+    end
     if GROUP == "ALL" || GROUP == "OTHER"
         @time include("other.jl")
     end
diff --git a/test/setup.jl b/test/setup.jl
index a6b07bbc0..c769bc430 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -287,7 +287,7 @@ function kitaev_model(; t = 1.0, mu = 1.0, Delta = 1.0, L = Inf)
     else
         lattice = fill(space(TB, 1), L)
         onsite_terms = ((i,) => CP for i in 1:L)
-        twosite_terms = ((i, i + 1) => TP + SC for i in 1:(L - 1))
+        twosite_terms = ((i, i + 1) => TB + SC for i in 1:(L - 1))
         terms = Iterators.flatten(twosite_terms, onsite_terms)
         return FiniteMPOHamiltonian(lattice, terms)
     end