diff --git a/docs/src/man/models.md b/docs/src/man/models.md
index 465ceb2c1..5539fb6f8 100644
--- a/docs/src/man/models.md
+++ b/docs/src/man/models.md
@@ -28,7 +28,7 @@ function MPSKitModels.transverse_field_ising(
     return LocalOperator(
         spaces,
         (neighbor => ZZ for neighbor in nearest_neighbours(lattice))...,
-        ((idx,) => X for idx in vertices(lattice))...,
+        ([idx,] => X for idx in vertices(lattice))...,
     )
 end
 ```
diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl
index 62e32177e..d1de2d3ad 100644
--- a/src/algorithms/contractions/bp_contractions.jl
+++ b/src/algorithms/contractions/bp_contractions.jl
@@ -60,13 +60,51 @@ function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::BPE
 end
 
 function contract_local_operator(
-        inds::NTuple{1, CartesianIndex{2}},
+        inds::Vector{CartesianIndex{2}},
+        O::AbstractTensorMap,
+        ket::InfinitePEPS,
+        bra::InfinitePEPS,
+        env::BPEnv,
+    )
+    length(inds) == 1 && return contract_local_operator1x1(only(inds), O, ket, bra, env)
+
+    if length(inds) == 2
+        ind_relative = inds[2] - inds[1]
+        if ind_relative == CartesianIndex(1, 0)
+            return contract_local_operator2x1(inds[1], O, ket, bra, env)
+        elseif ind_relative == CartesianIndex(0, 1)
+            return contract_local_operator1x2(inds[1], O, ket, bra, env)
+        end
+    end
+    error("No implementation for contractions for BP environments with $inds")
+end
+function contract_local_norm(
+        inds::Vector{CartesianIndex{2}},
+        ket::InfinitePEPS,
+        bra::InfinitePEPS,
+        env::BPEnv,
+    )
+    length(inds) == 1 && return contract_local_norm1x1(only(inds), ket, bra, env)
+
+    if length(inds) == 2
+        ind_relative = inds[2] - inds[1]
+        if ind_relative == CartesianIndex(1, 0)
+            return contract_local_norm2x1(inds[1], ket, bra, env)
+        elseif ind_relative == CartesianIndex(0, 1)
+            return contract_local_norm1x2(inds[1], ket, bra, env)
+        end
+    end
+    error("No implementation for contractions for BP environments with $inds")
+end
+
+function contract_local_operator1x1(
+        ind::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 1, 1},
         ket::InfinitePEPS,
         bra::InfinitePEPS,
         env::BPEnv,
     )
-    row, col = Tuple(only(inds))
+    row, col = Tuple(ind)
     M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)]
     M_east = env.messages[EAST, mod1(row, end), _next(col, end)]
     M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)]
@@ -83,10 +121,10 @@ function contract_local_operator(
     end
 end
 
-function contract_local_norm(
-        inds::NTuple{1, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
+function contract_local_norm1x1(
+        ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
-    row, col = Tuple(only(inds))
+    row, col = Tuple(ind)
     M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)]
     M_east = env.messages[EAST, mod1(row, end), _next(col, end)]
     M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)]
@@ -102,24 +140,7 @@ function contract_local_norm(
     end
 end
 
-function contract_local_operator(
-        inds::NTuple{2, CartesianIndex{2}},
-        O::AbstractTensorMap{<:Any, <:Any, 2, 2},
-        ket::InfinitePEPS,
-        bra::InfinitePEPS,
-        env::BPEnv,
-    )
-    ind_relative = inds[2] - inds[1]
-    return if ind_relative == CartesianIndex(1, 0)
-        contract_vertical_operator(inds[1], O, ket, bra, env)
-    elseif ind_relative == CartesianIndex(0, 1)
-        contract_horizontal_operator(inds[1], O, ket, bra, env)
-    else
-        error("Only contractions for nearest neighbor bonds are implemented.")
-    end
-end
-
-function contract_vertical_operator(
+function contract_local_operator2x1(
         coord::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 2, 2},
         ket::InfinitePEPS,
@@ -147,7 +168,7 @@ function contract_vertical_operator(
         O[dNb dSb; dNt dSt]
 end
 
-function contract_horizontal_operator(
+function contract_local_operator1x2(
         coord::CartesianIndex{2},
         O::AbstractTensorMap{<:Any, <:Any, 2, 2},
         ket::InfinitePEPS,
@@ -181,20 +202,7 @@ function contract_horizontal_operator(
     end
 end
 
-function contract_local_norm(
-        inds::NTuple{2, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
-    )
-    ind_relative = inds[2] - inds[1]
-    return if ind_relative == CartesianIndex(1, 0)
-        contract_vertical_norm(inds[1], ket, bra, env)
-    elseif ind_relative == CartesianIndex(0, 1)
-        contract_horizontal_norm(inds[1], ket, bra, env)
-    else
-        error("Only contractions for nearest neighbor bonds are implemented.")
-    end
-end
-
-function contract_vertical_norm(
+function contract_local_norm2x1(
         coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
     row, col = Tuple(coord)
@@ -217,7 +225,7 @@ function contract_vertical_norm(
         M_northwest[DNWb; DNWt]
 end
 
-function contract_horizontal_norm(
+function contract_local_norm1x2(
         coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv
     )
     row, col = Tuple(coord)
diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl
index f6add4e9e..ff734f00a 100644
--- a/src/algorithms/contractions/localoperator.jl
+++ b/src/algorithms/contractions/localoperator.jl
@@ -17,23 +17,24 @@ specified by `inds`. The `peps` is contracted with `O` from above and below, and
 sandwich is surrounded with the appropriate environment tensors.
 """
 function contract_local_operator(
-        inds::NTuple{N, CartesianIndex{2}},
-        O::AbstractTensorMap{T, S, N, N},
-        ket::InfinitePEPS, bra::InfinitePEPS,
-        env::CTMRGEnv,
-    ) where {T, S, N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap,
+        ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv,
+    )
+    static_inds = Tuple(Val.(inds))
     return _contract_local_operator(static_inds, O, (ket, bra), env)
 end
 function contract_local_operator(
-        inds::NTuple{N, Tuple{Int, Int}},
-        O::AbstractTensorMap{T, S, N, N},
-        ket::InfinitePEPS, bra::InfinitePEPS,
-        env::CTMRGEnv,
-    ) where {T, S, N}
+        inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap,
+        ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv,
+    )
     return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env)
 end
 
+Base.@deprecate(
+    contract_local_operator(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv),
+    contract_local_operator(collect(inds), ket, bra, env)
+)
+
 # This implements the contraction of an operator acting on sites `inds`.
 # The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance).
 
@@ -264,14 +265,14 @@ on a rectangular patch based on `inds` but replacing the operator with an identi
 that the PEPS norm is computed. (Note that this is not the physical norm of the state.)
 """
 function contract_local_norm(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
+    static_inds = Tuple(Val.(inds))
     return _contract_local_norm(static_inds, (ket, bra), env)
 end
 function contract_local_norm(
-        inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
     return contract_local_norm(CartesianIndex.(inds), ket, bra, env)
 end
 @generated function _contract_local_norm(
@@ -299,6 +300,11 @@ end
     return macroexpand(@__MODULE__, returnex)
 end
 
+Base.@deprecate(
+    contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv),
+    contract_local_norm(collect(inds), ket, bra, env)
+)
+
 @doc """
 $(SIGNATURES)
 
@@ -310,52 +316,84 @@ specified by `inds`. The result is normalized such that `tr(ρ) = 1`.
 """ reduced_densitymatrix
 
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    static_inds = Val.(inds)
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    )
+    length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env)
+
+    if length(inds) == 2
+        if inds[2] - inds[1] == CartesianIndex(1, 0)
+            return reduced_densitymatrix2x1(inds[1], ket, bra, env)
+        elseif inds[2] - inds[1] == CartesianIndex(0, 1)
+            return reduced_densitymatrix1x2(inds[1], ket, bra, env)
+        end
+    end
+
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (ket, bra), env)
 end
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
+    )
     size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported"))
-    static_inds = Val.(inds)
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (state,), env)
 end
 function reduced_densitymatrix(
-        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
+        inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
+    )
     size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions"))
     size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported"))
-    static_inds = Val.(inds)
+    static_inds = Tuple(Val.(inds))
     return _contract_densitymatrix(static_inds, (ket, bra), env)
 end
-function reduced_densitymatrix(
+reduced_densitymatrix(inds, ket::InfinitePEPS, env::CTMRGEnv) =
+    reduced_densitymatrix(inds, ket, ket, env)
+
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), state, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
+        inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
+    ) where {N},
+    reduced_densitymatrix(collect(inds), ket, bra, env)
+)
+
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), ket, bra, env)
-end
-function reduced_densitymatrix(
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), ket, bra, env)
-end
-function reduced_densitymatrix(inds, ket::InfinitePEPS, env::CTMRGEnv)
-    return reduced_densitymatrix(inds, ket, ket, env)
-end
-function reduced_densitymatrix(
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env)
+)
+Base.@deprecate(
+    reduced_densitymatrix(
         inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env::CTMRGEnv
-    ) where {N}
-    return reduced_densitymatrix(CartesianIndex.(inds), state, env)
-end
+    ) where {N},
+    reduced_densitymatrix(collect(CartesianIndex.(inds)), state, env)
+)
 
 # Special case 1x1 density matrix:
 # Keep contraction order but try to optimize intermediate permutations:
 # EE_SWA is largest object so keep largest legs to the front there
-function reduced_densitymatrix(
-        inds::Tuple{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
+function reduced_densitymatrix1x1(
+        inds::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv
     )
-    row, col = Tuple(inds[1])
+    row, col = Tuple(inds)
 
     # Unpack variables and absorb corners
     A = ket[mod1(row, end), mod1(col, end)]
diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl
index a7b9a9d01..324a8d604 100644
--- a/src/algorithms/time_evolution/gaugefix_su.jl
+++ b/src/algorithms/time_evolution/gaugefix_su.jl
@@ -28,7 +28,7 @@ function _trivial_gates(elt::Type{<:Number}, lattice::Matrix{S}) where {S <: Ele
         h = TensorKit.id(elt, V1 ⊗ V2)
         return [site1, site2] => h
     end
-    return TrotterGates(lattice, vec(gates))
+    return LocalCircuit(lattice, vec(gates))
 end
 
 """
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 5e10b2097..9d6842844 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -108,6 +108,7 @@ function _su_iter!(
     # rotate back
     A = _bond_rotation(A, bond[1], rev; inv = true)
     B = _bond_rotation(B, bond[1], rev; inv = true)
+    rev && (s = transpose(s))
     # remove environment weights
     siteA, siteB = map(sites) do site
         return CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))
@@ -127,12 +128,12 @@ end
 One iteration of simple update
 """
 function su_iter(
-        state::InfiniteState, gates::TrotterGates,
+        state::InfiniteState, circuit::LocalCircuit,
         alg::SimpleUpdate, env::SUWeight
     )
     Nr, Nc, = size(state)
     state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
-    for (sites, gate) in gates.terms
+    for (sites, gate) in circuit.gates
         if length(sites) == 1
             # 1-site gate
             # TODO: special treatment for bipartite state
diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl
index bd1c66e73..84a4882b5 100644
--- a/src/algorithms/time_evolution/trotter_gate.jl
+++ b/src/algorithms/time_evolution/trotter_gate.jl
@@ -1,141 +1,120 @@
 """
-    struct TrotterGates{T <: Vector}
+$(TYPEDEF)
 
-Collection of Trotter evolution gates and MPOs obtained from
-a Hamiltonian containing long-range or multi-site terms.
+Circuit consisting of local gates and MPOs.
 
 ## Fields
 
-- `lattice::Matrix{S}`: The lattice on which the gates acts.
-- `terms::T`: The vector of `sites => gate` pairs, where `sites` is a
-vector of `CartesianIndex`s storing the sites on which the `gate` acts.
+$(TYPEDFIELDS)
 """
-struct TrotterGates{T <: Vector, S}
+struct LocalCircuit{O, S}
+    "lattice of physical spaces on which the gates act"
     lattice::Matrix{S}
-    terms::T
-    # TODO: check physical spaces of terms match `lattice`
+
+    "list of `sites => gate` pairs that make up the circuit"
+    gates::Vector{Pair{Vector{CartesianIndex{2}}, O}}
+
+    LocalCircuit{O, S}(lattice::Matrix{S}) where {O, S} =
+        new{O, S}(lattice, Vector{Pair{Vector{CartesianIndex{2}}, O}}())
 end
 
-"""
-    physicalspace(gates::TrotterGates)
+LocalCircuit{O}(lattice::Matrix{<:ElementarySpace}) where {O} =
+    LocalCircuit{O, eltype(lattice)}(lattice)
+LocalCircuit{O}(lattice, gates::Pair...) where {O} = LocalCircuit{O}(lattice, gates)
 
-Return lattice of physical spaces on which the `TrotterGates` is defined.
-"""
-physicalspace(gates::TrotterGates) = gates.lattice
-Base.size(gates::TrotterGates) = size(physicalspace(gates))
+function LocalCircuit{O}(lattice, terms) where {O}
+    operator = LocalCircuit{O}(lattice)
+    for (inds, term) in terms
+        add_factor!(operator, inds, term)
+    end
+    return operator
+end
 
-const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2}
+# Default to Any for eltype: needs to be abstract anyways so not that much to gain
+LocalCircuit(lattice, terms) = LocalCircuit{Any}(lattice, terms)
+LocalCircuit(lattice, terms::Pair...) = LocalCircuit(lattice, terms)
+
+add_factor!(operator::LocalCircuit, inds::Tuple, term::AbstractTensorMap) = add_factor!(operator, collect(inds), term)
+add_factor!(operator::LocalCircuit, inds::Vector, term::AbstractTensorMap) = add_factor!(operator, map(CartesianIndex{2}, inds), term)
+function add_factor!(operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap)
+    # input checks
+    length(inds) == numin(term) == numout(term) || throw(ArgumentError("Incompatible number of indices and tensor legs"))
+    for (i, ind) in enumerate(inds)
+        ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator)))
+        physicalspace(operator, ind_translated) == domain(term)[i] == codomain(term)[i] ||
+            throw(SpaceMismatch("Incompatible physical spaces at $(ind)."))
+    end
 
-"""
-    is_equivalent_site(
-        site1::CartesianIndex{2}, site2::CartesianIndex{2},
-        (Nrow, Ncol)::NTuple{2, Int}
-    )
+    # permute input
+    if !issorted(inds)
+        I = sortperm(inds)
+        inds = inds[I]
+        term = permute(term, (Tuple(I), Tuple(I) .+ numout(term)))
+    end
 
-Check if two lattice sites are related by a (periodic) lattice translation.
-"""
-function is_equivalent_site(
-        site1::CartesianIndex{2}, site2::CartesianIndex{2},
-        (Nrow, Ncol)::NTuple{2, Int}
-    )
-    shift = site1 - site2
-    return mod(shift[1], Nrow) == 0 && mod(shift[2], Ncol) == 0
-end
+    # translate coordinates
+    I1 = first(inds)
+    I1_mod = CartesianIndex(mod1.(Tuple(I1), size(operator)))
+    inds .-= (I1 - I1_mod)
 
-"""
-    _get_site_term(ham::LocalOperator, site::CartesianIndex{2})
+    push!(operator.gates, inds => term)
 
-Get the sum of all 1-site terms at `site` in `ham`.
-If there are no such terms, return the zero operator at `site`.
-"""
-function _get_site_term(ham::LocalOperator, site::CartesianIndex{2})
-    r, c = mod1.(Tuple(site), size(ham))
-    V = physicalspace(ham)[r, c]
-    term = zeros(scalartype(ham), V ← V)
-    for (sites, op) in ham.terms
-        length(sites) != 1 && continue
-        if is_equivalent_site(sites[1], site, size(ham))
-            term = term + op
+    return operator
+end
+# for MPO term
+# TODO: consider directly use MPSKit.FiniteMPO
+function add_factor!(
+        operator::LocalCircuit, inds::Vector{CartesianIndex{2}}, term::Vector{M}
+    ) where {M <: AbstractTensorMap}
+    # input checks
+    length(inds) >= 2 || throw(ArgumentError("Gate MPO must act on 2 or more sites."))
+    length(inds) == length(term) || throw(ArgumentError("Incompatible number of indices and length of gate MPO."))
+    allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates."))
+    for (i, (ind, t)) in enumerate(zip(inds, term))
+        ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator)))
+        out_ax = (i == 1) ? 1 : 2
+        in_ax = (i == 1) ? 2 : 3
+        physicalspace(operator, ind_translated) == space(t, out_ax) == space(t, in_ax)' ||
+            throw(SpaceMismatch("Incompatible physical spaces at $(ind)."))
+        if i >= 2
+            ind_prev = inds[i - 1]
+            sum(Tuple(ind - ind_prev) .^ 2) == 1 || throw(ArgumentError("Two consecutive sites in `inds` must be nearest neighbours for MPO terms."))
         end
     end
-    return term
+    # for MPO term, `inds` should not be sorted
+    push!(operator.gates, inds => term)
+    return operator
 end
 
-"""
-    is_equivalent_bond(
-        bond1::NTuple{2, CartesianIndex{2}}, bond2::NTuple{2, CartesianIndex{2}},
-        (Nrow, Ncol)::NTuple{2, Int},
-    )
-
-Check if two 2-site bonds are related by a (periodic) lattice translation.
-"""
-function is_equivalent_bond(
-        bond1::NTuple{2, CartesianIndex{2}}, bond2::NTuple{2, CartesianIndex{2}},
-        (Nrow, Ncol)::NTuple{2, Int},
-    )
-    r1 = bond1[1] - bond1[2]
-    r2 = bond2[1] - bond2[2]
-    shift_row = bond1[1][1] - bond2[1][1]
-    shift_col = bond1[1][2] - bond2[1][2]
-    return r1 == r2 && mod(shift_row, Nrow) == 0 && mod(shift_col, Ncol) == 0
+function checklattice(::Type{Bool}, H1::LocalCircuit, H2::LocalCircuit)
+    return physicalspace(H1) == physicalspace(H2)
+end
+function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalCircuit)
+    return physicalspace(peps) == physicalspace(O)
+end
+function checklattice(::Type{Bool}, H::LocalCircuit, peps::InfinitePEPS)
+    return checklattice(Bool, peps, H)
+end
+function checklattice(::Type{Bool}, pepo::InfinitePEPO, O::LocalCircuit)
+    return size(pepo, 3) == 1 && physicalspace(pepo) == physicalspace(O)
+end
+function checklattice(::Type{Bool}, O::LocalCircuit, pepo::InfinitePEPO)
+    return checklattice(Bool, pepo, O)
 end
 
 """
-    _get_bond_term(ham::LocalOperator, bond::NTuple{2, CartesianIndex{2}})
+    physicalspace(gates::LocalCircuit)
 
-Get the sum of all 2-site terms on `bond` in `ham`.
-If there are no such terms, return the zero operator on `bond`.
+Return lattice of physical spaces on which the `LocalCircuit` is defined.
 """
-function _get_bond_term(ham::LocalOperator, bond::NTuple{2, CartesianIndex{2}})
-    # create zero operator
-    r1, c1 = mod1.(Tuple(bond[1]), size(ham))
-    r2, c2 = mod1.(Tuple(bond[2]), size(ham))
-    V1 = physicalspace(ham)[r1, c1]
-    V2 = physicalspace(ham)[r2, c2]
-    term = zeros(scalartype(ham), V1 ⊗ V2 ← V1 ⊗ V2)
-    for (sites, op) in ham.terms
-        length(sites) != 2 && continue
-        if is_equivalent_bond(sites, bond, size(ham))
-            term += op
-        elseif is_equivalent_bond(sites, reverse(bond), size(ham))
-            op′ = permute(op, ((2, 1), (4, 3)); copy = true)
-            term += op′
-        end
-    end
-    return term
-end
+physicalspace(gates::LocalCircuit) = gates.lattice
+physicalspace(gates::LocalCircuit, args...) = physicalspace(gates)[args...]
+Base.size(gates::LocalCircuit) = size(physicalspace(gates))
 
-"""
-Get coordinates of sites in the 3-site triangular cluster
-used in Trotter evolution with next-nearest neighbor gates,
-with southwest corner at `[row, col]`.
-```
-    NORTHWEST   NORTHEAST
-        2---1   3---2
-        |           |
-        3           1
-
-        1           3
-        |           |
-        2---3   1---2
-    SOUTHWEST   SOUTHEAST
-```
-"""
-function _nnn_cluster_sites(dir::Int, row::Int, col::Int)
-    @assert 1 <= dir <= 4
-    return if dir == NORTHWEST
-        map(CartesianIndex, [(row - 1, col + 1), (row - 1, col), (row, col)])
-    elseif dir == NORTHEAST
-        map(CartesianIndex, [(row, col + 1), (row - 1, col + 1), (row - 1, col)])
-    elseif dir == SOUTHEAST
-        map(CartesianIndex, [(row, col), (row, col + 1), (row - 1, col + 1)])
-    else # dir == SOUTHWEST
-        map(CartesianIndex, [(row - 1, col), (row, col), (row, col + 1)])
-    end
-end
+const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2}
 
 """
-Convert an N-site gate to MPO form by SVD, 
+Convert an N-site gate (N ≥ 2) to MPO by SVD, 
 in which the axes are ordered as
 ```
     site 1      mid sites      site N
@@ -147,8 +126,10 @@ in which the axes are ordered as
 ```
 """
 function gate_to_mpo(
-        gate::AbstractTensorMap{T, S, N, N}, trunc = trunctol(; atol = MPSKit.Defaults.tol)
-    ) where {T <: Number, S <: ElementarySpace, N}
+        gate::AbstractTensorMap{<:Any, <:Any, N, N};
+        trunc = trunctol(; atol = MPSKit.Defaults.tol)
+    ) where {N}
+    @assert N >= 2
     Os = MPSKit.decompose_localmpo(MPSKit.add_util_leg(gate), trunc)
     return map(1:N) do i
         if i == 1
@@ -192,11 +173,12 @@ end
 """
 Trotterize a trivial Hamiltonian `H` containing only 1-site terms.
 """
-function _trotterize_1site!(gates::Vector, H::LocalOperator, dt::Number; atol::Real)
-    for site in CartesianIndices(size(H))
-        gate = _get_site_term(H, site)
-        (norm(gate) <= atol) && continue
-        push!(gates, [site] => exp(-dt * gate))
+function _trotterize_1site!(gates::Vector, H::LocalOperator, dt::Number)
+    for x in CartesianIndices(size(H))
+        coord = [x]
+        haskey(H.terms, coord) || continue
+        gate = exp(H.terms[coord] * -dt)
+        push!(gates, coord => gate)
     end
     return gates
 end
@@ -204,51 +186,59 @@ end
 """
 Trotterize nearest neighbor terms (grouped with 1-site terms)
 in the Hamiltonian `H`.
-
-Gate order: `(d, c, r)`
-- d = 1: horizontal bond ((r, c), (r, c+1))
-- d = 2:   vertical bond ((r, c), (r-1, c))
 """
-function _trotterize_nn2site!(gates::Vector, H::LocalOperator, dt::Number; atol::Real, force_mpo::Bool = false)
-    Nr, Nc = size(H)
-    T = scalartype(H)
-    for (d, c, r) in Iterators.product(1:2, 1:Nc, 1:Nr)
-        site1 = CartesianIndex(r, c)
-        site2 = (d == 1) ? CartesianIndex(r, c + 1) : CartesianIndex(r - 1, c)
-        # group with 1-site terms
-        s1term = _get_site_term(H, site1)
-        unit1 = TensorKit.id(T, space(s1term, 1))
-        s2term = _get_site_term(H, site2)
-        unit2 = TensorKit.id(T, space(s2term, 1))
-        gate = _get_bond_term(H, (site1, site2))
-        gate = gate + (s1term ⊗ unit2 + unit1 ⊗ s2term) / 4
-        (norm(gate) <= atol) && continue
-        gate = exp(-dt * gate)
+function _trotterize_nn2site!(
+        gates::Vector, H::LocalOperator, dt::Number; force_mpo::Bool = false
+    )
+    vs = [CartesianIndex(0, 1), CartesianIndex(1, 0)]
+    for x in CartesianIndices(size(H)), v in vs
+        y = x + v
+        coord = [x, y]
+        haskey(H.terms, coord) || continue
+        gate = exp(H.terms[coord] * -dt)
         force_mpo && (gate = gate_to_mpo(gate))
-        push!(gates, [site1, site2] => gate)
+        push!(gates, coord => gate)
     end
     return gates
 end
 
 """
-Trotterize a next-nearest neighbor terms in a Hamiltonian.
+Trotterize next-nearest neighbor terms in a Hamiltonian,
+converting them to 3-site MPO gates. 
+For each gate, the order of sites is
+```
+    2---3   1---2
+    |           |
+    1           3
 
-Gate order: `(c, r, d)`
-- d = 1 (NORTHWEST), ..., 4 (SOUTHWEST) labels the triangular 3-site clusters.
+    1           3
+    |           |
+    2---3   1---2
+```
 """
-function _trotterize_nnn2site!(gates::Vector, H::LocalOperator, dt::Number; atol::Real)
-    Nr, Nc = size(H)
+function _trotterize_nnn2site!(gates::Vector, H::LocalOperator, dt::Number)
     T = scalartype(H)
-    for (c, r, d) in Iterators.product(1:Nc, 1:Nr, 1:4)
-        sites = _nnn_cluster_sites(d, r, c)
-        gate = _get_bond_term(H, (sites[1], sites[3]))
-        (norm(gate) <= atol) && continue
-        gate = exp(-(dt / 2) * gate) # account for double counting
-        # combine with identity at sites[2]
-        r2, c2 = mod1(sites[2][1], Nr), mod1(sites[2][2], Nc)
-        id_ = TensorKit.id(T, physicalspace(H)[r2, c2])
-        gate = permute(gate ⊗ id_, ((1, 3, 2), (4, 6, 5)))
-        push!(gates, sites => gate_to_mpo(gate))
+    vs = [
+        # ⌞ next-nearest-neighbour
+        (CartesianIndex(1, 0), CartesianIndex(1, 1)),
+        # ⌜ next-nearest-neighbour
+        (CartesianIndex(-1, 0), CartesianIndex(-1, 1)),
+        # ⌝ next-nearest-neighbour
+        (CartesianIndex(0, 1), CartesianIndex(1, 1)),
+        # ⌟ next-nearest-neighbour
+        (CartesianIndex(0, 1), CartesianIndex(-1, 1)),
+    ]
+    for x1 in CartesianIndices(size(H)), v in vs
+        x2, x3 = x1 + v[1], x1 + v[2]
+        coord = [x1, x3]
+        haskey(H.terms, coord) || continue
+        gate = gate_to_mpo(exp(H.terms[coord] * -dt / 2))
+        x2′ = CartesianIndex(mod1.(Tuple(x2), size(H)))
+        b = TensorKit.BraidingTensor{T}(
+            physicalspace(H, x2′), left_virtualspace(gate[2])
+        )
+        insert!(gate, 2, TensorMap(b))
+        push!(gates, [x1, x2, x3] => gate)
     end
     return gates
 end
@@ -271,21 +261,16 @@ function trotterize(
         symmetrize_gates::Bool = false, force_mpo::Bool = false
     )
     dist = _check_hamiltonian_for_trotter(H)
-    T = scalartype(H)
-    atol = eps(real(T))^(3 / 4)
+
     dt′ = symmetrize_gates ? (dt / 2) : dt
-    gates = Vector{Pair{Any, Any}}()
 
-    # TODO: order of gates is fixed for more tight control.
-    # Consider directly iterating over H.terms in the future.
+    gates = Vector{Pair{Vector{CartesianIndex{2}}, Any}}()
+
+    dist >= 0 && _trotterize_1site!(gates, H, dt′)
+    dist >= 1 && _trotterize_nn2site!(gates, H, dt′; force_mpo)
+    dist >= 2 && _trotterize_nnn2site!(gates, H, dt′)
 
-    # 1-site gates are only constructed when H only has 1-site terms
-    dist == 0 && _trotterize_1site!(gates, H, dt′; atol)
-    # 2-site NN gates grouped with 1-site terms
-    dist >= 1 && _trotterize_nn2site!(gates, H, dt′; atol, force_mpo)
-    # 3-site NNN gate MPOs
-    dist >= 2 && _trotterize_nnn2site!(gates, H, dt′; atol)
+    symmetrize_gates && append!(gates, reverse(gates))
 
-    symmetrize_gates && push!(gates, reverse(gates)...)
-    return TrotterGates(physicalspace(H), gates)
+    return LocalCircuit(physicalspace(H), gates)
 end
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 13d4652ae..647ad1221 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -12,7 +12,7 @@ function MPSKit.expectation_value(
         ket::Union{InfinitePEPS, InfinitePEPO}, env::CTMRGEnv
     )
     checklattice(bra, O, ket)
-    term_vals = dtmap([O.terms...]) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
+    term_vals = dtmap(collect(O.terms)) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
         ρ = reduced_densitymatrix(inds, ket, bra, env)
         return trmul(operator, ρ)
     end
@@ -23,7 +23,7 @@ function MPSKit.expectation_value(
         state::InfinitePEPO, O::LocalOperator, env::CTMRGEnv
     )
     checklattice(state, O)
-    term_vals = dtmap([O.terms...]) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
+    term_vals = dtmap(collect(O.terms)) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
         ρ = reduced_densitymatrix(inds, state, env)
         return trmul(operator, ρ)
     end
diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index 44ea66b49..251d0c5a5 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -3,67 +3,93 @@
 """
 $(TYPEDEF)
 
-A sum of local operators acting on a lattice. The lattice is stored as a matrix of vector spaces,
-and the terms are stored as a tuple of pairs of indices and operators.
+A sum of local operators acting on a lattice.
+The lattice is stored as a matrix of vector spaces, and the terms are stored as a `Dict` of indices mapping to operators.
 
 ## Fields
 
+$(TYPEDFIELDS)
 - `lattice::Matrix{S}`: The lattice on which the operator acts.
-- `terms::T`: The terms of the operator, stored as a tuple of pairs of indices and operators.
+- `terms::Dict{Vector{CartesianIndex{2}}, O}`: The terms of the operator, mapping coordinates to operators
 
 ## Constructors
 
     LocalOperator(lattice::Matrix{S}, terms::Pair...)
-    LocalOperator{T,S}(lattice::Matrix{S}, terms::T) where {T,S}
+    LocalOperator{T, S}(lattice::Matrix{S}, terms::T) where {T,S}
 
 ## Examples
 
 ```julia
 lattice = fill(ℂ^2, 1, 1) # single-site unitcell
-O1 = LocalOperator(lattice, ((1, 1),) => σx, ((1, 1), (1, 2)) => σx ⊗ σx, ((1, 1), (2, 1)) => σx ⊗ σx)
+O1 = LocalOperator(lattice, [(1, 1),] => σx, [(1, 1), (1, 2)] => σx ⊗ σx, [(1, 1), (2, 1)] => σx ⊗ σx)
 ```
 """
-struct LocalOperator{T <: Tuple, S}
+struct LocalOperator{O, S}
+    "lattice of physical spaces on which the gates act"
     lattice::Matrix{S}
-    terms::T
-    function LocalOperator{T, S}(lattice::Matrix{S}, terms::T) where {T, S}
-        plattice = PeriodicArray(lattice)
-        # Check if the indices of the operator are valid with themselves and the lattice
-        for (inds, operator) in terms
-            @assert operator isa AbstractTensorMap
-            @assert eltype(inds) <: CartesianIndex
-            @assert numout(operator) == numin(operator) == length(inds)
-            @assert spacetype(operator) == S
-
-            for i in 1:length(inds)
-                @assert space(operator, i) == plattice[inds[i]]
-            end
-        end
-        return new{T, S}(lattice, terms)
+
+    "list of `sites => term` pairs that make up the operator"
+    terms::Dict{Vector{CartesianIndex{2}}, O}
+
+    LocalOperator{O, S}(lattice::Matrix{S}) where {O, S} =
+        new{O, S}(lattice, Dict{Vector{CartesianIndex{2}}, O}())
+end
+
+LocalOperator{O}(lattice::Matrix{<:ElementarySpace}) where {O} =
+    LocalOperator{O, eltype(lattice)}(lattice)
+LocalOperator{O}(lattice, terms::Pair...) where {O} = LocalOperator{O}(lattice, terms)
+
+function LocalOperator{O}(lattice, terms) where {O}
+    operator = LocalOperator{O}(lattice)
+    for (inds, term) in terms
+        add_term!(operator, inds, term)
     end
+    return operator
 end
-function LocalOperator(
-        lattice::Matrix, terms::Pair...;
-        atol = maximum(x -> eps(real(scalartype(x[2])))^(3 / 4), terms),
+
+# Default to Any for eltype: needs to be abstract anyways so not that much to gain
+LocalOperator(lattice, terms) = LocalOperator{Any}(lattice, terms)
+LocalOperator(lattice, terms::Pair...) = LocalOperator(lattice, terms)
+# TODO: add terms beyond AbstractTensorMap
+# e.g. tensor product of 1-site operators, MPOs
+add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap) = add_term!(operator, collect(inds), term)
+add_term!(operator::LocalOperator, inds::Vector, term::AbstractTensorMap) = add_term!(operator, map(CartesianIndex{2}, inds), term)
+function add_term!(
+        operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap;
+        atol = zero(real(scalartype(term))),
     )
-    allinds = getindex.(terms, 1)
-    alloperators = getindex.(terms, 2)
-
-    relevant_terms = []
-    for inds in unique(allinds)
-        operator = sum(alloperators[findall(==(inds), allinds)])
-        cinds = if !(eltype(inds) <: CartesianIndex) # force indices to be CartesianIndices
-            map(CartesianIndex, inds)
-        else
-            inds
-        end
-        norm(operator) > atol && push!(relevant_terms, cinds => operator)
+    # input checks
+    length(inds) == numin(term) == numout(term) || throw(ArgumentError("Incompatible number of indices and tensor legs"))
+    allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates."))
+    for (i, ind) in enumerate(inds)
+        ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator)))
+        physicalspace(operator, ind_translated) == domain(term)[i] == codomain(term)[i] ||
+            throw(SpaceMismatch("Incompatible physical spaces"))
+    end
+    norm(term) <= atol && return operator # skip adding negligible terms
+
+    # permute input
+    if !issorted(inds)
+        I = sortperm(inds)
+        inds = inds[I]
+        term = permute(term, (Tuple(I), Tuple(I) .+ numout(term)))
+    end
+
+    # translate coordinates
+    I1 = first(inds)
+    I1_mod = CartesianIndex(mod1.(Tuple(I1), size(operator)))
+    inds .-= (I1 - I1_mod)
+
+    if haskey(operator.terms, inds)
+        operator.terms[inds] = VI.add!!(operator.terms[inds], term)
+    else
+        operator.terms[inds] = term
     end
 
-    terms_tuple = Tuple(relevant_terms)
-    return LocalOperator{typeof(terms_tuple), eltype(lattice)}(lattice, terms_tuple)
+    return operator
 end
 
+
 """
     checklattice(Bool, args...)
     checklattice(args...)
@@ -95,14 +121,24 @@ function checklattice(::Type{Bool}, O::LocalOperator, pepo::InfinitePEPO)
 end
 @non_differentiable checklattice(args...)
 
-function Base.repeat(O::LocalOperator, m::Int, n::Int)
-    lattice = repeat(O.lattice, m, n)
-    terms = []
-    for (inds, operator) in O.terms, i in 1:m, j in 1:n
-        offset = CartesianIndex((i - 1) * size(O.lattice, 1), (j - 1) * size(O.lattice, 2))
-        push!(terms, (inds .+ Ref(offset)) => operator)
+function Base.similar(operator::LocalOperator, lattice::Matrix{<:ElementarySpace})
+    return similar(operator, eltype(operator), lattice)
+end
+function Base.similar(
+        operator::LocalOperator, ::Type{O} = eltype(operator), lattice::Matrix{<:ElementarySpace} = physicalspace(operator)
+    ) where {O}
+    return LocalOperator{O}(lattice)
+end
+
+function Base.repeat(operator::LocalOperator, m::Int, n::Int)
+    operator_repeated = similar(operator, repeat(physicalspace(operator), m, n))
+    for i in 1:m, j in 1:n
+        offset = CartesianIndex((i - 1) * size(operator, 1), (j - 1) * size(operator, 2))
+        for (inds, term) in operator.terms
+            add_term!(operator_repeated, inds .+ offset, term)
+        end
     end
-    return LocalOperator(lattice, terms...)
+    return operator_repeated
 end
 
 """
@@ -110,11 +146,11 @@ end
 
 Return lattice of physical spaces on which the `LocalOperator` is defined.
 """
-function physicalspace(O::LocalOperator)
-    return O.lattice
-end
+physicalspace(O::LocalOperator) = O.lattice
+physicalspace(O::LocalOperator, args...) = physicalspace(O)[args...]
 
-Base.size(O::LocalOperator) = size(physicalspace(O))
+Base.size(O::LocalOperator, args...) = size(physicalspace(O), args...)
+Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O
 
 # Real and imaginary part
 # -----------------------
@@ -127,10 +163,8 @@ end
 
 # Linear Algebra
 # --------------
-function Base.:*(α::Number, O::LocalOperator)
-    scaled_terms = map(((inds, operator),) -> (inds => α * operator), O.terms)
-    return LocalOperator{typeof(scaled_terms), eltype(O.lattice)}(O.lattice, scaled_terms)
-end
+Base.:*(α::Number, O::LocalOperator) =
+    LocalOperator(physicalspace(O), inds => α * operator for (inds, operator) in O.terms)
 Base.:*(O::LocalOperator, α::Number) = α * O
 
 Base.:/(O::LocalOperator, α::Number) = O * inv(α)
@@ -138,7 +172,7 @@ Base.:\(α::Number, O::LocalOperator) = inv(α) * O
 
 function Base.:+(O1::LocalOperator, O2::LocalOperator)
     checklattice(O1, O2)
-    return LocalOperator(O1.lattice, O1.terms..., O2.terms...)
+    return LocalOperator(physicalspace(O1), mergewith(VI.add, O1.terms, O2.terms))
 end
 
 Base.:-(O::LocalOperator) = -1 * O
@@ -147,20 +181,17 @@ Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2)
 # VectorInterface
 # ---------------
 
-function VI.scalartype(::Type{<:LocalOperator{T}}) where {T}
-    return promote_type((scalartype(last(fieldtypes(p))) for p in fieldtypes(T))...)
+# Since we allow abstract types in T, value and type domain might not match
+function VI.scalartype(operator::LocalOperator)
+    return promote_type((scalartype(term[2]) for term in operator.terms)...)
 end
 
+
 # Equivalence
 # -----------
 
-function Base.:(==)(O₁::LocalOperator, O₂::LocalOperator)
-    lat = O₁.lattice == O₂.lattice
-    terms = all(zip(O₁.terms, O₂.terms)) do (t₁, t₂)
-        return t₁ == t₂
-    end
-    return lat && terms
-end
+Base.:(==)(O₁::LocalOperator, O₂::LocalOperator) =
+    physicalspace(O₁) == physicalspace(O₂) && O₁.terms == O₂.terms
 
 # Rotation
 # ----------------------
@@ -178,27 +209,27 @@ function siterot180(site::CartesianIndex{2}, unitcell::NTuple{2, Int})
 end
 
 function Base.rotr90(H::LocalOperator)
-    Hsize = size(H.lattice)
-    lattice2 = rotr90(H.lattice)
-    terms2 = ((Tuple(siterotr90(site, Hsize) for site in sites) => op) for (sites, op) in H.terms)
-    return LocalOperator(lattice2, terms2...)
+    Hsize = size(H)
+    lattice2 = rotr90(physicalspace(H))
+    terms2 = (siterotr90.(inds, Ref(Hsize)) => term for (inds, term) in H.terms)
+    return LocalOperator(lattice2, terms2)
 end
 function Base.rotl90(H::LocalOperator)
-    Hsize = size(H.lattice)
-    lattice2 = rotl90(H.lattice)
-    terms2 = ((Tuple(siterotl90(site, Hsize) for site in sites) => op) for (sites, op) in H.terms)
-    return LocalOperator(lattice2, terms2...)
+    Hsize = size(H)
+    lattice2 = rotl90(physicalspace(H))
+    terms2 = (siterotl90.(inds, Ref(Hsize)) => term for (inds, term) in H.terms)
+    return LocalOperator(lattice2, terms2)
 end
 function Base.rot180(H::LocalOperator)
-    Hsize = size(H.lattice)
-    lattice2 = rot180(H.lattice)
-    terms2 = ((Tuple(siterot180(site, Hsize) for site in sites) => op) for (sites, op) in H.terms)
-    return LocalOperator(lattice2, terms2...)
+    Hsize = size(H)
+    lattice2 = rot180(physicalspace(H))
+    terms2 = (siterot180.(inds, Ref(Hsize)) => term for (inds, term) in H.terms)
+    return LocalOperator(lattice2, terms2)
 end
 
 # Charge shifting
 # ---------------
-TensorKit.spacetype(::Type{T}) where {S, T <: LocalOperator{<:Any, S}} = S
+TensorKit.spacetype(::Type{<:LocalOperator{<:Any, S}}) where {S} = S
 
 @generated function _fuse_isomorphisms(
         op::AbstractTensorMap{<:Any, S, N, N}, fs::Vector{<:AbstractTensorMap{<:Any, S, 1, 2}}
@@ -244,7 +275,7 @@ Change the spaces of a `LocalOperator` by fusing in an auxiliary charge into the
 the operator on every site, according to a given matrix of 'auxiliary' physical charges.
 """
 function MPSKit.add_physical_charge(H::LocalOperator, charges::AbstractMatrix{<:Sector})
-    size(physicalspace(H)) == size(charges) ||
+    size(H) == size(charges) ||
         throw(ArgumentError("Incompatible lattice and auxiliary charge sizes"))
     sectortype(H) === eltype(charges) ||
         throw(SectorMismatch("Incompatible lattice and auxiliary charge sizes"))
@@ -252,16 +283,13 @@ function MPSKit.add_physical_charge(H::LocalOperator, charges::AbstractMatrix{<:
     # auxiliary spaces will be fused into codomain, so need to dualize the space to fuse
     # the charge into the domain as desired
     # also, make indexing periodic for convenience
-    Paux = PeriodicArray(map(c -> Vect[typeof(c)](c => 1)', charges))
+    Paux = PeriodicArray(map(c -> spacetype(H)(c => 1)', charges))
 
     # new physical spaces
     Pspaces = map(fuse, physicalspace(H), Paux)
 
-    new_terms = map(H.terms) do (sites, op)
-        Paux_slice = map(Base.Fix1(getindex, Paux), sites)
-        return sites => _fuse_ids(op, Paux_slice)
-    end
-    H´ = LocalOperator(Pspaces, new_terms...)
-
-    return H´
+    return LocalOperator(
+        Pspaces,
+        inds => _fuse_ids(op, Tuple(map(Base.Fix1(getindex, Paux), inds))) for (inds, op) in H.terms
+    )
 end
diff --git a/src/operators/models.jl b/src/operators/models.jl
index 37eceaa72..6e1d927b9 100644
--- a/src/operators/models.jl
+++ b/src/operators/models.jl
@@ -15,8 +15,8 @@ function nearest_neighbour_hamiltonian(
     for I in eachindex(IndexCartesian(), lattice)
         J1 = I + CartesianIndex(1, 0)
         J2 = I + CartesianIndex(0, 1)
-        push!(terms, (I, J1) => h)
-        push!(terms, (I, J2) => h)
+        push!(terms, [I, J1] => h)
+        push!(terms, [I, J2] => h)
     end
     return LocalOperator(lattice, terms...)
 end
@@ -35,7 +35,7 @@ function MPSKitModels.transverse_field_ising(
     return LocalOperator(
         spaces,
         (neighbor => ZZ for neighbor in nearest_neighbours(lattice))...,
-        ((idx,) => X for idx in vertices(lattice))...,
+        ([idx] => X for idx in vertices(lattice))...,
     )
 end
 
@@ -104,7 +104,7 @@ function MPSKitModels.bose_hubbard_model(
     H = LocalOperator(
         spaces,
         (neighbor => -t * hopping_term for neighbor in nearest_neighbours(lattice))...,
-        ((idx,) => U / 2 * interaction_term - mu * N for idx in vertices(lattice))...,
+        ([idx] => U / 2 * interaction_term - mu * N for idx in vertices(lattice))...,
     )
 
     if symmetry === Trivial
@@ -223,7 +223,7 @@ function pwave_superconductor(
     y_neighbors = filter(n -> n[2].I[1] > n[1].I[1], nearest_neighbours(lattice))
     return LocalOperator(
         spaces,
-        ((idx,) => h0 for idx in vertices(lattice))...,
+        ([idx] => h0 for idx in vertices(lattice))...,
         (neighbor => hx for neighbor in x_neighbors)...,
         (neighbor => hy for neighbor in y_neighbors)...,
     )
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 6abb4146b..ad258ca78 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -31,7 +31,7 @@ function heisenberg_XYZ_c4v(
         rmul!(S_zz(T, S; spin = spin), Jz)
     spaces = fill(domain(term)[1], (1, 1))
     return LocalOperator( # horizontal and vertical contributions are identical
-        spaces, (CartesianIndex(1, 1), CartesianIndex(1, 2)) => 2 * term
+        spaces, [CartesianIndex(1, 1), CartesianIndex(1, 2)] => 2 * term
     )
 end
 
@@ -45,15 +45,17 @@ end
     # optimize energy and compute correlation lengths
     peps, env, E, = fixedpoint(H, peps₀, env₀; optimizer_alg = (; tol = gradtol, maxiter = 25))
     ξ_h, ξ_v, = correlation_length(peps, env)
-
+    @info "Optimized energy = $E."
     @test E ≈ E_ref atol = 1.0e-2
+    @info "ξₕ = $(ξ_h)."
+    @info "ξᵥ = $(ξ_v)."
     @test all(@. ξ_h > 0 && ξ_v > 0)
 end
 
 @testset "C4v AD optimization with scalartype T=$T and projector_alg=$projector_alg" for (T, projector_alg) in
     Iterators.product([Float64, ComplexF64], [:c4v_eigh, :c4v_qr])
     # initialize symmetric states
-    Random.seed!(12345)
+    Random.seed!(1234567)
     symm = RotateReflect()
     H′ = heisenberg_XYZ_c4v(InfiniteSquare())
     H = T <: Real ? real(H′) : H′
@@ -69,8 +71,10 @@ end
         boundary_alg = (; alg = :c4v, projector_alg, maxiter = 500),
     )
     ξ_h, ξ_v, = correlation_length(peps, env)
-
+    @info "Optimized energy = $E."
     @test E ≈ E_ref atol = 1.0e-2
+    @info "ξₕ = $(ξ_h)."
+    @info "ξᵥ = $(ξ_v)."
     @test only(ξ_h) ≈ only(ξ_v)
 end
 
@@ -85,8 +89,10 @@ end
     # optimize energy and compute correlation lengths
     peps, env, E, = fixedpoint(H, peps₀, env₀; optimizer_alg = (; tol = gradtol, maxiter = 25))
     ξ_h, ξ_v, = correlation_length(peps, env)
-
+    @info "Optimized energy = $E."
     @test E ≈ 2 * E_ref atol = 1.0e-2
+    @info "ξₕ = $(ξ_h)."
+    @info "ξᵥ = $(ξ_v)."
     @test all(@. ξ_h > 0 && ξ_v > 0)
 end
 
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 591e11952..53fad0a68 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -91,7 +91,7 @@ end
         @test mpo_to_gate3(gs) ≈ gate
         for gate_ax in 1:2
             Ms2 = _flip_virtuals!(deepcopy(Ms1), flips)
-            PEPSKit._apply_gatempo!(Ms2, gs)
+            PEPSKit._apply_gatempo!(Ms2, gs; gate_ax)
             fid = fidelity_cluster(
                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms1],
                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms2]
@@ -120,17 +120,12 @@ end
     # applying 2-site gates decomposed to MPO or not,
     # resulting energy should be almost the same
     e_sites = map((true, false)) do force_mpo
-        dts = [1.0e-2, 1.0e-2]
-        tols = [1.0e-6, 1.0e-8]
         peps, wts = deepcopy(peps0), deepcopy(wts0)
-        for (n, (dt, tol)) in enumerate(zip(dts, tols))
-            trunc = truncerror(; atol = 1.0e-10) & truncrank(n == 1 ? 4 : 2)
-            alg = SimpleUpdate(; trunc, force_mpo)
-            peps, wts, = time_evolve(
-                peps, ham, dt, 10000, alg, wts;
-                tol, symmetrize_gates = true, check_interval = 1000
-            )
-        end
+        trunc = truncerror(; atol = 1.0e-10) & truncrank(4)
+        alg = SimpleUpdate(; trunc, force_mpo)
+        peps, wts, = time_evolve(
+            peps, ham, 0.01, 10000, alg, wts; tol = 1.0e-6, check_interval = 1000
+        )
         normalize!.(peps.A, Inf)
         env = CTMRGEnv(wts)
         for trunc in truncs_env
diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl
index 2b3385156..30067d3d1 100644
--- a/test/toolbox/densitymatrices.jl
+++ b/test/toolbox/densitymatrices.jl
@@ -57,8 +57,8 @@ end
     O_doubled_singlesite = LocalOperator(physicalspace(ρ_peps), (site,) => O_doubled)
     E2 = expectation_value(ρ_peps, O_doubled_singlesite, ρ_peps, env)
     @test E1 ≈ E2
-    val = contract_local_operator((site,), O_doubled, ρ_peps, ρ_peps, env)
-    nrm = contract_local_norm((site,), ρ_peps, ρ_peps, env)
+    val = contract_local_operator([site], O_doubled, ρ_peps, ρ_peps, env)
+    nrm = contract_local_norm([site], ρ_peps, ρ_peps, env)
     @test E1 ≈ val / nrm
 
     # two sites
@@ -71,8 +71,8 @@ end
         O_doubled_twosite = LocalOperator(physicalspace(ρ_peps), inds => O_doubled ⊗ O_doubled)
         E2 = expectation_value(ρ_peps, O_doubled_twosite, ρ_peps, env)
         @test E1 ≈ E2
-        val = contract_local_operator(inds, O_doubled ⊗ O_doubled, ρ_peps, ρ_peps, env)
-        nrm = contract_local_norm(inds, ρ_peps, ρ_peps, env)
+        val = contract_local_operator(collect(inds), O_doubled ⊗ O_doubled, ρ_peps, ρ_peps, env)
+        nrm = contract_local_norm(collect(inds), ρ_peps, ρ_peps, env)
         @test E1 ≈ val / nrm
     end
 end