diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index adc18395a..406d85613 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -24,6 +24,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Changed
 
+- In `correlator`, sites `js` are allowed on either side of the reference site `i` by adding right-to-left contractions (#399)
+
 ### Deprecated
 
 ### Removed
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index cd0974636..14f59873c 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -96,6 +96,7 @@ include("algorithms/contractions/bondenv/gaugefix.jl")
 include("algorithms/contractions/bondenv/als_solve.jl")
 include("algorithms/contractions/bondenv/benv_ctm.jl")
 include("algorithms/contractions/correlator/peps.jl")
+include("algorithms/contractions/correlator/pepo_purified.jl")
 include("algorithms/contractions/correlator/pepo_1layer.jl")
 
 include("algorithms/ctmrg/sparse_environments.jl")
@@ -128,6 +129,7 @@ include("algorithms/bp/gaugefix.jl")
 
 include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
+include("algorithms/correlator_adapters.jl")
 include("algorithms/correlators.jl")
 
 include("algorithms/optimization/fixed_point_differentiation.jl")
diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl
index 580f271d7..83924ff35 100644
--- a/src/algorithms/contractions/correlator/pepo_1layer.jl
+++ b/src/algorithms/contractions/correlator/pepo_1layer.jl
@@ -1,4 +1,6 @@
-function start_correlator(
+# -------- For left-to-right correlator contraction --------
+
+function start_correlator_left(
         i::CartesianIndex{2}, ρ::InfinitePEPO,
         O::PFTensor, env::CTMRGEnv
     )
@@ -25,8 +27,8 @@ function start_correlator(
     return Vn, Vo
 end
 
-function end_correlator_numerator(
-        j::CartesianIndex{2}, V::CTMRGEdgeTensor{T, S, 3},
+function end_correlator_right_numerator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1},
         ρ::InfinitePEPO, O::PFTensor, env::CTMRGEnv
     ) where {T, S}
     (size(ρ, 3) == 1) ||
@@ -44,8 +46,8 @@ function end_correlator_numerator(
         t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1]
 end
 
-function end_correlator_denominator(
-        j::CartesianIndex{2}, V::CTMRGEdgeTensor{T, S, 2}, env::CTMRGEnv
+function end_correlator_right_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 2, 1}, env::CTMRGEnv
     ) where {T, S}
     r, c = Tuple(j)
     C_northeast = corner(env, NORTHEAST, r - 1, c + 1)
@@ -54,3 +56,62 @@ function end_correlator_denominator(
     return @autoopt @tensor V[χS DE; χN] * C_northeast[χN; χNE] *
         E_east[χNE DE; χSE] * C_southeast[χSE; χS]
 end
+
+# -------- For right-to-left correlator contraction --------
+
+function start_correlator_right(
+        i::CartesianIndex{2}, ρ::InfinitePEPO,
+        O::PFTensor, env::CTMRGEnv
+    )
+    (size(ρ, 3) == 1) ||
+        throw(ArgumentError("The input PEPO ρ must have only one layer."))
+    r, c = Tuple(i)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_east = edge(env, EAST, r, c + 1)
+    E_south = edge(env, SOUTH, r + 1, c)
+    C_northeast = corner(env, NORTHEAST, r - 1, c + 1)
+    C_southeast = corner(env, SOUTHEAST, r + 1, c + 1)
+    t = twistdual(ρ[r, c], 1:2)
+    @autoopt @tensor Vn[χNW DW; χSW] :=
+        E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] *
+        E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] * t[d d; DN DE DS DW]
+    @autoopt @tensor Vo[χNW DW dstring; χSW] :=
+        E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] *
+        E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        removeunit(O, 1)[d2; d1 dstring] * t[d1 d2; DN DE DS DW]
+    return Vn, Vo
+end
+
+
+function end_correlator_left_numerator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1},
+        ρ::InfinitePEPO, O::PFTensor, env::CTMRGEnv
+    ) where {T, S}
+    (size(ρ, 3) == 1) ||
+        throw(ArgumentError("The input PEPO ρ must have only one layer."))
+    r, c = Tuple(j)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_south = edge(env, SOUTH, r + 1, c)
+    E_west = edge(env, WEST, r, c - 1)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+    t = twistdual(ρ[r, c], 1:2)
+    return @autoopt @tensor V[χNE DE dstring; χSE] *
+        E_south[χSE DS; χSW2] * C_southwest[χSW2; χSW] *
+        E_west[χSW DW; χNW] * C_northwest[χNW; χN] *
+        E_north[χN DN; χNE] *
+        t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1]
+end
+
+function end_correlator_left_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 2, 1}, env::CTMRGEnv
+    ) where {T, S}
+    r, c = Tuple(j)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    E_west = edge(env, WEST, r, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+    return @autoopt @tensor V[χN DE; χS] * C_southwest[χS; χSW] *
+        E_west[χSW DE; χNW] * C_northwest[χNW; χN]
+end
diff --git a/src/algorithms/contractions/correlator/pepo_purified.jl b/src/algorithms/contractions/correlator/pepo_purified.jl
new file mode 100644
index 000000000..362342b5f
--- /dev/null
+++ b/src/algorithms/contractions/correlator/pepo_purified.jl
@@ -0,0 +1,135 @@
+# -------- For left-to-right correlator contraction --------
+
+function start_correlator_left(
+        i::CartesianIndex{2},
+        below::InfinitePEPO,
+        O::MPOTensor,
+        above::InfinitePEPO,
+        env::CTMRGEnv,
+    )
+    r, c = Tuple(i)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_south = edge(env, SOUTH, r + 1, c)
+    E_west = edge(env, WEST, r, c - 1)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+    ket = twistdual(below[r, c], (1, 2))
+    bra = above[r, c]
+
+    @autoopt @tensor Vn[χSE Detop Debot; χNE] :=
+        E_south[χSE Dstop Dsbot; χSW2] *
+        C_southwest[χSW2; χSW] *
+        E_west[χSW Dwtop Dwbot; χNW] *
+        C_northwest[χNW; χN] *
+        conj(bra[d a; Dnbot Debot Dsbot Dwbot]) *
+        ket[d a; Dntop Detop Dstop Dwtop] *
+        E_north[χN Dntop Dnbot; χNE]
+
+    @autoopt @tensor Vo[χSE Detop dstring Debot; χNE] :=
+        E_south[χSE Dstop Dsbot; χSW2] *
+        C_southwest[χSW2; χSW] *
+        E_west[χSW Dwtop Dwbot; χNW] *
+        C_northwest[χNW; χN] *
+        conj(bra[d1 a; Dnbot Debot Dsbot Dwbot]) *
+        removeunit(O, 1)[d1; d2 dstring] *
+        ket[d2 a; Dntop Detop Dstop Dwtop] *
+        E_north[χN Dntop Dnbot; χNE]
+
+    return Vn, Vo
+end
+
+function end_correlator_right_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 4, 1},
+        above::InfinitePEPO,
+        O::MPOTensor,
+        below::InfinitePEPO,
+        env::CTMRGEnv,
+    ) where {T, S}
+    r, c = Tuple(j)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_east = edge(env, EAST, r, c + 1)
+    E_south = edge(env, SOUTH, r + 1, c)
+    C_northeast = corner(env, NORTHEAST, r - 1, c + 1)
+    C_southeast = corner(env, SOUTHEAST, r + 1, c + 1)
+    ket = twistdual(above[r, c], (1, 2))
+    bra = below[r, c]
+
+    return @autoopt @tensor V[χSW DWt dstring DWb; χNW] *
+        E_south[χSSE DSt DSb; χSW] *
+        E_east[χNEE DEt DEb; χSEE] *
+        E_north[χNW DNt DNb; χNNE] *
+        C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        conj(bra[db a; DNb DEb DSb DWb]) *
+        ket[dt a; DNt DEt DSt DWt] *
+        removeunit(O, 4)[dstring db; dt]
+end
+
+# -------- For right-to-left correlator contraction --------
+
+function start_correlator_right(
+        i::CartesianIndex{2},
+        below::InfinitePEPO,
+        O::MPOTensor,
+        above::InfinitePEPO,
+        env::CTMRGEnv,
+    )
+    r, c = Tuple(i)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_east = edge(env, EAST, r, c + 1)
+    E_south = edge(env, SOUTH, r + 1, c)
+    C_northeast = corner(env, NORTHEAST, r - 1, c + 1)
+    C_southeast = corner(env, SOUTHEAST, r + 1, c + 1)
+    ket = twistdual(below[r, c], (1, 2))
+    bra = above[r, c]
+
+    @autoopt @tensor Vn[χNW DWtop DWbot; χSW] :=
+        E_south[χSSE Dstop Dsbot; χSW] *
+        E_east[χNEE Detop Debot; χSEE] *
+        E_north[χNW Dntop Dnbot; χNNE] *
+        C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        conj(bra[d a; Dnbot Debot Dsbot DWbot]) *
+        ket[d a; Dntop Detop Dstop DWtop]
+
+    @autoopt @tensor Vo[χNW DWtop dstring DWbot; χSW] :=
+        E_south[χSSE Dstop Dsbot; χSW] *
+        E_east[χNEE Detop Debot; χSEE] *
+        E_north[χNW Dntop Dnbot; χNNE] *
+        C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        conj(bra[d1 a; Dnbot Debot Dsbot DWbot]) *
+        removeunit(O, 1)[d1; d2 dstring] *
+        ket[d2 a; Dntop Detop Dstop DWtop]
+
+    return Vn, Vo
+end
+
+function end_correlator_left_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 4, 1},
+        above::InfinitePEPO,
+        O::MPOTensor,
+        below::InfinitePEPO,
+        env::CTMRGEnv,
+    ) where {T, S}
+    r, c = Tuple(j)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_south = edge(env, SOUTH, r + 1, c)
+    E_west = edge(env, WEST, r, c - 1)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+    ket = twistdual(above[r, c], (1, 2))
+    bra = below[r, c]
+
+    return @autoopt @tensor V[χNE DEt dstring DEb; χSE] *
+        E_south[χSE DSt DSb; χSW2] *
+        C_southwest[χSW2; χSW] *
+        E_west[χSW DWt DWb; χNW] *
+        C_northwest[χNW; χN] *
+        E_north[χN DNt DNb; χNE] *
+        conj(bra[db a; DNb DEb DSb DWb]) *
+        ket[dt a; DNt DEt DSt DWt] *
+        removeunit(O, 4)[dstring db; dt]
+end
diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl
index fafdd938a..39392d709 100644
--- a/src/algorithms/contractions/correlator/peps.jl
+++ b/src/algorithms/contractions/correlator/peps.jl
@@ -1,4 +1,6 @@
-function start_correlator(
+# -------- For left-to-right correlator contraction --------
+
+function start_correlator_left(
         i::CartesianIndex{2},
         below::InfinitePEPS,
         O::MPOTensor,
@@ -37,7 +39,7 @@ function start_correlator(
     return Vn, Vo
 end
 
-function end_correlator_numerator(
+function end_correlator_right_numerator(
         j::CartesianIndex{2},
         V::AbstractTensorMap{T, S, 4, 1},
         above::InfinitePEPS,
@@ -64,7 +66,7 @@ function end_correlator_numerator(
         removeunit(O, 4)[dstring db; dt]
 end
 
-function end_correlator_denominator(
+function end_correlator_right_denominator(
         j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1},
         env::CTMRGEnv
     ) where {T, S}
@@ -78,3 +80,84 @@ function end_correlator_denominator(
         E_east[χNE DEt DEb; χSE] *
         C_southeast[χSE; χS]
 end
+
+# -------- For right-to-left correlator contraction --------
+
+function start_correlator_right(
+        i::CartesianIndex{2},
+        below::InfinitePEPS,
+        O::MPOTensor,
+        above::InfinitePEPS,
+        env::CTMRGEnv,
+    )
+    r, c = Tuple(i)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_east = edge(env, EAST, r, c + 1)
+    E_south = edge(env, SOUTH, r + 1, c)
+    C_northeast = corner(env, NORTHEAST, r - 1, c + 1)
+    C_southeast = corner(env, SOUTHEAST, r + 1, c + 1)
+    sandwich = (below[r, c], above[r, c])
+
+    @autoopt @tensor Vn[χNW DWtop DWbot; χSW] :=
+        E_south[χSSE Dstop Dsbot; χSW] *
+        E_east[χNEE Detop Debot; χSEE] *
+        E_north[χNW Dntop Dnbot; χNNE] *
+        C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        conj(bra(sandwich)[d; Dnbot Debot Dsbot DWbot]) *
+        ket(sandwich)[d; Dntop Detop Dstop DWtop]
+
+    @autoopt @tensor Vo[χNW DWtop dstring DWbot; χSW] :=
+        E_south[χSSE Dstop Dsbot; χSW] *
+        E_east[χNEE Detop Debot; χSEE] *
+        E_north[χNW Dntop Dnbot; χNNE] *
+        C_northeast[χNNE; χNEE] *
+        C_southeast[χSEE; χSSE] *
+        conj(bra(sandwich)[d1; Dnbot Debot Dsbot DWbot]) *
+        removeunit(O, 1)[d1; d2 dstring] *
+        ket(sandwich)[d2; Dntop Detop Dstop DWtop]
+
+    return Vn, Vo
+end
+
+function end_correlator_left_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 4, 1},
+        above::InfinitePEPS,
+        O::MPOTensor,
+        below::InfinitePEPS,
+        env::CTMRGEnv,
+    ) where {T, S}
+    r, c = Tuple(j)
+    E_north = edge(env, NORTH, r - 1, c)
+    E_south = edge(env, SOUTH, r + 1, c)
+    E_west = edge(env, WEST, r, c - 1)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+    sandwich = (above[r, c], below[r, c])
+
+    return @autoopt @tensor V[χNE DEt dstring DEb; χSE] *
+        E_south[χSE DSt DSb; χSW2] *
+        C_southwest[χSW2; χSW] *
+        E_west[χSW DWt DWb; χNW] *
+        C_northwest[χNW; χN] *
+        E_north[χN DNt DNb; χNE] *
+        conj(bra(sandwich)[db; DNb DEb DSb DWb]) *
+        ket(sandwich)[dt; DNt DEt DSt DWt] *
+        removeunit(O, 4)[dstring db; dt]
+end
+
+function end_correlator_left_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1},
+        env::CTMRGEnv
+    ) where {T, S}
+    r, c = Tuple(j)
+    C_northwest = corner(env, NORTHWEST, r - 1, c - 1)
+    E_west = edge(env, WEST, r, c - 1)
+    C_southwest = corner(env, SOUTHWEST, r + 1, c - 1)
+
+    return @autoopt @tensor V[χN DEt DEb; χS] *
+        C_southwest[χS; χSW] *
+        E_west[χSW DEt DEb; χNW] *
+        C_northwest[χNW; χN]
+end
diff --git a/src/algorithms/contractions/transfer.jl b/src/algorithms/contractions/transfer.jl
index 4a88cf0d0..d21241292 100644
--- a/src/algorithms/contractions/transfer.jl
+++ b/src/algorithms/contractions/transfer.jl
@@ -21,7 +21,7 @@ Apply an edge transfer matrix to the left.
         Etop::CTMRGEdgeTensor{<:Any, S, N₂},
         Ebot::CTMRGEdgeTensor{<:Any, S, N₂}
     ) where {S, N₁, N₂}
-    t_out = tensorexpr(:v, -1, -(2:(N₁ + 1)))
+    t_out = tensorexpr(:v′, -1, -(2:(N₁ + 1)))
     t_top = tensorexpr(:Etop, 2:(N₂ + 1), -(N₁ + 1))
     t_bot = tensorexpr(:Ebot, (-1, (3:(N₂ + 1))...), 1)
     t_in = tensorexpr(:v, 1, (-(2:N₁)..., 2))
@@ -47,7 +47,7 @@ Apply an edge transfer matrix to the right.
         Etop::CTMRGEdgeTensor{<:Any, S, N₂},
         Ebot::CTMRGEdgeTensor{<:Any, S, N₂}
     ) where {S, N₁, N₂}
-    t_out = tensorexpr(:v, -1, -(2:(N₁ + 1)))
+    t_out = tensorexpr(:v′, -1, -(2:(N₁ + 1)))
     t_top = tensorexpr(:Etop, (-1, (3:(N₂ + 1))...), 1)
     t_bot = tensorexpr(:Ebot, (2, (3:(N₂ + 1))...), -(N₁ + 1))
     t_in = tensorexpr(:v, 1, (-(2:N₁)..., 2))
@@ -70,33 +70,38 @@ Apply an edge transfer matrix to the left.
 ```
 """
 function edge_transfer_left(
-        v::CTMRGEdgeTensor{<:Any, S, 3},
-        O::PEPSSandwich,
-        Etop::CTMRGEdgeTensor{<:Any, S, 3},
-        Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PEPSSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
     ) where {S}
-    @autoopt @tensor v´[χ_SE D_E_above D_E_below; χ_NE] :=
+    return @autoopt @tensor v′[χ_SE D_E_above D_E_below; χ_NE] :=
         v[χ_SW D_W_above D_W_below; χ_NW] *
         Etop[χ_NW D_N_above D_N_below; χ_NE] *
         Ebot[χ_SE D_S_above D_S_below; χ_SW] *
         ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
         conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below])
-
-    return v´
 end
 function edge_transfer_left(
-        v::CTMRGEdgeTensor{<:Any, S, 2},
-        O::PFTensor,
-        Etop::CTMRGEdgeTensor{<:Any, S, 2},
-        Ebot::CTMRGEdgeTensor{<:Any, S, 2},
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PEPOPurifiedSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+    ) where {S}
+    ket_tensor = twistdual(ket(O), (1, 2))
+    bra_tensor = bra(O)
+    return @autoopt @tensor v′[χ_SE D_E_above D_E_below; χ_NE] :=
+        v[χ_SW D_W_above D_W_below; χ_NW] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_above D_S_below; χ_SW] *
+        ket_tensor[d a; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra_tensor[d a; D_N_below D_E_below D_S_below D_W_below])
+end
+function edge_transfer_left(
+        v::AbstractTensorMap{<:Any, S, 2, 1}, O::PFTensor,
+        Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2},
     ) where {S}
-    @autoopt @tensor v´[χ_SE D_E; χ_NE] :=
+    return @autoopt @tensor v′[χ_SE D_E; χ_NE] :=
         v[χ_SW D_W; χ_NW] *
         Etop[χ_NW D_N; χ_NE] *
         Ebot[χ_SE D_S; χ_SW] *
         O[D_W D_S; D_N D_E]
-
-    return v´
 end
 
 """
@@ -113,33 +118,38 @@ Apply an edge transfer matrix to the right.
 ```
 """
 function edge_transfer_right(
-        v::CTMRGEdgeTensor{<:Any, S, 3},
-        O::PEPSSandwich,
-        Etop::CTMRGEdgeTensor{<:Any, S, 3},
-        Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PEPSSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
     ) where {S}
-    @autoopt @tensor v′[χ_NW D_W_above D_W_below; χ_SW] :=
+    return @autoopt @tensor v′[χ_NW D_W_above D_W_below; χ_SW] :=
         v[χ_NE D_E_above D_E_below; χ_SE] *
         Etop[χ_NW D_N_above D_N_below; χ_NE] *
         Ebot[χ_SE D_S_above D_S_below; χ_SW] *
         ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
         conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below])
-
-    return v′
 end
 function edge_transfer_right(
-        v::CTMRGEdgeTensor{<:Any, S, 2},
-        O::PFTensor,
-        Etop::CTMRGEdgeTensor{<:Any, S, 2},
-        Ebot::CTMRGEdgeTensor{<:Any, S, 2},
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PEPOPurifiedSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+    ) where {S}
+    ket_tensor = twistdual(ket(O), (1, 2))
+    bra_tensor = bra(O)
+    return @autoopt @tensor v′[χ_NW D_W_above D_W_below; χ_SW] :=
+        v[χ_NE D_E_above D_E_below; χ_SE] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_above D_S_below; χ_SW] *
+        ket_tensor[d a; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra_tensor[d a; D_N_below D_E_below D_S_below D_W_below])
+end
+function edge_transfer_right(
+        v::AbstractTensorMap{<:Any, S, 2, 1}, O::PFTensor,
+        Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2},
     ) where {S}
     return @autoopt @tensor v′[χ_NW D_W; χ_SW] :=
         v[χ_NE D_E; χ_SE] *
         Etop[χ_NW D_N; χ_NE] *
         Ebot[χ_SE D_S; χ_SW] *
         O[D_W D_S; D_N D_E]
-
-    return v′
 end
 
 """
@@ -156,7 +166,7 @@ Apply an edge transfer matrix to the left on an excited vector.
 ```
 """
 function edge_transfer_left(
-        v::CTMRGEdgeTensor{<:Any, S, 4}, O::PEPSSandwich,
+        v::AbstractTensorMap{<:Any, S, 4, 1}, O::PEPSSandwich,
         Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
     ) where {S}
     return @autoopt @tensor v′[χ_SE D_E_above d_string D_E_below; χ_NE] :=
@@ -167,7 +177,20 @@ function edge_transfer_left(
         conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below])
 end
 function edge_transfer_left(
-        v::CTMRGEdgeTensor{<:Any, S, 3}, O::PFTensor,
+        v::AbstractTensorMap{<:Any, S, 4, 1}, O::PEPOPurifiedSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+    ) where {S}
+    ket_tensor = twistdual(ket(O), (1, 2))
+    bra_tensor = bra(O)
+    return @autoopt @tensor v′[χ_SE D_E_above d_string D_E_below; χ_NE] :=
+        v[χ_SW D_W_above d_string D_W_below; χ_NW] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_above D_S_below; χ_SW] *
+        ket_tensor[d a; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra_tensor[d a; D_N_below D_E_below D_S_below D_W_below])
+end
+function edge_transfer_left(
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PFTensor,
         Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2},
     ) where {S}
     return @autoopt @tensor v′[χ_SE D_E d_string; χ_NE] :=
@@ -176,3 +199,51 @@ function edge_transfer_left(
         Ebot[χ_SE D_S; χ_SW] *
         O[D_W D_S; D_N D_E]
 end
+
+"""
+    transfer_right(v, O, Et, Eb)
+    
+Apply an edge transfer matrix to the right on an excited vector..
+
+```
+──Et─┐
+  │  │
+──O──v-
+  │  │
+──qƎ─┘
+```
+"""
+function edge_transfer_right(
+        v::AbstractTensorMap{<:Any, S, 4, 1}, O::PEPSSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+    ) where {S}
+    return @autoopt @tensor v′[χ_NW D_W_above d_string D_W_below; χ_SW] :=
+        v[χ_NE D_E_above d_string D_E_below; χ_SE] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_above D_S_below; χ_SW] *
+        ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below])
+end
+function edge_transfer_right(
+        v::AbstractTensorMap{<:Any, S, 4, 1}, O::PEPOPurifiedSandwich,
+        Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3},
+    ) where {S}
+    ket_tensor = twistdual(ket(O), (1, 2))
+    bra_tensor = bra(O)
+    return @autoopt @tensor v′[χ_NW D_W_above d_string D_W_below; χ_SW] :=
+        v[χ_NE D_E_above d_string D_E_below; χ_SE] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_above D_S_below; χ_SW] *
+        ket_tensor[d a; D_N_above D_E_above D_S_above D_W_above] *
+        conj(bra_tensor[d a; D_N_below D_E_below D_S_below D_W_below])
+end
+function edge_transfer_right(
+        v::AbstractTensorMap{<:Any, S, 3, 1}, O::PFTensor,
+        Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2},
+    ) where {S}
+    return @autoopt @tensor v′[χ_NW D_W d_string; χ_SW] :=
+        v[χ_NE D_E d_string; χ_SE] *
+        Etop[χ_NW D_N; χ_NE] *
+        Ebot[χ_SE D_S; χ_SW] *
+        O[D_W D_S; D_N D_E]
+end
diff --git a/src/algorithms/correlator_adapters.jl b/src/algorithms/correlator_adapters.jl
new file mode 100644
index 000000000..44e5311dd
--- /dev/null
+++ b/src/algorithms/correlator_adapters.jl
@@ -0,0 +1,262 @@
+# Correlator adapters for InfinitePEPS / InfinitePEPO
+
+struct _PEPSCorrelator{B <: InfinitePEPS, K <: InfinitePEPS, E <: CTMRGEnv}
+    bra::B
+    ket::K
+    env::E
+
+    function _PEPSCorrelator(bra::B, ket::K, env::E) where {B, K, E}
+        size(ket) == size(bra) ||
+            throw(DimensionMismatch("The ket and bra must have the same unit cell."))
+        return new{B, K, E}(bra, ket, env)
+    end
+end
+
+struct _PEPOPurifiedCorrelator{B <: InfinitePEPO, K <: InfinitePEPO, E <: CTMRGEnv}
+    bra::B
+    ket::K
+    env::E
+
+    function _PEPOPurifiedCorrelator(bra::B, ket::K, env::E) where {B, K, E}
+        size(ket) == size(bra) ||
+            throw(DimensionMismatch("The ket and bra must have the same unit cell."))
+        size(ket, 3) == size(bra, 3) == 1 ||
+            throw(ArgumentError("Purified PEPO correlators require one-layer PEPOs."))
+        return new{B, K, E}(bra, ket, env)
+    end
+end
+
+const _BraketCorrelator = Union{_PEPSCorrelator, _PEPOPurifiedCorrelator}
+
+_braket_correlator(bra::InfinitePEPS, ket::InfinitePEPS, env::CTMRGEnv) =
+    _PEPSCorrelator(bra, ket, env)
+_braket_correlator(bra::InfinitePEPO, ket::InfinitePEPO, env::CTMRGEnv) =
+    _PEPOPurifiedCorrelator(bra, ket, env)
+
+struct _PEPOTraceCorrelator{P <: InfinitePEPO, E <: CTMRGEnv}
+    ρ::P
+    env::E
+
+    function _PEPOTraceCorrelator(ρ::P, env::E) where {P, E}
+        (size(ρ, 3) == 1) ||
+            throw(ArgumentError("The input PEPO ρ must have only one layer."))
+        return new{P, E}(ρ, env)
+    end
+end
+
+function _edge_transfermatrix(row::Int, col::Int, context::_BraketCorrelator)
+    return _edge_transfermatrix(row, col, context.bra, context.ket, context.env)
+end
+
+function _edge_transfermatrix(row::Int, col::Int, context::_PEPOTraceCorrelator)
+    return _edge_transfermatrix(row, col, context.ρ, context.env)
+end
+
+function _correlator_scalartype(context::_BraketCorrelator, O::FiniteMPO)
+    return TensorOperations.promote_contract(
+        scalartype(context.bra), scalartype(context.ket),
+        scalartype(context.env), scalartype.(O)...
+    )
+end
+
+function _correlator_scalartype(context::_PEPOTraceCorrelator, O::FiniteMPO)
+    return TensorOperations.promote_contract(
+        scalartype(context.ρ), scalartype(context.env), scalartype.(O)...
+    )
+end
+
+_correlator_unitcell(context::_BraketCorrelator) = size(context.bra)[1:2]
+_correlator_unitcell(context::_PEPOTraceCorrelator) = size(context.ρ)[1:2]
+
+function Base.rotl90(context::_PEPSCorrelator)
+    rotated_bra = rotl90(context.bra)
+    rotated_ket = context.bra === context.ket ? rotated_bra : rotl90(context.ket)
+    return _PEPSCorrelator(rotated_bra, rotated_ket, rotl90(context.env))
+end
+
+function Base.rotl90(context::_PEPOPurifiedCorrelator)
+    rotated_bra = rotl90(context.bra)
+    rotated_ket = context.bra === context.ket ? rotated_bra : rotl90(context.ket)
+    return _PEPOPurifiedCorrelator(rotated_bra, rotated_ket, rotl90(context.env))
+end
+
+Base.rotl90(context::_PEPOTraceCorrelator) =
+    _PEPOTraceCorrelator(rotl90(context.ρ), rotl90(context.env))
+
+# -------- For left-to-right correlator contraction --------
+
+function _start_correlator_left(
+        i::CartesianIndex{2}, context::_BraketCorrelator, O::MPOTensor
+    )
+    return start_correlator_left(i, context.bra, O, context.ket, context.env)
+end
+
+function _start_correlator_left(
+        i::CartesianIndex{2}, context::_PEPOTraceCorrelator, O::PFTensor
+    )
+    return start_correlator_left(i, context.ρ, O, context.env)
+end
+
+function _end_correlator_right_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 4, 1},
+        context::_BraketCorrelator,
+        O::MPOTensor,
+    ) where {T, S}
+    return end_correlator_right_numerator(j, V, context.bra, O, context.ket, context.env)
+end
+
+function _end_correlator_right_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 3, 1},
+        context::_PEPOTraceCorrelator,
+        O::PFTensor,
+    ) where {T, S}
+    return end_correlator_right_numerator(j, V, context.ρ, O, context.env)
+end
+
+function _end_correlator_right_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap, context
+    )
+    return end_correlator_right_denominator(j, V, context.env)
+end
+
+function _correlator_horizontal_right!(G, targets, context, O::FiniteMPO, i::CartesianIndex{2})
+    # left start for operator and norm contractions
+    c = i # current column being handled
+    Vn, Vo = _start_correlator_left(c, context, O[1])
+    j_last = last(targets)[2]
+    for (k, j) in targets
+        local numerator
+        while j > c
+            c += CartesianIndex(0, 1)
+            if c == j
+                numerator = _end_correlator_right_numerator(j, Vo, context, O[2])
+            end
+            T = _edge_transfermatrix(c[1], c[2], context)
+            c != j_last && (Vo = Vo * T)
+            Vn = Vn * T
+        end
+        # compute overlap without operator
+        denominator = _end_correlator_right_denominator(j, Vn, context)
+        G[k] = numerator / denominator
+    end
+    return G
+end
+
+# -------- For right-to-left correlator contraction --------
+
+function _start_correlator_right(
+        i::CartesianIndex{2}, context::_BraketCorrelator, O::MPOTensor
+    )
+    return start_correlator_right(i, context.bra, O, context.ket, context.env)
+end
+
+function _start_correlator_right(
+        i::CartesianIndex{2}, context::_PEPOTraceCorrelator, O::PFTensor
+    )
+    return start_correlator_right(i, context.ρ, O, context.env)
+end
+
+function _end_correlator_left_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 4, 1},
+        context::_BraketCorrelator,
+        O::MPOTensor,
+    ) where {T, S}
+    return end_correlator_left_numerator(j, V, context.bra, O, context.ket, context.env)
+end
+
+function _end_correlator_left_numerator(
+        j::CartesianIndex{2},
+        V::AbstractTensorMap{T, S, 3, 1},
+        context::_PEPOTraceCorrelator,
+        O::PFTensor,
+    ) where {T, S}
+    return end_correlator_left_numerator(j, V, context.ρ, O, context.env)
+end
+
+function _end_correlator_left_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap, context
+    )
+    return end_correlator_left_denominator(j, V, context.env)
+end
+
+function _correlator_horizontal_left!(G, targets, context, O::FiniteMPO, i::CartesianIndex{2})
+    # right start for operator and norm contractions
+    c = i # current column being handled
+    Vn, Vo = _start_correlator_right(c, context, O[1])
+    j_last = last(targets)[2]
+    for (k, j) in targets
+        local numerator
+        while j < c
+            c -= CartesianIndex(0, 1)
+            if c == j
+                numerator = _end_correlator_left_numerator(j, Vo, context, O[2])
+            end
+            T = _edge_transfermatrix(c[1], c[2], context)
+            c != j_last && (Vo = T * Vo)
+            Vn = T * Vn
+        end
+        # compute overlap without operator
+        denominator = _end_correlator_left_denominator(j, Vn, context)
+        G[k] = numerator / denominator
+    end
+    return G
+end
+
+# -------- main correlator implementation --------
+
+function _check_horizontal_correlator_sites(
+        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}
+    )
+    all(==(i[1]) ∘ first ∘ Tuple, js) ||
+        throw(ArgumentError("Not a horizontal correlation function"))
+    all(!=(i), js) ||
+        throw(ArgumentError("Correlator target sites must differ from the reference site"))
+    allunique(js) ||
+        throw(ArgumentError("Correlator target sites must be unique"))
+    return true
+end
+
+function _correlator_horizontal(
+        context, operator,
+        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}},
+    )
+    _check_horizontal_correlator_sites(i, js)
+    O = FiniteMPO(operator)
+    length(O) == 2 || throw(ArgumentError("Operator must act on two sites"))
+    # preallocate with correct scalartype
+    G = similar(js, _correlator_scalartype(context, O))
+    # calculate correlator on the right and left separately
+    right_targets = sort(
+        filter(((k, j),) -> j[2] > i[2], collect(enumerate(js)));
+        by = x -> x[2][2]
+    )
+    left_targets = sort(
+        filter(((k, j),) -> j[2] < i[2], collect(enumerate(js)));
+        by = x -> x[2][2], rev = true,
+    )
+    isempty(right_targets) || _correlator_horizontal_right!(G, right_targets, context, O, i)
+    isempty(left_targets) || _correlator_horizontal_left!(G, left_targets, context, O, i)
+    return G
+end
+
+function _correlator_vertical(context, operator, i::CartesianIndex{2}, js)
+    unitcell = _correlator_unitcell(context)
+    rotated_i = siterotl90(i, unitcell)
+    rotated_js = map(j -> siterotl90(j, unitcell), js)
+    return _correlator_horizontal(rotl90(context), operator, rotated_i, rotated_js)
+end
+
+function _correlator(
+        context, O, i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}
+    )
+    if all(==(i[1]) ∘ first ∘ Tuple, js)
+        return _correlator_horizontal(context, O, i, js)
+    elseif all(==(i[2]) ∘ last ∘ Tuple, js)
+        return _correlator_vertical(context, O, i, js)
+    else
+        error("Only horizontal or vertical correlators are implemented")
+    end
+end
diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl
index 0c95723b6..ab7c5641f 100644
--- a/src/algorithms/correlators.jl
+++ b/src/algorithms/correlators.jl
@@ -1,194 +1,40 @@
 const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianIndices{N}}
 
-# Correlators in InfinitePEPS
-
-function correlator_horizontal(
-        bra::InfinitePEPS,
-        operator,
-        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}},
-        ket::InfinitePEPS,
-        env::CTMRGEnv,
-    )
-    size(ket) == size(bra) ||
-        throw(DimensionMismatch("The ket and bra must have the same unit cell."))
-    all(==(i[1]) ∘ first ∘ Tuple, js) ||
-        throw(ArgumentError("Not a horizontal correlation function"))
-    issorted(vcat(i, js); by = last ∘ Tuple) ||
-        throw(ArgumentError("Not an increasing sequence of coordinates"))
-    O = FiniteMPO(operator)
-    length(O) == 2 || throw(ArgumentError("Operator must act on two sites"))
-    # preallocate with correct scalartype
-    G = similar(
-        js,
-        TensorOperations.promote_contract(
-            scalartype(bra), scalartype(ket), scalartype(env), scalartype.(O)...
-        ),
-    )
-    # left start for operator and norm contractions
-    Vn, Vo = start_correlator(i, bra, O[1], ket, env)
-    i += CartesianIndex(0, 1)
-    for (k, j) in enumerate(js)
-        # transfer until left of site j
-        while j > i
-            Etop = edge(env, NORTH, i[1] - 1, i[2])
-            Ebot = edge(env, SOUTH, i[1] + 1, i[2])
-            sandwich = (
-                ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)],
-            )
-            T = edge_transfermatrix(Etop, sandwich, Ebot)
-            Vo = Vo * T
-            Vn = Vn * T
-            i += CartesianIndex(0, 1)
-        end
-        # compute overlap with operator
-        numerator = end_correlator_numerator(j, Vo, bra, O[2], ket, env)
-        # transfer right of site j
-        Etop = edge(env, NORTH, i[1] - 1, i[2])
-        Ebot = edge(env, SOUTH, i[1] + 1, i[2])
-        sandwich = (ket[i[1], i[2]], bra[i[1], i[2]])
-        T = edge_transfermatrix(Etop, sandwich, Ebot)
-        if k < length(js)
-            Vo = Vo * T
-        end
-        Vn = Vn * T
-        i += CartesianIndex(0, 1)
-        # compute overlap without operator
-        denominator = end_correlator_denominator(j, Vn, env)
-        G[k] = numerator / denominator
-    end
-    return G
-end
-
-function correlator_vertical(
-        bra::InfinitePEPS,
-        operator,
-        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}},
-        ket::InfinitePEPS,
-        env::CTMRGEnv,
-    )
-    rotated_bra = rotl90(bra)
-    rotated_ket = bra === ket ? rotated_bra : rotl90(ket)
-    rotated_i = siterotl90(i, size(bra))
-    rotated_js = map(j -> siterotl90(j, size(bra)), js)
-    return correlator_horizontal(
-        rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env)
-    )
-end
+# Correlators in InfinitePEPS or purified InfinitePEPO
 
 function MPSKit.correlator(
-        bra::InfinitePEPS,
+        bra::S,
         O,
         i::CartesianIndex{2}, js::CoordCollection{2},
-        ket::InfinitePEPS,
+        ket::S,
         env::CTMRGEnv,
-    )
-    js = vec(js) # map CartesianIndices to actual Vector instead of Matrix
-
-    if all(==(i[1]) ∘ first ∘ Tuple, js)
-        return correlator_horizontal(bra, O, i, js, ket, env)
-    elseif all(==(i[2]) ∘ last ∘ Tuple, js)
-        return correlator_vertical(bra, O, i, js, ket, env)
-    else
-        error("Only horizontal or vertical correlators are implemented")
-    end
+    ) where {S <: InfiniteState}
+    return _correlator(_braket_correlator(bra, ket, env), O, i, vec(js))
 end
 
 function MPSKit.correlator(
-        bra::InfinitePEPS,
+        bra::S,
         O,
         i::CartesianIndex{2}, j::CartesianIndex{2},
-        ket::InfinitePEPS,
+        ket::S,
         env::CTMRGEnv,
-    )
+    ) where {S <: InfiniteState}
     return only(correlator(bra, O, i, j:j, ket, env))
 end
 
+## reserved for InfinitePEPS
 function MPSKit.correlator(state::InfinitePEPS, O, i::CartesianIndex{2}, j, env::CTMRGEnv)
     return MPSKit.correlator(state, O, i, j, state, env)
 end
 
 # Correlators in InfinitePEPO (tr(ρO))
 
-function correlator_horizontal(
-        ρ::InfinitePEPO, operator,
-        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}},
-        env::CTMRGEnv
-    )
-    (size(ρ, 3) == 1) ||
-        throw(ArgumentError("The input PEPO ρ must have only one layer."))
-    all(==(i[1]) ∘ first ∘ Tuple, js) ||
-        throw(ArgumentError("Not a horizontal correlation function"))
-    issorted(vcat(i, js); by = last ∘ Tuple) ||
-        throw(ArgumentError("Not an increasing sequence of coordinates"))
-    O = FiniteMPO(operator)
-    length(O) == 2 || throw(ArgumentError("Operator must act on two sites"))
-    # preallocate with correct scalartype
-    G = similar(
-        js,
-        TensorOperations.promote_contract(
-            scalartype(ρ), scalartype(env), scalartype.(O)...
-        ),
-    )
-    # left start for operator and norm contractions
-    Vn, Vo = start_correlator(i, ρ, O[1], env)
-    i += CartesianIndex(0, 1)
-    for (k, j) in enumerate(js)
-        # transfer until left of site j
-        while j > i
-            Etop = edge(env, NORTH, i[1] - 1, i[2])
-            Omid = trace_physicalspaces(ρ[i[1], i[2]])
-            Ebot = edge(env, SOUTH, i[1] + 1, i[2])
-            T = edge_transfermatrix(Etop, Omid, Ebot)
-            Vo = Vo * T
-            Vn = Vn * T
-            i += CartesianIndex(0, 1)
-        end
-        # compute overlap with operator
-        numerator = end_correlator_numerator(j, Vo, ρ, O[2], env)
-        # transfer right of site j
-        Etop = edge(env, NORTH, i[1] - 1, i[2])
-        Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)])
-        Ebot = edge(env, SOUTH, i[1] + 1, i[2])
-        T = edge_transfermatrix(Etop, Omid, Ebot)
-        if k < length(js)
-            Vo = Vo * T
-        end
-        Vn = Vn * T
-        i += CartesianIndex(0, 1)
-        # compute overlap without operator
-        denominator = end_correlator_denominator(j, Vn, env)
-        G[k] = numerator / denominator
-    end
-    return G
-end
-
-function correlator_vertical(
-        ρ::InfinitePEPO, operator,
-        i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}},
-        env::CTMRGEnv,
-    )
-    rotated_ρ = rotl90(ρ)
-    unitcell = size(ρ)[1:2]
-    rotated_i = siterotl90(i, unitcell)
-    rotated_js = map(j -> siterotl90(j, unitcell), js)
-    return correlator_horizontal(
-        rotated_ρ, operator, rotated_i, rotated_js, rotl90(env)
-    )
-end
-
 function MPSKit.correlator(
         ρ::InfinitePEPO, O,
         i::CartesianIndex{2}, js::CoordCollection{2},
         env::CTMRGEnv,
     )
-    js = vec(js) # map CartesianIndices to Vector instead of Matrix
-    if all(==(i[1]) ∘ first ∘ Tuple, js)
-        return correlator_horizontal(ρ, O, i, js, env)
-    elseif all(==(i[2]) ∘ last ∘ Tuple, js)
-        return correlator_vertical(ρ, O, i, js, env)
-    else
-        error("Only horizontal or vertical correlators are implemented")
-    end
+    return _correlator(_PEPOTraceCorrelator(ρ, env), O, i, vec(js))
 end
 
 function MPSKit.correlator(
@@ -198,5 +44,3 @@ function MPSKit.correlator(
     )
     return only(correlator(ρ, O, i, j:j, env))
 end
-
-# TODO: Correlators in InfinitePEPO (⟨ρ|O|ρ⟩)
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index c532d3e1d..b53e0f443 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -8,9 +8,8 @@ In the latter case the first signature corresponds to a single layer PEPO contra
 the second signature yields a bilayer contraction instead.
 """
 function MPSKit.expectation_value(
-        bra::Union{InfinitePEPS, InfinitePEPO}, O::LocalOperator,
-        ket::Union{InfinitePEPS, InfinitePEPO}, env::CTMRGEnv
-    )
+        bra::S, O::LocalOperator, ket::S, env::CTMRGEnv
+    ) where {S <: InfiniteState}
     checklattice(bra, O, ket)
     term_vals = dtmap(collect(O.terms)) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
         ρ = reduced_densitymatrix(inds, ket, bra, env)
diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl
index 354d3fe73..19e045e07 100644
--- a/src/algorithms/transfermatrix.jl
+++ b/src/algorithms/transfermatrix.jl
@@ -12,7 +12,7 @@ end
 
 Base.:*(tm1::T, tm2::T) where {T <: EdgeTransferMatrix} = ProductTransferMatrix([tm1, tm2])
 
-# TODO: really not sure it TensorKit.flip is the suitable method for this...
+# TODO: really not sure if TensorKit.flip is the suitable method for this...
 function TensorKit.flip(tm::EdgeTransferMatrix)
     return EdgeTransferMatrix(tm.top, tm.mid, tm.bot, !tm.isflipped)
 end
@@ -36,3 +36,21 @@ function edge_transfermatrix(a::AbstractVector, b, c::AbstractVector, isflipped
     tot = ProductTransferMatrix(convert(Vector, edge_transfermatrix.(a, b, c)))
     return isflipped ? flip(tot) : tot
 end
+
+function _edge_transfermatrix(
+        row::Int, col::Int, bra::S, ket::S, env::CTMRGEnv
+    ) where {S <: InfiniteState}
+    Etop = edge(env, NORTH, row - 1, col)
+    Ebot = edge(env, SOUTH, row + 1, col)
+    sandwich = (ket[row, col], bra[row, col])
+    return edge_transfermatrix(Etop, sandwich, Ebot)
+end
+
+function _edge_transfermatrix(
+        row::Int, col::Int, ρ::InfinitePEPO, env::CTMRGEnv
+    )
+    Etop = edge(env, NORTH, row - 1, col)
+    Ebot = edge(env, SOUTH, row + 1, col)
+    Omid = trace_physicalspaces(ρ[row, col])
+    return edge_transfermatrix(Etop, Omid, Ebot)
+end
diff --git a/src/networks/local_sandwich.jl b/src/networks/local_sandwich.jl
index 4c46cd9fe..90386e9fe 100644
--- a/src/networks/local_sandwich.jl
+++ b/src/networks/local_sandwich.jl
@@ -108,8 +108,11 @@ end
 
 ## PEPO
 
-const PEPOSandwich{N, T <: PEPSTensor, P <: PEPOTensor} = Tuple{T, T, Vararg{P, N}}
+const PEPOPurifiedSandwich{K <: PEPOTensor, B <: PEPOTensor} = Tuple{K, B}
+ket(O::PEPOPurifiedSandwich) = O[1]
+bra(O::PEPOPurifiedSandwich) = O[2]
 
+const PEPOSandwich{N, T <: PEPSTensor, P <: PEPOTensor} = Tuple{T, T, Vararg{P, N}}
 ket(O::PEPOSandwich) = O[1]
 bra(O::PEPOSandwich) = O[2]
 pepo(O::PEPOSandwich) = O[3:end]
diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl
index 766457c23..e5d4d2306 100644
--- a/test/utility/correlator.jl
+++ b/test/utility/correlator.jl
@@ -17,9 +17,8 @@ function get_spaces(sym::Type{<:Sector})
 end
 
 site0 = CartesianIndex(1, 1)
-maxsep = 6
-site1xs = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep)
-site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep)
+site1xs = collect(site0 + CartesianIndex(0, i) for i in [1, -1, 3, -2])
+site1ys = collect(site0 + CartesianIndex(i, 0) for i in [1, -1, 3, -2])
 
 @testset "Correlator in InfinitePEPS ($(sym))" for sym in syms
     Random.seed!(100)
@@ -40,6 +39,31 @@ site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep)
             @info vals2
             @test vals1 ≈ vals2
         end
+        @test_throws ArgumentError correlator(peps, op, site0, site0, env)
+    end
+end
+
+@testset "Correlator in purified InfinitePEPO ($(sym))" for sym in syms
+    Random.seed!(100)
+    Vphy, Venv, Nspaces, Espaces = get_spaces(sym)
+    # TODO: test dual physical space
+    for Vp in [Vphy]
+        op = randn(ComplexF64, Vp ⊗ Vp → Vp ⊗ Vp)
+        Pspaces = fill(Vp, size(Nspaces))
+        pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces)
+        peps = InfinitePEPS(pepo)
+        env = CTMRGEnv(randn, ComplexF64, peps, Venv)
+        for site1s in (site1xs, site1ys)
+            vals1 = correlator(pepo, op, site0, site1s, pepo, env)
+            vals2 = map(site1s) do site1
+                O = LocalOperator(Pspaces, (site0, site1) => op)
+                return expectation_value(pepo, O, pepo, env)
+            end
+            @info vals1
+            @info vals2
+            @test vals1 ≈ vals2
+        end
+        @test_throws ArgumentError correlator(pepo, op, site0, site0, pepo, env)
     end
 end
 
@@ -63,5 +87,6 @@ end
             @info vals2
             @test vals1 ≈ vals2
         end
+        @test_throws ArgumentError correlator(pepo, op, site0, site0, env)
     end
 end