diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl
index 76d32a37e..0e6a07aae 100644
--- a/src/algorithms/correlators.jl
+++ b/src/algorithms/correlators.jl
@@ -38,6 +38,7 @@ function correlator_horizontal(
             )
             T = TransferMatrix(Atop, sandwich, _dag(Abot))
             Vo = Vo * T
+            twistdual!(T.below, 2:numout(T.below))
             Vn = Vn * T
             i += CartesianIndex(0, 1)
         end
@@ -53,6 +54,7 @@ function correlator_horizontal(
         )
         T = TransferMatrix(Atop, sandwich, _dag(Abot))
         Vo = Vo * T
+        twistdual!(T.below, 2:numout(T.below))
         Vn = Vn * T
         i += CartesianIndex(0, 1)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 787964a9b..4f5f9ea56 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -84,6 +84,9 @@ end
         @time @safetestset "Rotation of InfiniteWeightPEPS" begin
             include("utility/iwpeps_rotation.jl")
         end
+        @time @safetestset "Correlators" begin
+            include("utility/correlator.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "EXAMPLES"
         @time @safetestset "Transverse Field Ising model" begin
diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl
new file mode 100644
index 000000000..18943dfee
--- /dev/null
+++ b/test/utility/correlator.jl
@@ -0,0 +1,26 @@
+using Test
+using Random
+using TensorKit
+using PEPSKit
+import MPSKitModels: TJOperators as tJ
+
+Pspace = tJ.tj_space(Trivial, Trivial)
+Vspace = Vect[FermionParity](0 => 2, 1 => 2)
+Espace = Vect[FermionParity](0 => 3, 1 => 3)
+Random.seed!(100)
+peps = InfinitePEPS(rand, ComplexF64, Pspace, Vspace; unitcell=(2, 2));
+env = CTMRGEnv(rand, ComplexF64, peps, Espace);
+lattice = collect(space(t, 1) for t in peps.A)
+
+site0 = CartesianIndex(1, 1)
+maxsep = 8
+site1s = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep)
+
+op = tJ.S_exchange(ComplexF64, Trivial, Trivial);
+
+vals1 = correlator(peps, op, site0, site1s, env)
+vals2 = collect(begin
+    O = LocalOperator(lattice, (site0, site1) => op)
+    val = expectation_value(peps, O, env)
+end for site1 in site1s)
+@test vals1 ≈ vals2