diff --git a/src/operators/densempo.jl b/src/operators/densempo.jl
index daf9f9f6..ca4c1584 100644
--- a/src/operators/densempo.jl
+++ b/src/operators/densempo.jl
@@ -83,12 +83,12 @@ function Base.convert(::Type{<:AbstractTensorMap}, mpo::FiniteMPO)
     U_right = Tensor(ones, scalartype(mpo), V_right')
 
     tensors = vcat(U_left, mpo.opp, U_right)
-    indices = [[i, -i, -(i + N), i + 1] for i in 1:length(mpo)]
+    indices = [[i, -i, -(2N - i + 1), i + 1] for i in 1:length(mpo)]
     pushfirst!(indices, [1])
     push!(indices, [N + 1])
     O = ncon(tensors, indices)
 
-    return transpose(O, (ntuple(identity, N), ntuple(i -> i + N, N)))
+    return transpose(O, (ntuple(identity, N), ntuple(i -> 2N - i + 1, N)))
 end
 
 # Linear Algebra
@@ -210,7 +210,8 @@ function Base.:*(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
         Fₗ = i != 1 ? Fᵣ :
              isomorphism(A, fuse(left_virtualspace(mpo2, i), left_virtualspace(mpo1, i)),
                          left_virtualspace(mpo2, i) * left_virtualspace(mpo1, i))
-        Fᵣ = isomorphism(A, fuse(right_virtualspace(mpo2, i), right_virtualspace(mpo1, i)),
+        Fᵣ = isomorphism(A,
+                         fuse(right_virtualspace(mpo2, i)', right_virtualspace(mpo1, i)'),
                          right_virtualspace(mpo2, i)' * right_virtualspace(mpo1, i)')
 
         @plansor O[i][-1 -2; -3 -4] := Fₗ[-1; 1 4] * mpo2[i][1 2; -3 3] *
@@ -232,7 +233,7 @@ function Base.:*(mpo::FiniteMPO, mps::FiniteMPS)
         Fₗ = i != 1 ? Fᵣ :
              isomorphism(TT, fuse(left_virtualspace(A[i]), left_virtualspace(mpo, i)),
                          left_virtualspace(A[i]) * left_virtualspace(mpo, i))
-        Fᵣ = isomorphism(TT, fuse(right_virtualspace(A[i]), right_virtualspace(mpo, i)),
+        Fᵣ = isomorphism(TT, fuse(right_virtualspace(A[i])', right_virtualspace(mpo, i)'),
                          right_virtualspace(A[i])' * right_virtualspace(mpo, i)')
         @plansor A[i][-1 -2; -3] := Fₗ[-1; 1 3] * A[i][1 2; 4] *
                                     mpo[i][3 -2; 2 5] *
diff --git a/test/operators.jl b/test/operators.jl
index 7c105fce..a0dac1d4 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -19,39 +19,43 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 // 2 => 10, 3 // 2 => 5, 5
 @testset "FiniteMPO" begin
     # start from random operators
     L = 4
-    O₁ = TensorMap(rand, ComplexF64, (ℂ^2)^L, (ℂ^2)^L)
-    O₂ = TensorMap(rand, ComplexF64, space(O₁))
+    T = ComplexF64
 
-    # create MPO and convert it back to see if it is the same
-    mpo₁ = FiniteMPO(O₁) # type-unstable for now!
-    mpo₂ = FiniteMPO(O₂)
-    @test convert(TensorMap, mpo₁) ≈ O₁
-    @test convert(TensorMap, -mpo₂) ≈ -O₂
+    for V in (ℂ^2, U1Space(0 => 1, 1 => 1))
+        O₁ = TensorMap(rand, T, V^L, V^L)
+        O₂ = TensorMap(rand, T, space(O₁))
 
-    # test scalar multiplication
-    α = rand(ComplexF64)
-    @test convert(TensorMap, α * mpo₁) ≈ α * O₁
-    @test convert(TensorMap, mpo₁ * α) ≈ O₁ * α
+        # create MPO and convert it back to see if it is the same
+        mpo₁ = FiniteMPO(O₁) # type-unstable for now!
+        mpo₂ = FiniteMPO(O₂)
+        @test convert(TensorMap, mpo₁) ≈ O₁
+        @test convert(TensorMap, -mpo₂) ≈ -O₂
 
-    # test addition and multiplication
-    @test convert(TensorMap, mpo₁ + mpo₂) ≈ O₁ + O₂
-    @test convert(TensorMap, mpo₁ * mpo₂) ≈ O₁ * O₂
+        # test scalar multiplication
+        α = rand(T)
+        @test convert(TensorMap, α * mpo₁) ≈ α * O₁
+        @test convert(TensorMap, mpo₁ * α) ≈ O₁ * α
 
-    # test application to a state
-    ψ₁ = Tensor(rand, ComplexF64, domain(O₁))
-    mps₁ = FiniteMPS(ψ₁)
+        # test addition and multiplication
+        @test convert(TensorMap, mpo₁ + mpo₂) ≈ O₁ + O₂
+        @test convert(TensorMap, mpo₁ * mpo₂) ≈ O₁ * O₂
 
-    @test convert(TensorMap, mpo₁ * mps₁) ≈ O₁ * ψ₁
+        # test application to a state
+        ψ₁ = Tensor(rand, T, domain(O₁))
+        mps₁ = FiniteMPS(ψ₁)
 
-    @test dot(mps₁, mpo₁, mps₁) ≈ dot(ψ₁, O₁, ψ₁)
-    @test dot(mps₁, mpo₁, mps₁) ≈ dot(mps₁, mpo₁ * mps₁)
-    # test conversion to and from mps
-    mpomps₁ = convert(FiniteMPS, mpo₁)
-    mpompsmpo₁ = convert(FiniteMPO, mpomps₁)
+        @test convert(TensorMap, mpo₁ * mps₁) ≈ O₁ * ψ₁
 
-    @test convert(FiniteMPO, mpomps₁) ≈ mpo₁ rtol = 1e-6
+        @test dot(mps₁, mpo₁, mps₁) ≈ dot(ψ₁, O₁, ψ₁)
+        @test dot(mps₁, mpo₁, mps₁) ≈ dot(mps₁, mpo₁ * mps₁)
+        # test conversion to and from mps
+        mpomps₁ = convert(FiniteMPS, mpo₁)
+        mpompsmpo₁ = convert(FiniteMPO, mpomps₁)
 
-    @test dot(mpomps₁, mpomps₁) ≈ dot(mpo₁, mpo₁)
+        @test convert(FiniteMPO, mpomps₁) ≈ mpo₁ rtol = 1e-6
+
+        @test dot(mpomps₁, mpomps₁) ≈ dot(mpo₁, mpo₁)
+    end
 end
 
 @testset "Finite MPOHamiltonian" begin