diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index c152c578a..3fb28fa29 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -390,3 +390,53 @@ function Base.isapprox(
     # don't take square roots to avoid precision loss
     return norm₁₂² ≤ max(atol^2, rtol^2 * max(norm₁², norm₂²))
 end
+
+@doc """
+    swap(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme)
+    swap!(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=SDD(), trscheme)
+
+Compose the mpo with a swap gate applied to indices `i` and `i + 1`, effectively creating an
+operator that acts on the Hilbert spaces with those factors swapped.
+The keyword arguments `alg` and `trscheme` can be used to control how the resulting tensor
+is truncated again.
+""" swap, swap!
+
+swap(mpo::FiniteMPO, i::Integer; kwargs...) = swap!(copy(mpo), i; kwargs...)
+function swap!(
+        mpo::FiniteMPO{<:MPOTensor}, i::Integer;
+        inv::Bool = false,
+        alg = SDD(), trscheme = truncbelow(eps(real(scalartype(mpo)))^(4 / 5))
+    )
+    O₁, O₂ = mpo[i], mpo[i + 1]
+
+    if inv
+        @plansor O₂₁[-1 -2 -3; -4 -5 -6] :=
+            τ'[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ[2 3; -1 -4]
+    else
+        @plansor O₂₁[-1 -2 -3; -4 -5 -6] :=
+            τ[-3 -6; 4 5] * O₁[-2 4; 2 1] * O₂[1 5; 3 -5] * τ'[2 3; -1 -4]
+    end
+
+    U, S, Vᴴ, = tsvd!(O₂₁; alg, trunc = trscheme)
+    sqrtS = sqrt(S)
+    @plansor mpo[i][-1 -2; -3 -4] := U[-3 -1 -2; 1] * sqrtS[1; -4]
+    @plansor mpo[i + 1][-1 -2; -3 -4] := sqrtS[-1; 1] * Vᴴ[1; -3 -4 -2]
+    return mpo
+end
+
+@doc """
+    multiply_neighbours(mpo::FiniteMPO, i::Integer)
+    multiply_neighbours!(mpo::FiniteMPO, i::Integer)
+
+Construct the mpo of length `length(mpo) - 1` which is formed by multiplying the operators
+on site `i` and `i + 1`.
+""" multiply_neighbours, multiply_neighbours!
+
+multiply_neighbours(mpo::FiniteMPO, i::Integer) = multiply_neighbours!(copy(mpo), i)
+function multiply_neighbours!(mpo::FiniteMPO{<:MPOTensor}, i::Integer)
+    1 <= i < length(mpo) || throw(BoundsError(mpo, i))
+    O₁ = mpo[i]
+    O₂ = popat!(parent(mpo), i + 1)
+    @plansor mpo[i][-1 -2; -3 -4] := O₁[-1 -2; 1 2] * τ[1 2; 3 4] * O₂[3 4; -3 -4]
+    return mpo
+end
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 476983658..3e04ed5fe 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -309,14 +309,21 @@ function instantiate_operator(state::AbstractMPS, O::Pair)
     return instantiate_operator(physicalspace(state), O)
 end
 function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O)::Pair)
-    inds = inds′ isa Int ? tuple(inds′) : inds′
-    mpo = O isa FiniteMPO ? O : FiniteMPO(O)
+    inds = inds′ isa Int ? [inds′] : inds′
+    mpo = O isa FiniteMPO ? copy(O) : FiniteMPO(O)
 
     # convert to linear index type
-    operators = parent(mpo)
-    indices = map(inds) do I
-        return Base._to_linear_index(lattice, Tuple(I)...) # this should mean all inds are valid...
+    indices = Vector{Int}(undef, length(inds))
+    for i in eachindex(indices)
+        indices[i] = Base._to_linear_index(lattice, Tuple(inds[i])...) # this should mean all inds are valid...
     end
+
+    # sort indices and deduplicate
+    indices, mpo = canonicalize_indices!(indices, mpo)
+    operators = parent(mpo)
+
+    @assert allunique(indices) && issorted(indices) "From here on we require unique and ascending indices\n$indices"
+
     T = eltype(mpo)
     local_mpo = Union{T, scalartype(T)}[]
     sites = Int[]
@@ -339,6 +346,25 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O
     return sites => local_mpo
 end
 
+function canonicalize_indices!(indices, mpo)
+    # swap non-sorted entries
+    for i in 2:length(indices)
+        for j in reverse(i:length(indices))
+            if indices[j] < indices[j - 1]
+                swap!(mpo, j - 1)
+                indices[j - 1], indices[j] = indices[j], indices[j - 1]
+            end
+        end
+    end
+    for i in length(indices):-1:2
+        if indices[i] == indices[i - 1]
+            multiply_neighbours!(mpo, i - 1)
+            popat!(indices, i)
+        end
+    end
+    return indices, mpo
+end
+
 # yields the promoted tensortype of all tensors
 function _find_tensortype(nonzero_operators::AbstractArray)
     return mapreduce(promote_type, nonzero_operators) do x
diff --git a/test/operators.jl b/test/operators.jl
index 61e26c67d..2b8bb7f7c 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -218,6 +218,43 @@ module TestOperators
         end
     end
 
+    @testset "Finite MPOHamiltonian repeated indices" begin
+        X = randn(ComplexF64, ℂ^2, ℂ^2)
+        X += X'
+        Y = randn(ComplexF64, ℂ^2, ℂ^2)
+        Y += Y'
+        L = 4
+        chain = fill(space(X, 1), 4)
+
+        H1 = FiniteMPOHamiltonian(chain, (1,) => (X * X * Y * Y))
+        H2 = FiniteMPOHamiltonian(chain, (1, 1, 1, 1) => (X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * Y) ⊗ (X * Y)))
+        H2 = FiniteMPOHamiltonian(chain, (1, 2, 1, 2) => (X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * X * Y) ⊗ Y))
+        H2 = FiniteMPOHamiltonian(chain, (1, 1, 1, 2) => (X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2) => ((X * Y * Y) ⊗ X))
+        H2 = FiniteMPOHamiltonian(chain, (1, 2, 1, 1) => (X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO((X * X) ⊗ Y ⊗ Y))
+        H2 = FiniteMPOHamiltonian(chain, (1, 1, 2, 3) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO((Y * Y) ⊗ X ⊗ X))
+        H2 = FiniteMPOHamiltonian(chain, (2, 3, 1, 1) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+
+        H1 = FiniteMPOHamiltonian(chain, (1, 2, 3) => FiniteMPO(X ⊗ (X * Y) ⊗ Y))
+        H2 = FiniteMPOHamiltonian(chain, (1, 2, 2, 3) => FiniteMPO(X ⊗ X ⊗ Y ⊗ Y))
+        @test convert(TensorMap, H1) ≈ convert(TensorMap, H2)
+    end
+
     @testset "InfiniteMPOHamiltonian $(sectortype(pspace))" for (pspace, Dspace) in zip(pspaces, vspaces)
         # generate a 1-2-3 body interaction
         operators = ntuple(3) do i