diff --git a/src/treetensornetworks/solvers/contract_mpo_mps.jl b/src/treetensornetworks/solvers/contract_mpo_mps.jl
index ccabbc3c..e8e102aa 100644
--- a/src/treetensornetworks/solvers/contract_mpo_mps.jl
+++ b/src/treetensornetworks/solvers/contract_mpo_mps.jl
@@ -30,16 +30,20 @@ function ITensors.contract(
   init_mps = deepcopy(init_mps)
   init_mps = sim(linkinds, init_mps)
   Ai = siteinds(A)
-  ti = Vector{Index}(undef, n)
+  init_mpsi = siteinds(init_mps)
   for j in 1:n
+    ti = nothing
     for i in Ai[j]
       if !hasind(psi0[j], i)
-        ti[j] = i
+        ti = i
         break
       end
     end
+    if ti !== nothing
+      ci = commoninds(init_mpsi[j], A[j])[1]
+      replaceind!(init_mps[j], ci=>ti)
+    end
   end
-  replace_siteinds!(init_mps, ti)
 
   t = Inf
   reverse_step = false
diff --git a/test/test_treetensornetworks/test_solvers/test_contract_mpo.jl b/test/test_treetensornetworks/test_solvers/test_contract_mpo.jl
index dc3564d7..937fe7e9 100644
--- a/test/test_treetensornetworks/test_solvers/test_contract_mpo.jl
+++ b/test/test_treetensornetworks/test_solvers/test_contract_mpo.jl
@@ -51,4 +51,29 @@ using Test
   @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4
 end
 
+function asMPS(M::MPO, sites)
+  M_ = MPS(length(sites))
+  for n in eachindex(sites)
+    M_[n] = M[n]
+  end
+  return M_
+end
+
+@testset "Contract MPO-MPO" begin
+  nbit = 5
+  sites = siteinds("Qubit", nbit)
+  M1 = randomMPO(sites) + randomMPO(sites)
+  M2 = randomMPO(sites) + randomMPO(sites)
+
+  # The function `apply` does not work correctly with the mapping-MPO-to-MPS trick.
+  M1 = replaceprime(M1, 1=>2, 0=>1)
+
+  M2_ = asMPS(M2, sites)
+
+  M12_ref = asMPS(contract(M1, M2; alg="naive"), sites)
+  M12 = contract(M1, M2_; alg="fit")
+
+  @test M12_ref ≈ M12
+end
+
 nothing