diff --git a/.gitignore b/.gitignore
index b0acbe33b..8e2b6fec5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@ Manifest.toml
 .vscode
 .DS_Store
 **/dev/
+benchmark/results
diff --git a/benchmark/MPSKitBenchmarks/MPSKitBenchmarks.jl b/benchmark/MPSKitBenchmarks/MPSKitBenchmarks.jl
new file mode 100644
index 000000000..c70a01f70
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/MPSKitBenchmarks.jl
@@ -0,0 +1,59 @@
+module MPSKitBenchmarks
+
+using BenchmarkTools
+using MPSKit
+using TOML
+
+include("utils/BenchUtils.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 10000
+BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.15
+BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = 0.01
+
+const PARAMS_PATH = joinpath(@__DIR__, "etc", "params.json")
+const SUITE = BenchmarkGroup()
+const MODULES = Dict{String, Symbol}(
+    "derivatives" => :DerivativesBenchmarks
+)
+
+include("derivatives/DerivativesBenchmarks.jl")
+
+load!(id::AbstractString; kwargs...) = load!(SUITE, id; kwargs...)
+
+function load!(group::BenchmarkGroup, id::AbstractString; tune::Bool = false)
+    modsym = MODULES[id]
+    modpath = joinpath(dirname(@__FILE__), id, "$(modsym).jl")
+    Core.eval(@__MODULE__, :(include($modpath)))
+    mod = Core.eval(@__MODULE__, modsym)
+    modsuite = @invokelatest getglobal(mod, :SUITE)
+    group[id] = modsuite
+    if tune
+        results = BenchmarkTools.load(PARAMS_PATH)[1]
+        haskey(results, id) && loadparams!(modsuite, results[id], :evals)
+    end
+    return group
+end
+
+loadall!(; kwargs...) = loadall!(SUITE; kwargs...)
+
+function loadall!(group::BenchmarkGroup; verbose::Bool = true, tune::Bool = false)
+    for id in keys(MODULES)
+        if verbose
+            print("loading group $(repr(id))... ")
+            time = @elapsed load!(group, id, tune = false)
+            println("done (took $time seconds)")
+        else
+            load!(group, id; tune = false)
+        end
+    end
+    if tune
+        results = BenchmarkTools.load(PARAMS_PATH)[1]
+        for (id, suite) in group
+            haskey(results, id) && loadparams!(suite, results[id], :evals)
+        end
+    end
+    return group
+end
+
+end
diff --git a/benchmark/MPSKitBenchmarks/derivatives/AC2_benchmarks.jl b/benchmark/MPSKitBenchmarks/derivatives/AC2_benchmarks.jl
new file mode 100644
index 000000000..b7aeb8e68
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/AC2_benchmarks.jl
@@ -0,0 +1,107 @@
+struct AC2Spec{S <: ElementarySpace} # <: BenchmarkSpec
+    physicalspaces::NTuple{2, S}
+    mps_virtualspaces::NTuple{3, S}
+    mpo_virtualspaces::NTuple{3, SumSpace{S}}
+    nonzero_keys::NTuple{2, Vector{Tuple{Int, Int}}}
+end
+
+function AC2Spec(mps, mpo; site = length(mps) ÷ 2)
+    physicalspaces = (physicalspace(mps, site), physicalspace(mps, site + 1))
+    mps_virtualspaces = (left_virtualspace(mps, site), right_virtualspace(mps, site), right_virtualspace(mps, site + 1))
+    mpo_virtualspaces = (left_virtualspace(mpo, site), right_virtualspace(mpo, site), right_virtualspace(mpo, site + 1))
+    ks = (
+        map(x -> (x.I[1], x.I[4]), nonzero_keys(mpo[site])),
+        map(x -> (x.I[1], x.I[4]), nonzero_keys(mpo[site])),
+    )
+    return AC2Spec(physicalspaces, mps_virtualspaces, mpo_virtualspaces, ks)
+end
+
+benchname(spec::AC2Spec) = dim(spec.mps_virtualspaces[1]), dim(spec.mpo_virtualspaces[1])
+
+# Benchmarks
+# ----------
+function MPSKit.AC2_hamiltonian(spec::AC2Spec{S}; T::Type = Float64) where {S}
+    GL = randn(T, spec.mps_virtualspaces[1] ⊗ spec.mpo_virtualspaces[1]' ← spec.mps_virtualspaces[1])
+    GR = randn(T, spec.mps_virtualspaces[3] ⊗ spec.mpo_virtualspaces[3] ← spec.mps_virtualspaces[3])
+
+    W1 = MPSKit.JordanMPOTensor{T, S}(undef, spec.mpo_virtualspaces[1] ⊗ spec.physicalspaces[1] ← spec.physicalspaces[1] ⊗ spec.mpo_virtualspaces[2])
+    for (r, c) in spec.nonzero_keys[1]
+        r == c == 1 && continue
+        r == size(W1, 1) && c == size(W1, 4) && continue
+        W1[r, 1, 1, c] = randn!(W1[r, 1, 1, c])
+    end
+    W2 = MPSKit.JordanMPOTensor{T, S}(undef, spec.mpo_virtualspaces[2] ⊗ spec.physicalspaces[2] ← spec.physicalspaces[2] ⊗ spec.mpo_virtualspaces[3])
+    for (r, c) in spec.nonzero_keys[2]
+        r == c == 1 && continue
+        r == size(W2, 1) && c == size(W2, 4) && continue
+        W2[r, 1, 1, c] = randn!(W2[r, 1, 1, c])
+    end
+    H = InfiniteMPOHamiltonian(PeriodicVector([W1, W2]))
+
+    A = PeriodicVector(
+        [
+            randn(T, spec.mps_virtualspaces[1] ⊗ spec.physicalspaces[1] ← spec.mps_virtualspaces[2]),
+            randn(T, spec.mps_virtualspaces[2] ⊗ spec.physicalspaces[1] ← spec.mps_virtualspaces[3]),
+        ]
+    )
+    C = PeriodicVector(
+        [
+            randn(T, spec.mps_virtualspaces[2] ← spec.mps_virtualspaces[2]),
+            randn(T, spec.mps_virtualspaces[3] ← spec.mps_virtualspaces[3]),
+        ]
+    )
+    psi = InfiniteMPS{eltype(A), eltype(C)}(A, A, C, A)
+
+    GLs, GRs = MPSKit.initialize_environments(psi, H, psi)
+    envs = MPSKit.InfiniteEnvironments(GLs, GRs)
+
+    return MPSKit.AC2_hamiltonian(1, psi, H, psi, envs)
+end
+
+function contraction_benchmark(spec::AC2Spec; T::Type = Float64)
+    V = spec.mps_virtualspaces[1] ⊗ spec.physicalspaces[1] ← spec.mps_virtualspaces[3] ⊗ spec.physicalspaces[2]'
+    H_eff = MPSKit.AC2_hamiltonian(spec; T)
+    return @benchmarkable $H_eff * x setup = x = randn($T, $V)
+end
+
+function preparation_benchmark(spec::AC2Spec; T::Type = Float64)
+    V = spec.mps_virtualspaces[1] ⊗ spec.physicalspaces[1] ← spec.mps_virtualspaces[3] ⊗ spec.physicalspaces[2]'
+    H_eff = MPSKit.AC2_hamiltonian(spec; T)
+    return @benchmarkable MPSKit.prepare_operator!!($H_eff)
+end
+
+function prepared_benchmark(spec::AC2Spec; T::Type = Float64)
+    V = spec.mps_virtualspaces[1] ⊗ spec.physicalspaces[1] ← spec.mps_virtualspaces[3] ⊗ spec.physicalspaces[2]'
+    H_eff = MPSKit.AC2_hamiltonian(spec; T)
+    H_prep = MPSKit.prepare_operator!!(H_eff)
+    return @benchmarkable $H_prep * x setup = x = randn($T, $V)
+end
+
+# Converters
+# ----------
+function tomlify(spec::AC2Spec)
+    return Dict(
+        "physicalspaces" => collect(tomlify.(spec.physicalspaces)),
+        "mps_virtualspaces" => collect(tomlify.(spec.mps_virtualspaces)),
+        "mpo_virtualspaces" => collect(tomlify.(spec.mpo_virtualspaces)),
+        "nonzero_keys" => collect(map(Base.Fix1(map, collect), spec.nonzero_keys))
+    )
+end
+
+function untomlify(::Type{AC2Spec}, x)
+    to_space = Base.Fix1(untomlify, VectorSpace)
+    physicalspaces = Tuple(map(to_space, x["physicalspaces"]))
+    mps_virtualspaces = Tuple(map(to_space, x["mps_virtualspaces"]))
+    mpo_virtualspaces = Tuple(map(to_space, x["mpo_virtualspaces"]))
+    nonzero_keys = Tuple(map(Base.Fix1(map, Tuple), x["nonzero_keys"]))
+    return AC2Spec(physicalspaces, mps_virtualspaces, mpo_virtualspaces, nonzero_keys)
+end
+
+function Base.convert(::Type{AC2Spec{S₁}}, spec::AC2Spec{S₂}) where {S₁, S₂}
+    return S₁ === S₂ ? spec : AC2Spec(
+            S₁.(spec.physicalspaces),
+            S₁.(spec.mps_virtualspaces),
+            SumSpace{S₁}.(spec.mpo_virtualspaces),
+            spec.nonzero_keys
+        )
+end
diff --git a/benchmark/MPSKitBenchmarks/derivatives/DerivativesBenchmarks.jl b/benchmark/MPSKitBenchmarks/derivatives/DerivativesBenchmarks.jl
new file mode 100644
index 000000000..0b423dcf2
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/DerivativesBenchmarks.jl
@@ -0,0 +1,48 @@
+module DerivativesBenchmarks
+
+export AC2Spec
+
+using BenchmarkTools
+using TOML
+using TensorKit
+using BlockTensorKit
+using MPSKit
+using ..BenchUtils
+import ..BenchUtils: tomlify, untomlify
+
+const SUITE = BenchmarkGroup()
+
+const allparams = Dict(
+    "heisenberg_nn" => TOML.parsefile(joinpath(@__DIR__, "heisenberg_nn.toml")),
+    "heisenberg_nnn" => TOML.parsefile(joinpath(@__DIR__, "heisenberg_nnn.toml")),
+    "heisenberg_cylinder" => TOML.parsefile(joinpath(@__DIR__, "heisenberg_cylinder.toml")),
+    "heisenberg_coulomb" => TOML.parsefile(joinpath(@__DIR__, "heisenberg_coulomb.toml"))
+)
+
+include("AC2_benchmarks.jl")
+
+T = Float64
+
+suite_init = addgroup!(SUITE, "AC2_preparation")
+suite_apply = addgroup!(SUITE, "AC2_contraction")
+suite_prepped = addgroup!(SUITE, "AC2_prepared")
+
+for (model, params) in allparams
+    g_prep = addgroup!(suite_init, model)
+    g_prepped = addgroup!(suite_prepped, model)
+    g_contract = addgroup!(suite_apply, model)
+    for (symmetry, specs) in params
+        g_prep_sym = addgroup!(g_prep, symmetry)
+        g_contract_sym = addgroup!(g_contract, symmetry)
+        g_prepped_sym = addgroup!(g_prepped, symmetry)
+        for spec_dict in specs
+            spec = untomlify(AC2Spec, spec_dict)
+            name = benchname(spec)
+            g_contract_sym[name] = contraction_benchmark(spec; T)
+            g_prep_sym[name] = preparation_benchmark(spec; T)
+            g_prepped_sym[name] = prepared_benchmark(spec; T)
+        end
+    end
+end
+
+end
diff --git a/benchmark/MPSKitBenchmarks/derivatives/heisenberg_coulomb.toml b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_coulomb.toml
new file mode 100644
index 000000000..0a945aced
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_coulomb.toml
@@ -0,0 +1,132 @@
+[[Trivial]]
+mps_virtualspaces = ["ℂ^1", "ℂ^3", "ℂ^1"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^4", "ℂ^4", "ℂ^4"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^7", "ℂ^7", "ℂ^7"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^19", "ℂ^19", "ℂ^19"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^40", "ℂ^40", "ℂ^40"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^76", "ℂ^76", "ℂ^76"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^140", "ℂ^143", "ℂ^140"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^272", "ℂ^274", "ℂ^272"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^498 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 2, 1 => 1, -1 => 1)", "Rep[U₁](0 => 2, 1 => 1, -1 => 1)", "Rep[U₁](0 => 2, 1 => 1, -1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 7, 1 => 5, -1 => 5, 2 => 1, -2 => 1)", "Rep[U₁](0 => 7, 1 => 5, -1 => 5, 2 => 1, -2 => 1)", "Rep[U₁](0 => 7, 1 => 5, -1 => 5, 2 => 1, -2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 14, 1 => 10, -1 => 10, 2 => 3, -2 => 3)", "Rep[U₁](0 => 14, 1 => 10, -1 => 10, 2 => 3, -2 => 3)", "Rep[U₁](0 => 14, 1 => 10, -1 => 10, 2 => 3, -2 => 3)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 24, 1 => 18, -1 => 18, 2 => 7, -2 => 7, 3 => 1, -3 => 1)", "Rep[U₁](0 => 24, 1 => 18, -1 => 18, 2 => 7, -2 => 7, 3 => 1, -3 => 1)", "Rep[U₁](0 => 24, 1 => 18, -1 => 18, 2 => 7, -2 => 7, 3 => 1, -3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 42, 1 => 32, -1 => 32, 2 => 14, -2 => 14, 3 => 3, -3 => 3)", "Rep[U₁](0 => 43, 1 => 33, -1 => 33, 2 => 14, -2 => 14, 3 => 3, -3 => 3)", "Rep[U₁](0 => 42, 1 => 32, -1 => 32, 2 => 14, -2 => 14, 3 => 3, -3 => 3)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 80, 1 => 62, -1 => 62, 2 => 28, -2 => 28, 3 => 6, -3 => 6)", "Rep[U₁](0 => 80, 1 => 63, -1 => 63, 2 => 28, -2 => 28, 3 => 6, -3 => 6)", "Rep[U₁](0 => 80, 1 => 62, -1 => 62, 2 => 28, -2 => 28, 3 => 6, -3 => 6)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>166, 1=>166, -1=>166) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1)", "Rep[SU₂](1 => 1)", "Rep[SU₂](0 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1, 1 => 1)", "Rep[SU₂](0 => 1, 1 => 1)", "Rep[SU₂](0 => 1, 1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 2, 1 => 4, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 4, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 4, 2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 4, 1 => 7, 2 => 3)", "Rep[SU₂](0 => 4, 1 => 7, 2 => 3)", "Rep[SU₂](0 => 4, 1 => 7, 2 => 3)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 6, 1 => 11, 2 => 6, 3 => 1)", "Rep[SU₂](0 => 6, 1 => 11, 2 => 6, 3 => 1)", "Rep[SU₂](0 => 6, 1 => 11, 2 => 6, 3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 10, 1 => 18, 2 => 11, 3 => 3)", "Rep[SU₂](0 => 10, 1 => 19, 2 => 11, 3 => 3)", "Rep[SU₂](0 => 10, 1 => 18, 2 => 11, 3 => 3)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 18, 1 => 34, 2 => 22, 3 => 6)", "Rep[SU₂](0 => 17, 1 => 35, 2 => 22, 3 => 6)", "Rep[SU₂](0 => 18, 1 => 34, 2 => 22, 3 => 6)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 28, 1 => 58, 2 => 43, 3 => 14, 4 => 1)", "Rep[SU₂](0 => 30, 1 => 59, 2 => 43, 3 => 13, 4 => 1)", "Rep[SU₂](0 => 28, 1 => 58, 2 => 43, 3 => 14, 4 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 41, 1 => 84, 2 => 64, 3 => 26, 4 => 11, 5 => 6, 6 => 3)", "Rep[SU₂](0 => 34, 1 => 71, 2 => 63, 3 => 30, 4 => 12, 5 => 6, 6 => 3, 7 => 1)", "Rep[SU₂](0 => 40, 1 => 83, 2 => 62, 3 => 29, 4 => 12, 5 => 7, 6 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]], [[1, 1], [3, 3], [2, 2], [2, 3], [1, 2], [1, 3]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>166) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
diff --git a/benchmark/MPSKitBenchmarks/derivatives/heisenberg_cylinder.toml b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_cylinder.toml
new file mode 100644
index 000000000..bd1e55b79
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_cylinder.toml
@@ -0,0 +1,132 @@
+[[Trivial]]
+mps_virtualspaces = ["ℂ^3", "ℂ^1", "ℂ^3"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^3", "ℂ^4", "ℂ^3"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^16", "ℂ^12", "ℂ^16"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^43", "ℂ^41", "ℂ^45"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^110", "ℂ^104", "ℂ^110"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^269", "ℂ^267", "ℂ^266"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^668", "ℂ^669", "ℂ^669"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1 ⊞ ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 2, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 4, 1 => 3, -1 => 3, 2 => 2, -2 => 2, 3 => 1, -3 => 1)", "Rep[U₁](0 => 4, 1 => 3, -1 => 3, 2 => 1, -2 => 1)", "Rep[U₁](0 => 4, 1 => 3, -1 => 3, 2 => 2, -2 => 2, 3 => 1, -3 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 13, 1 => 9, -1 => 9, 2 => 5, -2 => 5, 3 => 1, -3 => 1)", "Rep[U₁](0 => 11, 1 => 10, -1 => 10, 2 => 4, -2 => 4, 3 => 1, -3 => 1)", "Rep[U₁](0 => 13, 1 => 9, -1 => 9, 2 => 6, -2 => 6, 3 => 1, -3 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 26, 1 => 22, -1 => 22, 2 => 13, -2 => 13, 3 => 6, -3 => 6, 4 => 1, -4 => 1)", "Rep[U₁](0 => 26, 1 => 21, -1 => 21, 2 => 12, -2 => 12, 3 => 5, -3 => 5, 4 => 1, -4 => 1)", "Rep[U₁](0 => 26, 1 => 22, -1 => 22, 2 => 13, -2 => 13, 3 => 6, -3 => 6, 4 => 1, -4 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 59, 1 => 53, -1 => 53, 2 => 32, -2 => 32, 3 => 14, -3 => 14, 4 => 5, -4 => 5, 5 => 1, -5 => 1)", "Rep[U₁](0 => 61, 1 => 51, -1 => 51, 2 => 35, -2 => 35, 3 => 12, -3 => 12, 4 => 4, -4 => 4, 5 => 1, -5 => 1)", "Rep[U₁](0 => 58, 1 => 52, -1 => 52, 2 => 32, -2 => 32, 3 => 14, -3 => 14, 4 => 5, -4 => 5, 5 => 1, -5 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 140, 1 => 121, -1 => 121, 2 => 83, -2 => 83, 3 => 41, -3 => 41, 4 => 14, -4 => 14, 5 => 4, -5 => 4, 6 => 1, -6 => 1)", "Rep[U₁](0 => 145, 1 => 131, -1 => 131, 2 => 83, -2 => 83, 3 => 38, -3 => 38, 4 => 9, -4 => 9, 5 => 1, -5 => 1)", "Rep[U₁](0 => 141, 1 => 122, -1 => 122, 2 => 83, -2 => 83, 3 => 40, -3 => 40, 4 => 14, -4 => 14, 5 => 4, -5 => 4, 6 => 1, -6 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 328, 1 => 293, -1 => 293, 2 => 204, -2 => 204, 3 => 107, -3 => 107, 4 => 43, -4 => 43, 5 => 11, -5 => 11, 6 => 1, -6 => 1)", "Rep[U₁](0 => 344, 1 => 298, -1 => 298, 2 => 206, -2 => 206, 3 => 102, -3 => 102, 4 => 37, -4 => 37, 5 => 8, -5 => 8, 6 => 1, -6 => 1)", "Rep[U₁](0 => 329, 1 => 293, -1 => 293, 2 => 203, -2 => 203, 3 => 108, -3 => 108, 4 => 43, -4 => 43, 5 => 11, -5 => 11, 6 => 1, -6 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 771, 1 => 692, -1 => 692, 2 => 495, -2 => 495, 3 => 283, -3 => 283, 4 => 124, -4 => 124, 5 => 41, -5 => 41, 6 => 8, -6 => 8, 7 => 1, -7 => 1)", "Rep[U₁](0 => 787, 1 => 704, -1 => 704, 2 => 493, -2 => 493, 3 => 275, -3 => 275, 4 => 118, -4 => 118, 5 => 36, -5 => 36, 6 => 8, -6 => 8, 7 => 1, -7 => 1)", "Rep[U₁](0 => 771, 1 => 692, -1 => 692, 2 => 495, -2 => 495, 3 => 282, -3 => 282, 4 => 124, -4 => 124, 5 => 41, -5 => 41, 6 => 8, -6 => 8, 7 => 1, -7 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](0 => 1)", "Rep[SU₂](1 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](0 => 1, 1 => 1)", "Rep[SU₂](1 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1, 1 => 1, 2 => 1, 3 => 1)", "Rep[SU₂](0 => 1, 1 => 2, 2 => 1)", "Rep[SU₂](0 => 1, 1 => 1, 2 => 1, 3 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 4, 1 => 4, 2 => 4, 3 => 1)", "Rep[SU₂](0 => 1, 1 => 6, 2 => 3, 3 => 1)", "Rep[SU₂](0 => 4, 1 => 3, 2 => 5, 3 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 4, 1 => 9, 2 => 7, 3 => 5, 4 => 1)", "Rep[SU₂](0 => 5, 1 => 9, 2 => 7, 3 => 4, 4 => 1)", "Rep[SU₂](0 => 4, 1 => 9, 2 => 7, 3 => 5, 4 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 6, 1 => 21, 2 => 18, 3 => 9, 4 => 4, 5 => 1)", "Rep[SU₂](0 => 10, 1 => 16, 2 => 23, 3 => 8, 4 => 3, 5 => 1)", "Rep[SU₂](0 => 6, 1 => 20, 2 => 18, 3 => 9, 4 => 4, 5 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 19, 1 => 38, 2 => 42, 3 => 27, 4 => 10, 5 => 3, 6 => 1)", "Rep[SU₂](0 => 14, 1 => 48, 2 => 45, 3 => 29, 4 => 8, 5 => 1)", "Rep[SU₂](0 => 19, 1 => 39, 2 => 43, 3 => 26, 4 => 10, 5 => 3, 6 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 35, 1 => 89, 2 => 97, 3 => 64, 4 => 32, 5 => 10, 6 => 1)", "Rep[SU₂](0 => 46, 1 => 92, 2 => 104, 3 => 65, 4 => 29, 5 => 7, 6 => 1)", "Rep[SU₂](0 => 36, 1 => 90, 2 => 95, 3 => 65, 4 => 32, 5 => 10, 6 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 79, 1 => 197, 2 => 212, 3 => 159, 4 => 83, 5 => 33, 6 => 7, 7 => 1)", "Rep[SU₂](0 => 83, 1 => 211, 2 => 218, 3 => 157, 4 => 82, 5 => 28, 6 => 7, 7 => 1)", "Rep[SU₂](0 => 79, 1 => 197, 2 => 213, 3 => 158, 4 => 83, 5 => 33, 6 => 7, 7 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 169, 1 => 427, 2 => 488, 3 => 380, 4 => 217, 5 => 92, 6 => 29, 7 => 6, 8 => 1)", "Rep[SU₂](0 => 185, 1 => 432, 2 => 510, 3 => 380, 4 => 210, 5 => 89, 6 => 26, 7 => 6)", "Rep[SU₂](0 => 168, 1 => 426, 2 => 489, 3 => 380, 4 => 217, 5 => 92, 6 => 29, 7 => 6, 8 => 1)"]
+nonzero_keys = [[[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]], [[1, 1], [10, 10], [7, 7], [6, 6], [9, 9], [8, 8], [5, 5], [2, 10], [3, 10], [4, 10], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
diff --git a/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nn.toml b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nn.toml
new file mode 100644
index 000000000..1e508e53a
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nn.toml
@@ -0,0 +1,132 @@
+[[Trivial]]
+mps_virtualspaces = ["ℂ^3", "ℂ^1", "ℂ^3"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^7", "ℂ^7", "ℂ^7"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^16", "ℂ^16", "ℂ^16"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^44", "ℂ^44", "ℂ^44"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^108", "ℂ^108", "ℂ^108"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^271", "ℂ^271", "ℂ^271"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^661", "ℂ^667", "ℂ^669"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)", "Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)", "Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)", "Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)", "Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)", "Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)", "Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)", "Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)", "Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 153, 1 => 130, -1 => 130, 2 => 79, -2 => 79, 3 => 33, -3 => 33, 4 => 10, -4 => 10, 5 => 2, -5 => 2)", "Rep[U₁](0 => 157, 1 => 133, -1 => 133, 2 => 80, -2 => 80, 3 => 33, -3 => 33, 4 => 8, -4 => 8, 5 => 1, -5 => 1)", "Rep[U₁](0 => 159, 1 => 134, -1 => 134, 2 => 80, -2 => 80, 3 => 32, -3 => 32, 4 => 8, -4 => 8, 5 => 1, -5 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 307, 1 => 273, -1 => 273, 2 => 194, -2 => 194, 3 => 114, -3 => 114, 4 => 58, -4 => 58, 5 => 24, -5 => 24, 6 => 6, -6 => 6, 7 => 1, -7 => 1)", "Rep[U₁](0 => 317, 1 => 280, -1 => 280, 2 => 194, -2 => 194, 3 => 110, -3 => 110, 4 => 53, -4 => 53, 5 => 20, -5 => 20, 6 => 5, -6 => 5, 7 => 1, -7 => 1)", "Rep[U₁](0 => 313, 1 => 278, -1 => 278, 2 => 196, -2 => 196, 3 => 114, -3 => 114, 4 => 56, -4 => 56, 5 => 20, -5 => 20, 6 => 4, -6 => 4)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 987, 1 => 830, -1 => 830, 2 => 432, -2 => 432, 3 => 161, -3 => 161, 4 => 77, -4 => 77, 5 => 28, -5 => 28, 6 => 7, -6 => 7, 7 => 1, -7 => 1)", "Rep[U₁](0 => 701, 1 => 635, -1 => 635, 2 => 478, -2 => 478, 3 => 303, -3 => 303, 4 => 162, -4 => 162, 5 => 70, -5 => 70, 6 => 22, -6 => 22, 7 => 4, -7 => 4)", "Rep[U₁](0 => 987, 1 => 830, -1 => 830, 2 => 431, -2 => 431, 3 => 161, -3 => 161, 4 => 77, -4 => 77, 5 => 28, -5 => 28, 6 => 7, -6 => 7, 7 => 1, -7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](0 => 1)", "Rep[SU₂](1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 2, 1 => 3, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 3, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 3, 2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)", "Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)", "Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)", "Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)", "Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)", "Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)", "Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 23, 1 => 51, 2 => 46, 3 => 23, 4 => 8, 5 => 2)", "Rep[SU₂](0 => 24, 1 => 53, 2 => 47, 3 => 25, 4 => 7, 5 => 1)", "Rep[SU₂](0 => 25, 1 => 54, 2 => 48, 3 => 24, 4 => 7, 5 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 34, 1 => 79, 2 => 80, 3 => 56, 4 => 34, 5 => 18, 6 => 5, 7 => 1)", "Rep[SU₂](0 => 37, 1 => 86, 2 => 84, 3 => 57, 4 => 33, 5 => 15, 6 => 4, 7 => 1)", "Rep[SU₂](0 => 35, 1 => 82, 2 => 82, 3 => 58, 4 => 36, 5 => 16, 6 => 4)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 157, 1 => 398, 2 => 271, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)", "Rep[SU₂](0 => 66, 1 => 157, 2 => 175, 3 => 141, 4 => 92, 5 => 48, 6 => 18, 7 => 4)", "Rep[SU₂](0 => 157, 1 => 399, 2 => 270, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1196, 1 => 1678, 2 => 483, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)", "Rep[SU₂](0 => 603, 1 => 1213, 2 => 280, 3 => 237, 4 => 153, 5 => 76, 6 => 28, 7 => 7, 8 => 1)", "Rep[SU₂](0 => 1198, 1 => 1646, 2 => 502, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
diff --git a/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nnn.toml b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nnn.toml
new file mode 100644
index 000000000..1e508e53a
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/derivatives/heisenberg_nnn.toml
@@ -0,0 +1,132 @@
+[[Trivial]]
+mps_virtualspaces = ["ℂ^3", "ℂ^1", "ℂ^3"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^7", "ℂ^7", "ℂ^7"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^16", "ℂ^16", "ℂ^16"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^44", "ℂ^44", "ℂ^44"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^108", "ℂ^108", "ℂ^108"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^271", "ℂ^271", "ℂ^271"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+[[Trivial]]
+mps_virtualspaces = ["ℂ^661", "ℂ^667", "ℂ^669"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)", "(ℂ^1 ⊞ ℂ^3 ⊞ ℂ^1)"]
+physicalspaces = ["ℂ^3", "ℂ^3"]
+
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)", "Rep[U₁](0 => 3, 1 => 2, -1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)", "Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)", "Rep[U₁](0 => 6, 1 => 4, -1 => 4, 2 => 1, -2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)", "Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)", "Rep[U₁](0 => 12, 1 => 10, -1 => 10, 2 => 5, -2 => 5, 3 => 1, -3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)", "Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)", "Rep[U₁](0 => 32, 1 => 25, -1 => 25, 2 => 11, -2 => 11, 3 => 2, -3 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)", "Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)", "Rep[U₁](0 => 71, 1 => 58, -1 => 58, 2 => 31, -2 => 31, 3 => 10, -3 => 10, 4 => 1, -4 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 153, 1 => 130, -1 => 130, 2 => 79, -2 => 79, 3 => 33, -3 => 33, 4 => 10, -4 => 10, 5 => 2, -5 => 2)", "Rep[U₁](0 => 157, 1 => 133, -1 => 133, 2 => 80, -2 => 80, 3 => 33, -3 => 33, 4 => 8, -4 => 8, 5 => 1, -5 => 1)", "Rep[U₁](0 => 159, 1 => 134, -1 => 134, 2 => 80, -2 => 80, 3 => 32, -3 => 32, 4 => 8, -4 => 8, 5 => 1, -5 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 307, 1 => 273, -1 => 273, 2 => 194, -2 => 194, 3 => 114, -3 => 114, 4 => 58, -4 => 58, 5 => 24, -5 => 24, 6 => 6, -6 => 6, 7 => 1, -7 => 1)", "Rep[U₁](0 => 317, 1 => 280, -1 => 280, 2 => 194, -2 => 194, 3 => 110, -3 => 110, 4 => 53, -4 => 53, 5 => 20, -5 => 20, 6 => 5, -6 => 5, 7 => 1, -7 => 1)", "Rep[U₁](0 => 313, 1 => 278, -1 => 278, 2 => 196, -2 => 196, 3 => 114, -3 => 114, 4 => 56, -4 => 56, 5 => 20, -5 => 20, 6 => 4, -6 => 4)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+[["Irrep[U₁]"]]
+mps_virtualspaces = ["Rep[U₁](0 => 987, 1 => 830, -1 => 830, 2 => 432, -2 => 432, 3 => 161, -3 => 161, 4 => 77, -4 => 77, 5 => 28, -5 => 28, 6 => 7, -6 => 7, 7 => 1, -7 => 1)", "Rep[U₁](0 => 701, 1 => 635, -1 => 635, 2 => 478, -2 => 478, 3 => 303, -3 => 303, 4 => 162, -4 => 162, 5 => 70, -5 => 70, 6 => 22, -6 => 22, 7 => 4, -7 => 4)", "Rep[U₁](0 => 987, 1 => 830, -1 => 830, 2 => 431, -2 => 431, 3 => 161, -3 => 161, 4 => 77, -4 => 77, 5 => 28, -5 => 28, 6 => 7, -6 => 7, 7 => 1, -7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))", "(Rep[U₁](0=>1) ⊞ Rep[U₁](0=>1, 1=>1, -1=>1) ⊞ Rep[U₁](0=>1))"]
+physicalspaces = ["Rep[U₁](0 => 1, 1 => 1, -1 => 1)", "Rep[U₁](0 => 1, 1 => 1, -1 => 1)"]
+
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](0 => 1)", "Rep[SU₂](1 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)", "Rep[SU₂](0 => 1, 1 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 2, 1 => 3, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 3, 2 => 1)", "Rep[SU₂](0 => 2, 1 => 3, 2 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)", "Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)", "Rep[SU₂](0 => 2, 1 => 5, 2 => 4, 3 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)", "Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)", "Rep[SU₂](0 => 7, 1 => 14, 2 => 9, 3 => 2)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)", "Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)", "Rep[SU₂](0 => 13, 1 => 27, 2 => 21, 3 => 9, 4 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 23, 1 => 51, 2 => 46, 3 => 23, 4 => 8, 5 => 2)", "Rep[SU₂](0 => 24, 1 => 53, 2 => 47, 3 => 25, 4 => 7, 5 => 1)", "Rep[SU₂](0 => 25, 1 => 54, 2 => 48, 3 => 24, 4 => 7, 5 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 34, 1 => 79, 2 => 80, 3 => 56, 4 => 34, 5 => 18, 6 => 5, 7 => 1)", "Rep[SU₂](0 => 37, 1 => 86, 2 => 84, 3 => 57, 4 => 33, 5 => 15, 6 => 4, 7 => 1)", "Rep[SU₂](0 => 35, 1 => 82, 2 => 82, 3 => 58, 4 => 36, 5 => 16, 6 => 4)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 157, 1 => 398, 2 => 271, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)", "Rep[SU₂](0 => 66, 1 => 157, 2 => 175, 3 => 141, 4 => 92, 5 => 48, 6 => 18, 7 => 4)", "Rep[SU₂](0 => 157, 1 => 399, 2 => 270, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
+[["Irrep[SU₂]"]]
+mps_virtualspaces = ["Rep[SU₂](0 => 1196, 1 => 1678, 2 => 483, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)", "Rep[SU₂](0 => 603, 1 => 1213, 2 => 280, 3 => 237, 4 => 153, 5 => 76, 6 => 28, 7 => 7, 8 => 1)", "Rep[SU₂](0 => 1198, 1 => 1646, 2 => 502, 3 => 84, 4 => 49, 5 => 21, 6 => 6, 7 => 1)"]
+nonzero_keys = [[[1, 1], [3, 3], [2, 3], [1, 2]], [[1, 1], [3, 3], [2, 3], [1, 2]]]
+mpo_virtualspaces = ["(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))", "(Rep[SU₂](0=>1) ⊞ Rep[SU₂](1=>1) ⊞ Rep[SU₂](0=>1))"]
+physicalspaces = ["Rep[SU₂](1 => 1)", "Rep[SU₂](1 => 1)"]
diff --git a/benchmark/MPSKitBenchmarks/scripts/setup_heisenberg.jl b/benchmark/MPSKitBenchmarks/scripts/setup_heisenberg.jl
new file mode 100644
index 000000000..e1587648e
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/scripts/setup_heisenberg.jl
@@ -0,0 +1,207 @@
+using TensorKit
+using TensorKit: type_repr
+using MPSKit
+using BlockTensorKit
+using BlockTensorKit: nonzero_keys
+using MPSKitModels
+using TOML
+
+include(joinpath(@__DIR__, "..", "MPSKitBenchmarks.jl"))
+using .MPSKitBenchmarks
+# MPSKitBenchmarks.load!("derivatives")
+using .MPSKitBenchmarks.DerivativesBenchmarks: AC2Spec
+using .MPSKitBenchmarks.BenchUtils: tomlify
+
+# needed to make parsing correct
+BlockTensorKit.SUMSPACE_SHOW_LIMIT[] = typemax(Int)
+
+# Utility functions
+# ----------------
+# setup "product state" -> ⨂^N |↑↓ - ↓↑⟩
+function initial_state(H)
+    @assert iseven(length(H))
+    @assert allequal(physicalspace(H))
+
+    pspace = physicalspace(H, 1)
+    A = rand(oneunit(pspace) ⊗ pspace^2 ← oneunit(pspace))
+    As = MPSKit.decompose_localmps(A)
+    return FiniteMPS(repeat(As, length(H) ÷ 2))
+end
+
+function generate_spaces(H, alg; D_min = 2, D_steps = 5)
+    # compute maximal spaces
+    psi_init = initial_state(H)
+    psi, = find_groundstate(psi_init, H, alg)
+
+    D_max = maximum(dim, left_virtualspace(psi))
+    Ds = round.(Int, logrange(D_min, D_max, D_steps))
+
+    return map(Ds) do D
+        mps = changebonds(psi, SvdCut(; trscheme = truncrank(D)))
+        return AC2Spec(mps, H)
+    end
+end
+
+# Parameters
+# ----------
+T = Float64
+spin = 1
+
+D_max = 10_000
+
+D_steps = 10
+D_min = 3
+
+
+# NN
+# --
+L = 100
+lattice = FiniteChain(L)
+
+symmetry = SU2Irrep
+alg = DMRG2(;
+    maxiter = 10, tol = 1.0e-12,
+    alg_eigsolve = (; tol = 1.0e-5, dynamic_tols = false, maxiter = 3),
+    trscheme = truncrank(D_max)
+)
+
+H = heisenberg_XXX(T, symmetry, lattice; spin);
+specs_su2 = generate_spaces(H, alg; D_min, D_steps)
+
+specs_triv = filter!(convert.(AC2Spec{ComplexSpace}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 1000
+end
+
+specs_u1 = filter!(convert.(AC2Spec{U1Space}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 5000
+end
+
+output_file = joinpath(@__DIR__, "..", "derivatives", "heisenberg_nn.toml")
+open(output_file, "w") do io
+    TOML.print(
+        io, Dict(
+            type_repr(Trivial) => tomlify.(specs_triv),
+            type_repr(U1Irrep) => tomlify.(specs_u1),
+            type_repr(SU2Irrep) => tomlify.(specs_su2)
+        )
+    )
+end
+@info("Spaces written to $output_file")
+
+# NNN
+# ---
+L = 100
+lattice = FiniteChain(L)
+
+symmetry = SU2Irrep
+alg = DMRG2(;
+    maxiter = 10, tol = 1.0e-12,
+    alg_eigsolve = (; tol = 1.0e-5, dynamic_tols = false, maxiter = 3),
+    trscheme = truncrank(D_max)
+)
+
+SS = S_exchange(T, symmetry; spin)
+H = @mpoham sum(next_nearest_neighbours(lattice)) do (i, j)
+    return SS{i, j}
+end
+H = heisenberg_XXX(T, symmetry, lattice; spin);
+specs_su2 = generate_spaces(H, alg; D_min, D_steps)
+
+specs_triv = filter!(convert.(AC2Spec{ComplexSpace}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 1000
+end
+
+specs_u1 = filter!(convert.(AC2Spec{U1Space}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 5000
+end
+
+output_file = joinpath(@__DIR__, "..", "derivatives", "heisenberg_nnn.toml")
+open(output_file, "w") do io
+    TOML.print(
+        io, Dict(
+            type_repr(Trivial) => tomlify.(specs_triv),
+            type_repr(U1Irrep) => tomlify.(specs_u1),
+            type_repr(SU2Irrep) => tomlify.(specs_su2)
+        )
+    )
+end
+@info("Spaces written to $output_file")
+
+# FiniteCylinder
+# --------------
+rows = 6
+cols = 12
+lattice = FiniteCylinder(rows, rows * cols)
+
+symmetry = SU2Irrep
+alg = DMRG2(;
+    maxiter = 10, tol = 1.0e-12,
+    alg_eigsolve = (; tol = 1.0e-5, dynamic_tols = false, maxiter = 3),
+    trscheme = truncrank(D_max)
+)
+
+H = heisenberg_XXX(T, symmetry, lattice; spin);
+specs_su2 = generate_spaces(H, alg; D_min, D_steps)
+
+specs_triv = filter!(convert.(AC2Spec{ComplexSpace}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 1000
+end
+
+specs_u1 = filter!(convert.(AC2Spec{U1Space}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 5000
+end
+
+output_file = joinpath(@__DIR__, "..", "derivatives", "heisenberg_cylinder.toml")
+open(output_file, "w") do io
+    TOML.print(
+        io, Dict(
+            type_repr(Trivial) => tomlify.(specs_triv),
+            type_repr(U1Irrep) => tomlify.(specs_u1),
+            type_repr(SU2Irrep) => tomlify.(specs_su2)
+        )
+    )
+end
+@info("Spaces written to $output_file")
+
+# Coulomb
+# -------
+L = 30
+symmetry = SU2Irrep
+SS = S_exchange(T, symmetry; spin)
+lattice = fill(space(SS, 1), L)
+terms = []
+for i in 1:(L - 1), j in (i + 1):L
+    push!(terms, (i, j) => SS / abs(i - j))
+end
+H = FiniteMPOHamiltonian(lattice, terms...);
+H = changebonds(H, SvdCut(; trscheme = truncrank(500)));
+
+D_max = 1_000
+symmetry = SU2Irrep
+alg = DMRG2(;
+    maxiter = 10, tol = 1.0e-12,
+    alg_eigsolve = (; tol = 1.0e-5, dynamic_tols = false, maxiter = 3),
+    trscheme = truncrank(D_max)
+)
+
+specs_su2 = generate_spaces(H, alg; D_min, D_steps)
+
+specs_triv = filter!(convert.(AC2Spec{ComplexSpace}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 500
+end
+
+specs_u1 = filter!(convert.(AC2Spec{U1Space}, specs_su2)) do spec
+    dim(spec.mps_virtualspaces[1]) < 500
+end
+
+output_file = joinpath(@__DIR__, "..", "derivatives", "heisenberg_coulomb.toml")
+open(output_file, "w") do io
+    TOML.print(
+        io, Dict(
+            type_repr(Trivial) => tomlify.(specs_triv),
+            type_repr(U1Irrep) => tomlify.(specs_u1),
+            type_repr(SU2Irrep) => tomlify.(specs_su2)
+        )
+    )
+end
+@info("Spaces written to $output_file")
diff --git a/benchmark/MPSKitBenchmarks/utils/BenchUtils.jl b/benchmark/MPSKitBenchmarks/utils/BenchUtils.jl
new file mode 100644
index 000000000..3db7f8592
--- /dev/null
+++ b/benchmark/MPSKitBenchmarks/utils/BenchUtils.jl
@@ -0,0 +1,27 @@
+module BenchUtils
+
+export tomlify, untomlify
+
+using TensorKit
+using BlockTensorKit
+using TOML
+
+tomlify(x::VectorSpace) = sprint(show, x; context = :limited => false)
+untomlify(::Type{<:VectorSpace}, x) = eval(Meta.parse(x))
+
+
+# Type piracy but oh well
+TensorKit.ComplexSpace(V::ElementarySpace) = ComplexSpace(dim(V), isdual(V))
+
+function TensorKit.U1Space(V::SU2Space)
+    dims = TensorKit.SectorDict{U1Irrep, Int}()
+    for c in sectors(V), m in (-c.j):(c.j)
+        u1 = U1Irrep(m)
+        dims[u1] = get(dims, u1, 0) + dim(V, c)
+    end
+    return U1Space(dims; dual = isdual(V))
+end
+
+BlockTensorKit.SumSpace{S}(V::SumSpace) where {S} = SumSpace(map(S, V.spaces))
+
+end
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 000000000..c1f0ff408
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
+MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684"
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
new file mode 100644
index 000000000..1e89b2161
--- /dev/null
+++ b/benchmark/benchmarks.jl
@@ -0,0 +1,17 @@
+# Load benchmark code
+include("MPSKitBenchmarks/MPSKitBenchmarks.jl")
+const SUITE = MPSKitBenchmarks.SUITE
+
+# Populate benchmarks
+# Detect if user supplied extra arguments to load only specific modules
+# e.g. julia benchmarks.jl --modules=linalg,tensornetworks
+modules_pattern = r"(?:--modules=)(\w+)"
+arg_id = findfirst(contains(modules_pattern), ARGS)
+if isnothing(arg_id)
+    MPSKitBenchmarks.loadall!()
+else
+    modules = split(only(match(modules_pattern, ARGS[arg_id]).captures[1]), ",")
+    for m in modules
+        MPSKitBenchmarks.load!(m)
+    end
+end
diff --git a/benchmark/plot_results.jl b/benchmark/plot_results.jl
new file mode 100644
index 000000000..ce3857d62
--- /dev/null
+++ b/benchmark/plot_results.jl
@@ -0,0 +1,302 @@
+using JSON
+using DataFrames
+using CairoMakie
+using Statistics
+
+# Loading in the data
+# -------------------
+resultdir = joinpath(@__DIR__, "results")
+resultfile(i) = "results_MPSKit@bench$i.json"
+
+versions = [0, 1, 2, 3, 5]
+
+df_contract = let df = DataFrame(
+        :version => Int[], :model => String[], :symmetry => String[],
+        :D => Int[], :V => Int[], :memory => Int[], :allocs => Int[], :times => Vector{Int}[]
+    )
+
+    for version in versions
+        result = JSON.parsefile(joinpath(resultdir, resultfile(version)))
+        for (model, model_res) in result.data.derivatives.data.AC2_contraction.data
+            for (symmetry, sym_res) in model_res.data
+                for (DV, bench) in sym_res.data
+                    D, V = eval(Meta.parse(DV))::Tuple{Int, Int}
+
+                    push!(
+                        df,
+                        (version, model, symmetry, D, V, bench.memory, bench.allocs, collect(Int, bench.times))
+                    )
+                end
+            end
+        end
+    end
+    df
+end
+df_prep = let df = DataFrame(
+        :version => Int[], :model => String[], :symmetry => String[],
+        :D => Int[], :V => Int[], :memory => Int[], :allocs => Int[], :times => Vector{Int}[]
+    )
+
+    for version in versions
+        result = JSON.parsefile(joinpath(resultdir, resultfile(version)))
+        for (model, model_res) in result.data.derivatives.data.AC2_preparation.data
+            for (symmetry, sym_res) in model_res.data
+                for (DV, bench) in sym_res.data
+                    D, V = eval(Meta.parse(DV))::Tuple{Int, Int}
+
+                    push!(
+                        df,
+                        (version, model, symmetry, D, V, bench.memory, bench.allocs, collect(Int, bench.times))
+                    )
+                end
+            end
+        end
+    end
+    df
+end
+
+# Plotting the results
+# --------------------
+fontsize = 20
+estimator = median
+
+f_times = let f = Figure(; size = (1400, 1400))
+    models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+    symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+    df_model = groupby(df_contract, [:model, :symmetry])
+    for row in eachindex(models), col in eachindex(symmetries)
+        df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+        ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "Δt (μs)", yscale = log10)
+        @assert !isnothing(df_data)
+        for (k, v) in pairs(groupby(df_data, :version))
+            Ds = v[!, :D]
+            times = estimator.(v[!, :times]) ./ 1.0e3
+            I = sortperm(Ds)
+            scatterlines!(ax, Ds[I], times[I]; label = "v$(k.version)")
+        end
+        axislegend(ax, position = :lt)
+    end
+
+    Label(f[0, 0], "times"; fontsize)
+    for (row, model) in enumerate(models)
+        Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+    end
+    for (col, symmetry) in enumerate(symmetries)
+        Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+    end
+
+    f
+end
+save(joinpath(resultdir, "bench_times.png"), f_times)
+
+f_times_relative = let f = Figure(; size = (1400, 1400))
+    models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+    symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+    df_model = groupby(df_contract, [:model, :symmetry])
+    for row in eachindex(models), col in eachindex(symmetries)
+        df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+        ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "Δt / Δt₀")
+        hlines!([1], color = :red)
+        @assert !isnothing(df_data)
+
+        df_v = groupby(df_data, :version)
+
+        v = get(df_v, (; version = 0), nothing)
+        Ds = v[!, :D]
+        times = estimator.(v[!, :times])
+        I = sortperm(Ds)
+        times₀ = times[I]
+
+        for (k, v) in pairs(groupby(df_data, :version))
+            k.version == 0 && continue
+            Ds = v[!, :D]
+            I = sortperm(Ds)
+            times = estimator.(v[!, :times])[I]
+            scatterlines!(ax, Ds[I], times ./ times₀; label = "v$(k.version)")
+        end
+        axislegend(ax, position = :lt)
+    end
+
+    Label(f[0, 0], "times"; fontsize)
+    for (row, model) in enumerate(models)
+        Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+    end
+    for (col, symmetry) in enumerate(symmetries)
+        Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+    end
+
+    f
+end
+save(joinpath(resultdir, "bench_times_relative.png"), f_times_relative)
+
+f_allocs = let f = Figure(; size = (1400, 1400))
+    models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+    symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+    df_model = groupby(df_contract, [:model, :symmetry])
+    for row in eachindex(models), col in eachindex(symmetries)
+        df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+        ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "allocs", yscale = log10)
+        @assert !isnothing(df_data)
+        for (k, v) in pairs(groupby(df_data, :version))
+            Ds = v[!, :D]
+            allocs = estimator.(v[!, :allocs])
+            I = sortperm(Ds)
+            scatterlines!(ax, Ds[I], allocs[I]; label = "v$(k.version)")
+        end
+        axislegend(ax, position = :lt)
+    end
+
+    Label(f[0, 0], "allocs"; fontsize)
+    for (row, model) in enumerate(models)
+        Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+    end
+    for (col, symmetry) in enumerate(symmetries)
+        Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+    end
+
+    f
+end
+save(joinpath(resultdir, "bench_allocs.png"), f_allocs)
+
+f_memory = let f = Figure(; size = (1400, 1400))
+    models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+    symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+    df_model = groupby(df_contract, [:model, :symmetry])
+    for row in eachindex(models), col in eachindex(symmetries)
+        df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+        ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "memory (KiB)", yscale = log10)
+        @assert !isnothing(df_data)
+        for (k, v) in pairs(groupby(df_data, :version))
+            Ds = v[!, :D]
+            memory = estimator.(v[!, :memory]) ./ (2^10)
+            I = sortperm(Ds)
+            scatterlines!(ax, Ds[I], memory[I]; label = "v$(k.version)")
+        end
+        axislegend(ax, position = :lt)
+    end
+
+    Label(f[0, 0], "memory"; fontsize)
+    for (row, model) in enumerate(models)
+        Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+    end
+    for (col, symmetry) in enumerate(symmetries)
+        Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+    end
+
+    f
+end
+save(joinpath(resultdir, "bench_memory.png"), f_allocs)
+
+f_memory_relative = let f = Figure(; size = (1400, 1400))
+    models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+    symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+    df_model = groupby(df_contract, [:model, :symmetry])
+    for row in eachindex(models), col in eachindex(symmetries)
+        df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+        ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "memory / memory₀")
+        hlines!([1], color = :red)
+        @assert !isnothing(df_data)
+
+        df_v = groupby(df_data, :version)
+
+        v = get(df_v, (; version = 0), nothing)
+        Ds = v[!, :D]
+        times = estimator.(v[!, :memory])
+        I = sortperm(Ds)
+        times₀ = times[I]
+
+        for (k, v) in pairs(groupby(df_data, :version))
+            k.version == 0 && continue
+            Ds = v[!, :D]
+            I = sortperm(Ds)
+            times = estimator.(v[!, :memory])[I]
+            scatterlines!(ax, Ds[I], times ./ times₀; label = "v$(k.version)")
+        end
+        axislegend(ax, position = :lt)
+    end
+
+    Label(f[0, 0], "memory (relative)"; fontsize)
+    for (row, model) in enumerate(models)
+        Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+    end
+    for (col, symmetry) in enumerate(symmetries)
+        Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+    end
+
+    f
+end
+save(joinpath(resultdir, "bench_memory_relative.png"), f_memory_relative)
+
+
+# Including preparation times
+# ---------------------------
+for n_applications in [3, 10, 30]
+    f_times_relative = let f = Figure(; size = (1400, 1400))
+        models = ["heisenberg_nn", "heisenberg_nnn", "heisenberg_cylinder", "heisenberg_coulomb"]
+        symmetries = ["Trivial", "Irrep[U₁]", "Irrep[SU₂]"]
+
+
+        df_model = groupby(df_contract, [:model, :symmetry])
+        dfp_model = groupby(df_prep, [:model, :symmetry])
+        for row in eachindex(models), col in eachindex(symmetries)
+            df_data = get(df_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+            dfp_data = get(dfp_model, (; model = models[row], symmetry = symmetries[col]), nothing)
+            ax = Axis(f[row, col], xscale = log10, xlabel = "D", ylabel = "Δt / Δt₀")
+            hlines!([1], color = :red)
+            @assert !isnothing(df_data) && !isnothing(dfp_data)
+
+            df_v = groupby(df_data, :version)
+            dfp_v = groupby(dfp_data, :version)
+
+            v = get(df_v, (; version = 0), nothing)
+            Ds = v[!, :D]
+            times = estimator.(v[!, :times])
+            I = sortperm(Ds)
+            times₀ = n_applications .* times[I]
+
+            vp = get(dfp_v, (; version = 0), nothing)
+            Ds = vp[!, :D]
+            times = estimator.(vp[!, :times])
+            I = sortperm(Ds)
+            times₀ .+= times[I]
+
+            df_data_v = groupby(dfp_data, :version)
+            for (k, v) in pairs(groupby(df_data, :version))
+                k.version == 0 && continue
+                Ds = v[!, :D]
+                I = sortperm(Ds)
+                times = n_applications .* estimator.(v[!, :times])[I]
+
+                vp = get(df_data_v, (; k.version), nothing)
+                @assert !isnothing(vp)
+                Ds = vp[!, :D]
+                I = sortperm(Ds)
+                times .+= estimator.(vp[!, :times][I])
+
+                scatterlines!(ax, Ds[I], times ./ times₀; label = "v$(k.version)")
+            end
+            axislegend(ax, position = :lt)
+        end
+
+        Label(f[0, 0], "times"; fontsize)
+        for (row, model) in enumerate(models)
+            Label(f[row, 0], model; rotation = pi / 2, fontsize, tellheight = false, tellwidth = false)
+        end
+        for (col, symmetry) in enumerate(symmetries)
+            Label(f[0, col], symmetry; fontsize, tellheight = false, tellwidth = false)
+        end
+
+        f
+    end
+    save(joinpath(resultdir, "bench_prep_times_relative_n=$n_applications.png"), f_times_relative)
+end
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index cd465643d..96ed23cb4 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -60,11 +60,13 @@ using Compat: @compat
 # -------
 using TensorKit
 using TensorKit: BraidingTensor
+using TensorKit: TupleTools as TT
 using MatrixAlgebraKit
 using MatrixAlgebraKit: TruncationStrategy, PolarViaSVD, LAPACK_SVDAlgorithm
 using BlockTensorKit
 using BlockTensorKit: TensorMapSumSpace
 using TensorOperations
+using TensorOperations: AbstractBackend, DefaultBackend, DefaultAllocator
 using KrylovKit
 using KrylovKit: KrylovAlgorithm
 using OptimKit
@@ -97,6 +99,7 @@ include("utility/logging.jl")
 using .IterativeLoggers
 include("utility/iterativesolvers.jl")
 
+include("utility/allocator.jl")
 include("utility/styles.jl")
 include("utility/periodicarray.jl")
 include("utility/windowarray.jl")
diff --git a/src/algorithms/derivatives/derivatives.jl b/src/algorithms/derivatives/derivatives.jl
index 1f3980e37..4eac14c40 100644
--- a/src/algorithms/derivatives/derivatives.jl
+++ b/src/algorithms/derivatives/derivatives.jl
@@ -213,3 +213,19 @@ const DerivativeOrMultiplied{D <: DerivativeOperator} = Union{MultipliedOperator
 (x::LazySum{<:DerivativeOrMultiplied})(y, t::Number) = sum(O -> O(y, t), x)
 (x::LazySum{<:DerivativeOrMultiplied})(y) = sum(O -> O(y), x)
 Base.:*(h::LazySum{<:Union{DerivativeOrMultiplied}}, v) = h(v)
+
+# Operator preparation
+# --------------------
+"""
+    prepare_operator!!(O, [backend], [allocator]) -> O′
+
+Given an operator and vector, try to construct a more efficient representation of that operator for repeated application.
+This should always be used in conjunction with [`unprepare_operator!!`](@ref).
+"""
+prepare_operator!!(O, backend::AbstractBackend = DefaultBackend(), allocator = GrowingBuffer()) = O
+
+# to make benchmark scripts run
+prepare_operator!!(O, x::AbstractTensorMap, backend::AbstractBackend = DefaultBackend(), allocator = GrowingBuffer()) =
+    prepare_operator!!(O, backend, allocator), x
+unprepare_operator!!(y, O, x, backend::AbstractBackend = DefaultBackend(), allocator = GrowingBuffer()) =
+    y
diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl
index b922205f8..5015575e5 100644
--- a/src/algorithms/derivatives/hamiltonian_derivatives.jl
+++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl
@@ -1,3 +1,11 @@
+const _HAM_MPS_TYPES = Union{
+    FiniteMPS{<:MPSTensor},
+    WindowMPS{<:MPSTensor},
+    InfiniteMPS{<:MPSTensor},
+}
+
+# Single site derivative
+# ----------------------
 """
     JordanMPO_AC_Hamiltonian{O1,O2,O3}
 
@@ -10,292 +18,326 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator
     E::Union{O1, Missing} # finished
     C::Union{O2, Missing} # starting
     B::Union{O2, Missing} # ending
-    A::O3                # continuing
-    function JordanMPO_AC_Hamiltonian(
-            onsite, not_started, finished, starting, ending, continuing
-        )
-        # obtaining storagetype of environments since these should have already mixed
-        # the types of the operator and state
-        gl = continuing[1]
-        S = spacetype(gl)
-        M = storagetype(gl)
-        O1 = tensormaptype(S, 1, 1, M)
-        O2 = tensormaptype(S, 2, 2, M)
-        return new{O1, O2, typeof(continuing)}(
-            onsite, not_started, finished, starting, ending, continuing
-        )
-    end
-    function JordanMPO_AC_Hamiltonian{O1, O2, O3}(
-            onsite, not_started, finished, starting, ending, continuing
-        ) where {O1, O2, O3}
-        return new{O1, O2, O3}(onsite, not_started, finished, starting, ending, continuing)
-    end
+    A::Union{O3, Missing} # continuing
 end
 
-"""
-    JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4}
-
-Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
-In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
-"""
-struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator
-    II::Union{O1, Missing} # not_started
-    IC::Union{O2, Missing} # starting right
-    ID::Union{O1, Missing} # onsite right
-    CB::Union{O2, Missing} # starting left - ending right
-    CA::Union{O3, Missing} # starting left - continuing right
-    AB::Union{O3, Missing} # continuing left - ending right
-    AA::O4                 # continuing left - continuing right
-    BE::Union{O2, Missing} # ending left
-    DE::Union{O1, Missing} # onsite left
-    EE::Union{O1, Missing} # finished
-    function JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
-        # obtaining storagetype of environments since these should have already mixed
-        # the types of the operator and state
-        gl = AA[1]
-        S = spacetype(gl)
-        M = storagetype(gl)
-        O1 = tensormaptype(S, 1, 1, M)
-        O2 = tensormaptype(S, 2, 2, M)
-        O3 = tensormaptype(S, 3, 3, M)
-        return new{O1, O2, O3, typeof(AA)}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
-    end
-    function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(
-            II, IC, ID, CB, CA, AB, AA, BE, DE, EE
-        ) where {O1, O2, O3, O4}
-        return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
-    end
-end
-
-# Constructors
-# ------------
-const _HAM_MPS_TYPES = Union{
-    FiniteMPS{<:MPSTensor},
-    WindowMPS{<:MPSTensor},
-    InfiniteMPS{<:MPSTensor},
-}
-
 function AC_hamiltonian(
-        site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor},
-        above::_HAM_MPS_TYPES, envs
+        site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs
     )
+    @assert below === above "JordanMPO assumptions break"
     GL = leftenv(envs, site, below)
     GR = rightenv(envs, site, below)
     W = operator[site]
+    return JordanMPO_AC_Hamiltonian(GL, W, GR)
+end
+
+function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor)
+    # onsite
+    D = nonzero_length(W.D) > 0 ? only(W.D) : missing
+
+    # not started
+    I = size(W, 4) == 1 ? missing : removeunit(GR[1], 2)
+
+    # finished
+    E = size(W, 1) == 1 ? missing : removeunit(GL[end], 2)
 
     # starting
-    if nonzero_length(W.C) > 0
+    C = if nonzero_length(W.C) > 0
         GR_2 = GR[2:(end - 1)]
-        @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2]
-        starting = only(starting_)
+        @plansor starting[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2]
+        only(starting)
     else
-        starting = missing
+        missing
     end
 
     # ending
-    if nonzero_length(W.B) > 0
+    B = if nonzero_length(W.B) > 0
         GL_2 = GL[2:(end - 1)]
-        @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4]
-        ending = only(ending_)
+        @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4]
+        only(ending)
     else
-        ending = missing
+        missing
     end
 
+    # continuing
+    A = MPO_AC_Hamiltonian(GL[2:(end - 1)], W.A, GR[2:(end - 1)])
+
+    # obtaining storagetype of environments since these should have already mixed
+    # the types of the operator and state
+    S = spacetype(GL)
+    M = storagetype(GL)
+    O1 = tensormaptype(S, 1, 1, M)
+    O2 = tensormaptype(S, 2, 2, M)
+    O3 = typeof(A)
+
+    # specialization for nearest neighbours
+    nonzero_length(W.A) == 0 && (A = missing)
+
+    return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A)
+end
+
+function prepare_operator!!(
+        H::JordanMPO_AC_Hamiltonian{O1, O2, O3}, backend::AbstractBackend, allocator
+    ) where {O1, O2, O3}
+    C = H.C
+    B = H.B
+
     # onsite
-    if nonzero_length(W.D) > 0
-        if !ismissing(starting)
-            @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2]
-            onsite = missing
-        elseif !ismissing(ending)
-            @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4]
-            onsite = missing
-        else
-            onsite = W.D
-        end
+    D = if !ismissing(C)
+        Id = TensorKit.id(storagetype(W.D), space(C, 2))
+        @plansor C[-1 -2; -3 -4] += W.D[-1; -3] * Id[-4; -2]
+        missing
+    elseif !ismissing(B)
+        Id = TensorKit.id(storagetype(W.D), space(B, 1))
+        @plansor B[-1 -2; -3 -4] += Id[-1; -3] * W.D[-2; -4]
+        missing
     else
-        onsite = missing
+        W.D
     end
 
     # not_started
-    if isfinite(operator) && site == length(operator)
-        not_started = missing
-    elseif !ismissing(starting)
-        I = id(storagetype(GR[1]), physicalspace(W))
-        @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2]
-        not_started = missing
+    I = if !ismissing(C)
+        Id = id(storagetype(W.I), space(C, 1))
+        @plansor C[-1 -2; -3 -4] += I[-1; -3] * H.I[-4; -2]
+        missing
     else
-        not_started = removeunit(GR[1], 2)
+        H.I
     end
 
     # finished
-    if isfinite(operator) && site == 1
-        finished = missing
-    elseif !ismissing(ending)
-        I = id(storagetype(GL[end]), physicalspace(W))
-        @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4]
-        finished = missing
+    E = if !ismissing(B)
+        Id = id(storagetype(W.I), space(B, 2))
+        @plansor B[-1 -2; -3 -4] += H.E[-1; -3] * Id[-2; -4]
+        missing
     else
-        finished = removeunit(GL[end], 2)
+        H.E
     end
 
-    # continuing
-    A = W.A
-    continuing = (GL[2:(end - 1)], A, GR[2:(end - 1)])
+    O3′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O3, typeof(backend), typeof(allocator)})
+    A = ismissing(H.A) ? H.A : prepare_operator!!(H.A, backend, allocator)
 
-    return JordanMPO_AC_Hamiltonian(
-        onsite, not_started, finished, starting, ending, continuing
-    )
+    return JordanMPO_AC_Hamiltonian{O1, O2, O3′}(D, I, E, C, B, A)
+end
+
+
+# Two site derivative
+# -------------------
+"""
+    JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4}
+
+Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
+In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
+"""
+struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator
+    II::Union{O1, Missing} # not_started
+    IC::Union{O2, Missing} # starting right
+    ID::Union{O1, Missing} # onsite right
+    CB::Union{O2, Missing} # starting left - ending right
+    CA::Union{O3, Missing} # starting left - continuing right
+    AB::Union{O3, Missing} # continuing left - ending right
+    AA::Union{O4, Missing} # continuing left - continuing right
+    BE::Union{O2, Missing} # ending left
+    DE::Union{O1, Missing} # onsite left
+    EE::Union{O1, Missing} # finished
 end
 
 function AC2_hamiltonian(
-        site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor},
-        above::_HAM_MPS_TYPES, envs
+        site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs
     )
+    @assert below === above "JordanMPO assumptions break"
     GL = leftenv(envs, site, below)
     GR = rightenv(envs, site + 1, below)
-    W1 = operator[site]
-    W2 = operator[site + 1]
+    W1, W2 = operator[site], operator[site + 1]
+    return JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR)
+end
+
+function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor)
+    # not started
+    II = size(W2, 4) == 1 ? missing : transpose(removeunit(GR[1], 2))
+
+    # finished
+    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]
+        only(IC_)
+    else
+        missing
+    end
+
+    # onsite left
+    DE = nonzero_length(W1.D) > 0 ? only(W1.D) : missing
+
+    # onsite right
+    ID = nonzero_length(W2.D) > 0 ? only(W2.D) : 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]
+        only(CB_)
+    else
+        missing
+    end
 
     # starting left - continuing right
-    if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0
+    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] *
             GR[2:(end - 1)][-6 1; -3]
-        CA = only(CA_)
+        only(CA_)
     else
-        CA = missing
+        missing
     end
 
     # continuing left - ending right
-    if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0
+    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 = only(AB_)
+        only(AB_)
     else
-        AB = missing
+        missing
     end
 
-    # middle
-    if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0
-        if !ismissing(CA)
-            @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] *
-                removeunit(GR[end], 2)[-6; -3]
-            CB = missing
-        elseif !ismissing(AB)
-            @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] *
-                W1.C[-2; -5 1] * W2.B[1 -3; -6]
-            CB = missing
-        else
-            @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4]
-            CB = only(CB_)
+    # 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]
+        only(BE_)
+    else
+        missing
+    end
+
+    # continuing - continuing
+    AA = MPO_AC2_Hamiltonian(GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)])
+
+    S = spacetype(GL)
+    M = storagetype(GL)
+    O1 = tensormaptype(S, 1, 1, M)
+    O2 = tensormaptype(S, 2, 2, M)
+    O3 = tensormaptype(S, 3, 3, M)
+    O4 = typeof(AA)
+
+    if nonzero_length(W1.A) == 0 && nonzero_length(W2.A) == 0
+        AA = missing
+    else
+        mask1 = falses(size(W1.A, 1), size(W1.A, 4))
+        for I in nonzero_keys(W1.A)
+            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[I[1], I[4]] = true
         end
+
+        mask_left = transpose(mask1) * trues(size(mask1, 1))
+        mask_right = mask2 * trues(size(mask2, 2))
+        all(iszero, mask_left .* mask_right) && (AA = missing)
+    end
+
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
+
+end
+
+function prepare_operator!!(
+        H::JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}, backend::AbstractBackend, allocator
+    ) where {O1, O2, O3, O4}
+
+    CA = H.CA
+    AB = H.AB
+
+    CB = if !ismissing(CA) && !ismissing(H.CB)
+        Id = TensorKit.id(storagetype(H.CB), space(CA, 3))
+        @plansor CA[-1 -2 -3; -4 -5 -6] += H.CB[-1 -2; -4 -5] * Id[-3; -6]
+        missing
+    elseif !ismissing(AB) && !ismissing(H.CB)
+
     else
-        CB = missing
+        H.CB
     end
 
     # starting right
-    if nonzero_length(W2.C) > 0
-        if !ismissing(CA)
-            I = id(storagetype(GR[1]), physicalspace(W1))
-            @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) *
-                GR[2:(end - 1)][-6 1; -3]
-            IC = missing
-        else
-            @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2]
-        end
+    IC = if !ismissing(CA) && !ismissing(H.IC)
+        Id = TensorKit.id(storagetype(H.IC), space(CA, 1))
+        @plansor CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * H.IC[ -2 -3; -5 -6]
+        missing
     else
-        IC = missing
+        H.IC
     end
 
     # ending left
-    if nonzero_length(W1.B) > 0
-        if !ismissing(AB)
-            I = id(storagetype(GL[end]), physicalspace(W2))
-            @plansor AB[-1 -2 -3; -4 -5 -6] += GL[2:(end - 1)][-1 1; -4] *
-                (W1.B[1 -2; -5] * I[-3; -6])
-            BE = missing
-        else
-            @plansor BE[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4]
-        end
+    BE = if !ismissing(AB) && !ismissing(H.BE)
+        Id = TensorKit.id(storagetype(H.BE), space(AB, 3))
+        @plansor AB[-1 -2 -3; -4 -5 -6] += H.BE[-1 -2; -4 -5] * Id[-3; -6]
+        missing
     else
-        BE = missing
+        H.BE
     end
 
     # onsite left
-    if nonzero_length(W1.D) > 0
-        if !ismissing(BE)
-            @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4]
-            DE = missing
-        elseif !ismissing(AB)
-            I = id(storagetype(GL[end]), physicalspace(W2))
-            @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] *
-                (W1.D[-2; -5] * I[-3; -6])
-            DE = missing
-            # TODO: could also try in CA?
-        else
-            DE = only(W1.D)
-        end
+    DE = if !ismissing(BE) && !ismissing(H.DE)
+        Id = TensorKit.id(storagetype(H.DE), space(BE, 1))
+        @plansor BE[-1 -2; -3 -4] += Id[-1; -3] * H.DE[-2; -4]
+        missing
+    elseif !ismissing(AB) && !ismissing(H.DE)
+        Id1 = id(storagetype(H.DE), space(AB, 1))
+        Id2 = id(storagetype(H.DE), space(AB, 3))
+        @plansor AB[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.DE[-2; -5] * Id2[-3; -6]
+        missing
+        # TODO: could also try in CA?
     else
-        DE = missing
+        H.DE
     end
 
     # onsite right
-    if nonzero_length(W2.D) > 0
-        if !ismissing(IC)
-            @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2]
-            ID = missing
-        elseif !ismissing(CA)
-            I = id(storagetype(GR[1]), physicalspace(W1))
-            @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) *
-                removeunit(GR[end], 2)[-6; -3]
-            ID = missing
-        else
-            ID = only(W2.D)
-        end
+    ID = if !ismissing(IC) && !ismissing(H.ID)
+        Id = TensorKit.id(storagetype(H.ID), space(IC, 2))
+        @plansor IC[-1 -2; -3 -4] += H.ID[-1; -3] * Id[-2; -4]
+        missing
+    elseif !ismissing(CA) && !ismissing(H.ID)
+        Id1 = TensorKit.id(storagetype(H.ID), space(CA, 1))
+        Id2 = TensorKit.id(storagetype(H.ID), space(CA, 3))
+        @plansor CA[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.ID[-2; -5] * Id2[-3; -6]
+        missing
     else
-        ID = missing
+        H.ID
     end
 
     # finished
-    if isfinite(operator) && site + 1 == length(operator)
-        II = missing
-    elseif !ismissing(IC)
-        I = id(storagetype(GR[1]), physicalspace(W2))
-        @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2]
+    II = if !ismissing(IC) && !ismissing(H.II)
+        I = id(storagetype(H.II), space(IC, 1))
+        @plansor IC[-1 -2; -3 -4] += I[-1; -3] * H.II[-2; -4]
         II = missing
-    elseif !ismissing(CA)
-        I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2))
-        @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3]
+    elseif !ismissing(CA) && !ismissing(H.II)
+        I = id(storagetype(H.II), space(CA, 1) ⊗ space(CA, 2))
+        @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * H.II[-3; -6]
         II = missing
     else
-        II = transpose(removeunit(GR[1], 2))
+        H.II
     end
 
     # unstarted
-    if isfinite(operator) && site == 1
+    EE = if !ismissing(BE) && !ismissing(H.EE)
+        I = id(storagetype(H.EE), space(BE, 2))
+        @plansor BE[-1 -2; -3 -4] += H.EE[-1; -3] * I[-2; -4]
         EE = missing
-    elseif !ismissing(BE)
-        I = id(storagetype(GL[end]), physicalspace(W1))
-        @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4]
-        EE = missing
-    elseif !ismissing(AB)
-        I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2))
-        @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6]
+    elseif !ismissing(AB) && !ismissing(H.EE)
+        I = id(storagetype(H.EE), space(AB, 2) ⊗ space(AB, 3))
+        @plansor AB[-1 -2 -3; -4 -5 -6] += H.EE[-1; -4] * I[-2 -3; -5 -6]
         EE = missing
     else
-        EE = removeunit(GL[end], 2)
+        H.EE
     end
 
-    # continuing - continuing
-    # TODO: MPODerivativeOperator code reuse + optimization
-    AA = (GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)])
+    O4′ = Core.Compiler.return_type(prepare_operator!!, Tuple{O4, typeof(backend), typeof(allocator)})
+    AA = prepare_operator!!(H.AA, backend, allocator)
 
-    return JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4′}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
 end
 
 # Actions
 # -------
 function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor)
-    y = zerovector(x)
+    y = ismissing(H.A) ? zerovector(x) : H.A(x)
 
     ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1]
     ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3]
@@ -303,16 +345,12 @@ function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor)
     ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1]
     ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3]
 
-    GL, A, GR = H.A
-    if nonzero_length(A) > 0
-        @plansor y[-1 -2; -3] += GL[-1 5; 4] * x[4 2; 1] * A[5 -2; 2 3] * GR[1 3; -3]
-    end
-
     return y
 end
 
 function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor)
-    y = zerovector(x)
+    y = ismissing(H.AA) ? zerovector(x) : H.AA(x)
+
     ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1]
     ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1]
     ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1]
@@ -323,12 +361,5 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor)
     ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1]
     ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1]
 
-    GL, A1, A2, GR = H.AA
-    if nonzero_length(A1) > 0 && nonzero_length(A2) > 0
-        # TODO: there are too many entries here, this could be further optimized
-        @plansor y[-1 -2; -3 -4] += GL[-1 7; 6] * x[6 5; 1 3] * A1[7 -2; 5 4] *
-            A2[4 -4; 3 2] * GR[1 2; -3]
-    end
-
     return y
 end
diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 12b030b1e..30c8e1f38 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -105,3 +105,90 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen
         h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3]
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
+
+# prepared operators
+# ------------------
+struct PrecomputedDerivative{
+        T <: Number, S <: ElementarySpace, N₁, N₂, N₃, N₄,
+        T1 <: AbstractTensorMap{T, S, N₁, N₂}, T2 <: AbstractTensorMap{T, S, N₃, N₄},
+        B <: AbstractBackend, A,
+    } <: DerivativeOperator
+    leftenv::T1
+    rightenv::T2
+    backend::B
+    allocator::A
+end
+
+const PrecomputedACDerivative{T, S} = PrecomputedDerivative{T, S, 3, 2, 2, 1}
+const PrecomputedAC2Derivative{T, S} = PrecomputedDerivative{T, S, 3, 2, 3, 2}
+
+function prepare_operator!!(
+        H::MPO_C_Hamiltonian{<:MPSTensor, <:MPSTensor},
+        backend::AbstractBackend, allocator
+    )
+    leftenv = _transpose_tail(TensorMap(H.leftenv))
+    rightenv = TensorMap(H.rightenv)
+    return PrecomputedDerivative(leftenv, rightenv, backend, allocator)
+end
+function prepare_operator!!(
+        H::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor},
+        backend::AbstractBackend, allocator
+    )
+    cp = checkpoint(allocator)
+    @plansor backend = backend allocator = allocator begin
+        GL_O[-1 -2; -4 -5 -3] := H.leftenv[-1 1; -4] * H.operators[1][1 -2; -5 -3]
+    end
+    reset!(allocator, cp)
+    leftenv = fuse_legs(TensorMap(GL_O), 0, 2)
+    rightenv = TensorMap(H.rightenv)
+
+    return PrecomputedDerivative(leftenv, rightenv, backend, allocator)
+end
+
+function prepare_operator!!(
+        H::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor},
+        backend::AbstractBackend, allocator
+    )
+    cp = checkpoint(allocator)
+    @plansor backend = backend allocator = allocator begin
+        GL_O[-1 -2; -4 -5 -3] := H.leftenv[-1 1; -4] * H.operators[1][1 -2; -5 -3]
+        O_GR[-1 -2 -3; -4 -5] := H.operators[2][-3 -5; -2 1] * H.rightenv[-1 1; -4]
+    end
+    reset!(allocator, cp)
+
+    leftenv = fuse_legs(GL_O isa TensorMap ? GL_O : TensorMap(GL_O), 0, 2)
+    rightenv = fuse_legs(O_GR isa TensorMap ? O_GR : TensorMap(O_GR), 2, 0)
+    return PrecomputedDerivative(leftenv, rightenv, backend, allocator)
+end
+
+
+function (H::PrecomputedDerivative)(x::AbstractTensorMap)
+    allocator = H.allocator
+    cp = checkpoint(allocator)
+
+    R_fused = fuse_legs(H.rightenv, 0, numin(x))
+    x_fused = fuse_legs(x, numout(x), numin(x))
+
+    TC = TensorOperations.promote_contract(scalartype(x_fused), scalartype(R_fused))
+    xR = TensorOperations.tensoralloc_contract(
+        TC, x_fused, ((1,), (2,)), false, R_fused, ((1,), (2, 3)), false, ((1, 2), (3,)), Val(true), allocator
+    )
+
+    mul_front!(xR, x_fused, R_fused, One(), Zero(), H.backend, allocator)
+
+    LxR = H.leftenv * xR
+    TensorOperations.tensorfree!(xR, allocator)
+
+    reset!(allocator, cp)
+    return TensorMap{scalartype(LxR)}(LxR.data, codomain(H.leftenv) ← domain(H.rightenv))
+end
+
+const _ToPrepare = Union{
+    MPO_C_Hamiltonian{<:MPSTensor, <:MPSTensor},
+    MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor},
+    MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor},
+}
+
+function prepare_operator!!(H::Multiline{<:_ToPrepare}, backend::AbstractBackend, allocator)
+    return Multiline(map(x -> prepare_operator!!(x, backend, allocator), parent(H)))
+end
diff --git a/src/algorithms/fixedpoint.jl b/src/algorithms/fixedpoint.jl
index 81e7c69d9..fa29f80b2 100644
--- a/src/algorithms/fixedpoint.jl
+++ b/src/algorithms/fixedpoint.jl
@@ -10,9 +10,8 @@ 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
+    info.converged == 0 &&
         @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
-    end
 
     return vals[1], vecs[1]
 end
@@ -20,12 +19,9 @@ end
 function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi)
     TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg)
 
-    if info.converged == 0
+    info.converged == 0 &&
         @warnv 1 "fixedpoint 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"
-    end
+    size(TT, 2) > 1 && TT[2, 1] != 0 && @warnv 1 "non-unique fixedpoint detected"
 
     return vals[1], vecs[1]
 end
diff --git a/src/utility/allocator.jl b/src/utility/allocator.jl
new file mode 100644
index 000000000..f9dc1a482
--- /dev/null
+++ b/src/utility/allocator.jl
@@ -0,0 +1,67 @@
+@static if isdefined(Core, :Memory)
+    BufType = Memory{UInt8}
+else
+    BufType = Vector{UInt8}
+end
+
+const DEFAULT_SIZEHINT = 2^20 # 1MB
+
+mutable struct GrowingBuffer
+    buffer::BufType
+    offset::UInt
+    function GrowingBuffer(; sizehint = DEFAULT_SIZEHINT)
+        buffer = BufType(undef, sizehint)
+        return new(buffer, zero(UInt))
+    end
+end
+
+Base.length(buffer::GrowingBuffer) = length(buffer.buffer)
+Base.pointer(buffer::GrowingBuffer) = pointer(buffer.buffer) + buffer.offset
+
+function Base.sizehint!(buffer::GrowingBuffer, n::Integer; shrink::Bool = false)
+    n > 0 || throw(ArgumentError("invalid new buffer size"))
+    buffer.offset == 0 || error("cannot resize a buffer that is not fully reset")
+
+    n = shrink ? max(n, length(buffer)) : n
+    n = Int(Base.nextpow(2, n))
+
+    @static if isdefined(Core, :Memory)
+        buffer.buffer = BufType(undef, n)
+    else
+        sizehint!(buffer.buffer, n)
+    end
+    return buffer
+end
+
+checkpoint(buffer) = zero(UInt)
+reset!(buffer, checkpoint::UInt = zero(UInt)) = buffer
+
+checkpoint(buffer::GrowingBuffer) = buffer.offset
+
+function reset!(buffer::GrowingBuffer, checkpoint::UInt = zero(UInt))
+    if iszero(checkpoint) && buffer.offset > length(buffer)
+        # full reset - check for need to grow
+        newlength = Base.nextpow(2, buffer.offset) # round to nearest larger power of 2
+        buffer.offset = checkpoint
+        sizehint!(buffer, newlength)
+    else
+        buffer.offset = checkpoint
+    end
+    return buffer
+end
+
+# Allocating
+# ----------
+function TensorOperations.tensoralloc(
+        ::Type{A}, structure, ::Val{istemp}, buffer::GrowingBuffer
+    ) where {A <: AbstractArray, istemp}
+    T = eltype(A)
+    if istemp
+        ptr = convert(Ptr{T}, pointer(buffer))
+        buffer.offset += prod(structure) * sizeof(T)
+        buffer.offset < length(buffer) &&
+            return Base.unsafe_wrap(Array, ptr, structure)
+    end
+    return A(undef, structure)
+end
+TensorOperations.tensorfree!(::AbstractArray, ::GrowingBuffer) = nothing
diff --git a/src/utility/multiline.jl b/src/utility/multiline.jl
index a1c861ff1..ebcff438f 100644
--- a/src/utility/multiline.jl
+++ b/src/utility/multiline.jl
@@ -12,7 +12,7 @@ See also: [`MultilineMPS`](@ref) and [`MultilineMPO`](@ref)
 struct Multiline{T}
     data::PeriodicArray{T, 1}
     function Multiline{T}(data::AbstractVector{T}) where {T}
-        @assert allequal(length.(data)) "All lines must have the same length"
+        # @assert allequal(length.(data)) "All lines must have the same length"
         return new{T}(data)
     end
 end
@@ -22,7 +22,7 @@ Multiline(data::AbstractVector{T}) where {T} = Multiline{T}(data)
 # -----------------------
 Base.parent(m::Multiline) = m.data
 Base.size(m::Multiline) = (length(parent(m)), length(parent(m)[1]))
-Base.size(m::Multiline, i::Int) = getindex(size(m), i)
+Base.size(m::Multiline, i::Int) = i == 1 ? length(parent(m)) : i == 2 ? length(parent(m)[1]) : error()
 Base.length(m::Multiline) = prod(size(m))
 function Base.axes(m::Multiline, i::Int)
     return i == 1 ? axes(parent(m), 1) :
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index 49cb96648..d4ce0f6cf 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -8,8 +8,8 @@ function _transpose_as(t1::AbstractTensorMap, t2::AbstractTensorMap; copy::Bool
     return repartition(t1, numout(t2), numin(t2); copy)
 end
 
-_mul_front(C, A) = _transpose_front(C * _transpose_tail(A))
-_mul_tail(A, C) = A * C
+_mul_front(C, A) = mul_front(C, A) # _transpose_front(C * _transpose_tail(A))
+_mul_tail(A, C) = mul_tail(A, C) # A * C
 
 function _similar_tail(A::AbstractTensorMap)
     cod = _firstspace(A)
@@ -149,3 +149,206 @@ function check_unambiguous_braiding(V::VectorSpace)
     return check_unambiguous_braiding(Bool, V) ||
         throw(ArgumentError("cannot unambiguously braid $V"))
 end
+
+"""
+    matrix_contract(
+        A::AbstractTensorMap, B::AbstractTensorMap{T, S, 1, 1}, i::Int,
+        α::Number = One(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator();
+        transpose::Bool = false
+    )
+
+Compute the tensor contraction `α * A * B`, where the (1, 1) - tensor `B` is attached to index `i` of `A`.
+Whenever `transpose = true`, this contraction (lazily) uses `transpose(B)` instead.
+
+See also [`matrix_contract!`](@ref).
+"""
+function matrix_contract(
+        A::AbstractTensorMap, B::AbstractTensorMap{<:Any, <:Any, 1, 1}, i::Int,
+        α::Number = One(), backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator();
+        transpose::Bool = false
+    )
+    if i <= numout(A)
+        cod = ProductSpace(TT.setindex(codomain(A).spaces, space(B, transpose ? 1 : 2), i))
+        dom = domain(A)
+    else
+        cod = codomain(A)
+        dom = ProductSpace(TT.setindex(domain(A).spaces, space(B, transpose ? 1 : 2)', i - numout(A)))
+    end
+    T = TensorOperations.promote_contract(scalartype(A), scalartype(B), scalartype(α))
+    C = similar(A, T, cod ← dom)
+    return matrix_contract!(C, A, B, i, α, Zero(), backend, allocator; transpose)
+end
+
+"""
+    matrix_contract!(
+        C::AbstractTensorMap, A::AbstractTensorMap, B::AbstractTensorMap{T, S, 1, 1}, i::Int,
+        α::Number = One(), β::Number = Zero(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator();
+        transpose::Bool = false
+    )
+
+Compute the tensor contraction `C ← β * C + α * A * B`, where the (1, 1) - tensor `B` is attached to index `i` of `A`,
+and the result is added into `C`. Whenever `transpose = true`, this contraction (lazily) uses `transpose(B)` instead.
+
+See also [`matrix_contract`](@ref).
+"""
+function matrix_contract!(
+        C::AbstractTensorMap, A::AbstractTensorMap, B::AbstractTensorMap{<:Any, <:Any, 1, 1}, i::Int,
+        α::Number = One(), β::Number = Zero(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator();
+        transpose::Bool = false
+    )
+
+    @boundscheck for k in 1:numind(C)
+        numin(C) == numin(A) && numout(C) == numout(A) || throw(ArgumentError("Invalid number of dimensions"))
+        if k == i
+            space(C, k) == space(B, transpose ? 1 : 2) || throw(SpaceMismatch())
+            space(A, k) == space(B, transpose ? 2 : 1)' || throw(SpaceMismatch())
+        else
+            space(C, k) == space(A, k) || throw(SpaceMismatch())
+        end
+    end
+
+    N, N₁ = numind(C), numout(C)
+    pA = (TT.deleteat(ntuple(identity, N), i), (i,))
+    pB = transpose ? ((2,), (1,)) : ((1,), (2,))
+    pAB = TensorKit._canonicalize(TT.insertafter(ntuple(identity, N - 1), i - 1, (N,)), C)
+
+    Bblocks = blocks(B)
+    for ((f₁, f₂), c) in subblocks(C)
+        uncoupled_i = i <= N₁ ? f₁.uncoupled[i] : f₂.uncoupled[i - N₁]
+        transpose && (uncoupled_i = dual(uncoupled_i))
+        if TensorKit.hasblock(B, uncoupled_i)
+            a = A[f₁, f₂]
+            b = Bblocks[uncoupled_i]
+            TensorOperations.tensorcontract!(c, a, pA, false, b, pB, false, pAB, α, β, backend, allocator)
+        else
+            scale!(c, β)
+        end
+    end
+
+    return C
+end
+
+function mul_front(
+        A::AbstractTensorMap, B::AbstractTensorMap,
+        α::Number = One(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    cod = prod(i -> i == 1 ? space(A, 1) : space(B, i), 1:numout(B))
+    dom = domain(B)
+    T = TensorOperations.promote_contract(scalartype(A), scalartype(B), scalartype(α))
+    C = similar(A, T, cod ← dom)
+    return mul_front!(C, A, B, α, Zero(), backend, allocator)
+end
+
+function mul_front!(
+        C::AbstractTensorMap, A::AbstractTensorMap, B::AbstractTensorMap,
+        α::Number = One(), β::Number = Zero(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    (numin(C) == numin(B) && numout(C) == numout(B) && numin(A) == numout(A) == 1) ||
+        throw(SpaceMismatch())
+
+    numout(B) == 1 && return mul!(C, A, B, α, β)
+
+    cp = checkpoint(allocator)
+
+    Ablocks = blocks(A)
+    Bstructure = TensorKit.fusionblockstructure(space(B))
+    for ((f₁, f₂), c) in subblocks(C)
+        # fetch A block
+        u = first(f₁.uncoupled)
+        a = Ablocks[u]
+        isempty(a) && (scale!(c, β); continue)
+
+        # fetch B block
+        haskey(Bstructure.fusiontreeindices, (f₁, f₂)) || (scale!(c, β); continue)
+        b = B[f₁, f₂]
+
+        tensorcontract!(
+            c,
+            a, ((1,), (2,)), false,
+            b, ((1,), ntuple(i -> i + 1, numind(B) - 1)), false,
+            (ntuple(identity, numout(C)), ntuple(i -> i + numout(C), numin(C))),
+            α, β, backend, allocator
+        )
+    end
+    reset!(allocator, cp)
+
+    return C
+end
+
+function mul_tail(
+        A::AbstractTensorMap, B::AbstractTensorMap,
+        α::Number = One(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    cod = codomain(A)
+    dom = prod(i -> i == 1 ? domain(B)[1] : domain(A)[i], 1:numin(A))
+    T = TensorOperations.promote_contract(scalartype(A), scalartype(B), scalartype(α))
+    C = similar(A, T, cod ← dom)
+    return mul_tail!(C, A, B, α, Zero(), backend, allocator)
+end
+
+function mul_tail!(
+        C::AbstractTensorMap, A::AbstractTensorMap, B::AbstractTensorMap,
+        α::Number = One(), β::Number = Zero(),
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    (numin(C) == numin(A) && numout(C) == numout(A) && numin(B) == numout(B) == 1) ||
+        throw(SpaceMismatch())
+
+    numin(A) == 1 && return mul!(C, A, B, α, β)
+
+    cp = checkpoint(allocator)
+
+    Astructure = TensorKit.fusionblockstructure(space(A))
+    Bblocks = blocks(B)
+    for ((f₁, f₂), c) in subblocks(C)
+        # fetch B block
+        u = first(f₂.uncoupled)
+        b = Bblocks[u]
+        isempty(b) && (scale!(c, β); continue)
+
+        # fetch A block
+        haskey(Astructure.fusiontreeindices, (f₁, f₂)) || (scale!(c, β); continue)
+        a = A[f₁, f₂]
+
+        tensorcontract!(
+            c,
+            a, (ntuple(identity, numind(A) - 1), (1,)), false,
+            b, ((1,), (2,)), false,
+            (ntuple(identity, numout(C)), ntuple(i -> i + numout(C), numin(C))),
+            α, β, backend, allocator
+        )
+    end
+
+    reset!(allocator, cp)
+    return C
+end
+
+@inline fuse_legs(x::TensorMap, N₁::Int, N₂::Int) = fuse_legs(x, Val(N₁), Val(N₂))
+function fuse_legs(x::TensorMap, ::Val{N₁}, ::Val{N₂}) where {N₁, N₂}
+    ((0 <= N₁ <= numout(x)) && (0 <= N₂ <= numin(x))) || throw(ArgumentError("invalid fusing scheme"))
+    init = one(spacetype(x))
+
+    cod = if N₁ > 1
+        cod_spaces = codomain(x).spaces
+        fuse(prod(TT.getindices(cod_spaces, ntuple(identity, N₁)))) ⊗
+            prod(TT.getindices(cod_spaces, ntuple(i -> i + N₁, numout(x) - N₁)); init)
+    else
+        codomain(x)
+    end
+
+    dom = if N₂ > 1
+        dom_spaces = domain(x).spaces
+        dom = fuse(prod(TT.getindices(dom_spaces, ntuple(identity, N₂)); init)) ⊗
+            prod(TT.getindices(domain(x).spaces, ntuple(i -> i + N₂, numin(x) - N₂)); init)
+    else
+        domain(x)
+    end
+
+    return TensorMap{scalartype(x)}(x.data, cod ← dom)
+end
diff --git a/test/algorithms.jl b/test/algorithms.jl
index 207857951..5e41213f9 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -338,7 +338,7 @@ module TestAlgorithms
         H = force_planar(heisenberg_XXX(Trivial, Float64; spin = 1 // 2, L))
         ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)
         E = expectation_value(ψ, H)
-        ψ₀, = find_groundstate(ψ, H)
+        ψ₀, = find_groundstate(ψ, H; verbosity = verbosity_conv)
         E₀ = expectation_value(ψ₀, H)
 
         @testset "Finite $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs