Add correlator for InfinitePEPO (mixed state)

src/PEPSKit.jl

@@ -14,7 +14,7 @@ using ChainRulesCore, Zygote
 using LoggingExtras
 
 using MPSKit
-using MPSKit: MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix
+using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix
 import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 import MPSKit: infinite_temperature_density_matrix
 
@@ -57,6 +57,8 @@ include("algorithms/contractions/bondenv/benv_tools.jl")
 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_1layer.jl")
 
 include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")

src/algorithms/contractions/correlator/pepo_1layer.jl

@@ -0,0 +1,70 @@
+function MPSKit.transfer_left(
+        vec::AbstractTensor{T, S, 3}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}
+    ) where {T, S}
+    return @tensor y[-1 -2 -3] := vec[1 2 4] *
+        A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1])
+end
+
+function MPSKit.transfer_left(
+        v::AbstractTensor{T, S, 4}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}
+    ) where {T, S}
+    return @tensor t[d_string -1 -2 -3] := v[d_string 1 2 4] *
+        A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1])
+end
+
+function start_correlator(
+        i::CartesianIndex{2}, ρ::InfinitePEPO,
+        O::MPOTensor, env::CTMRGEnv
+    )
+    (size(ρ, 3) == 1) ||
+        throw(ArgumentError("The input PEPO ρ must have only one layer."))
+    r, c = Tuple(i)
+    E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)]
+    E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)]
+    E_west = env.edges[WEST, mod1(r, end), _prev(c, end)]
+    C_northwest = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)]
+    C_southwest = env.corners[SOUTHWEST, _next(r, end), _prev(c, end)]
+    t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2)
+    # TODO: part of these contractions is duplicated between the two output tensors,
+    # so could be optimized further
+    @autoopt @tensor Vn[χSE De χNE] :=
+        E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] *
+        E_west[χSW Dw; χNW] * C_northwest[χNW; χN] *
+        t[d d; Dn De Ds Dw] * E_north[χN Dn; χNE]
+    @autoopt @tensor Vo[dstring χSE De χNE] :=
+        E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] *
+        E_west[χSW Dw; χNW] * C_northwest[χNW; χN] *
+        removeunit(O, 1)[d2; d1 dstring] *
+        t[d1 d2; Dn De Ds Dw] * E_north[χN Dn; χNE]
+    return Vn, Vo
+end
+
+function end_correlator_numerator(
+        j::CartesianIndex{2}, V::AbstractTensor{T, S, 4},
+        ρ::InfinitePEPO, O::MPOTensor, env::CTMRGEnv
+    ) where {T, S}
+    (size(ρ, 3) == 1) ||
+        throw(ArgumentError("The input PEPO ρ must have only one layer."))
+    r, c = Tuple(j)
+    E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)]
+    E_east = env.edges[EAST, mod1(r, end), _next(c, end)]
+    E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)]
+    C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)]
+    t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2)
+    return @autoopt @tensor V[dstring χSW DW χNW] *
+        E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] *
+        C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] *
+        t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1]
+end
+
+function end_correlator_denominator(
+        j::CartesianIndex{2}, V::AbstractTensor{T, S, 3}, env::CTMRGEnv
+    ) where {T, S}
+    r, c = Tuple(j)
+    C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)]
+    E_east = env.edges[EAST, mod1(r, end), _next(c, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)]
+    return @autoopt @tensor V[χS DE χN] * C_northeast[χN; χNE] *
+        E_east[χNE DE; χSE] * C_southeast[χSE; χS]
+end

src/algorithms/contractions/correlator/peps.jl

@@ -0,0 +1,93 @@
+# transfer with excited GL
+function MPSKit.transfer_left(
+        GL::GenericMPSTensor{S, 4}, O::PEPSSandwich,
+        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3},
+    ) where {S}
+    return @autoopt @tensor GL′[χ_SE D_E_above d_string D_E_below; χ_NE] :=
+        GL[χ_SW D_W_above d_string D_W_below; χ_NW] *
+        conj(Ā[χ_SW D_S_above D_S_below; χ_SE]) *
+        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]) *
+        A[χ_NW D_N_above D_N_below; χ_NE]
+end
+
+function start_correlator(
+        i::CartesianIndex{2},
+        below::InfinitePEPS,
+        O::MPOTensor,
+        above::InfinitePEPS,
+        env::CTMRGEnv,
+    )
+    r, c = Tuple(i)
+    E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)]
+    E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)]
+    E_west = env.edges[WEST, mod1(r, end), _prev(c, end)]
+    C_northwest = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)]
+    C_southwest = env.corners[SOUTHWEST, _next(r, end), _prev(c, end)]
+    sandwich = (below[mod1(r, end), mod1(c, end)], above[mod1(r, end), mod1(c, end)])
+
+    # TODO: part of these contractions is duplicated between the two output tensors,
+    # so could be optimized further
+    @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(sandwich)[d; Dnbot Debot Dsbot Dwbot]) *
+        ket(sandwich)[d; 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(sandwich)[d1; Dnbot Debot Dsbot Dwbot]) *
+        removeunit(O, 1)[d1; d2 dstring] *
+        ket(sandwich)[d2; Dntop Detop Dstop Dwtop] *
+        E_north[χN Dntop Dnbot; χNE]
+
+    return Vn, Vo
+end
+
+function end_correlator_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 = env.edges[NORTH, _prev(r, end), mod1(c, end)]
+    E_east = env.edges[EAST, mod1(r, end), _next(c, end)]
+    E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)]
+    C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)]
+    sandwich = (above[mod1(r, end), mod1(c, end)], below[mod1(r, end), mod1(c, end)])
+
+    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(sandwich)[db; DNb DEb DSb DWb]) *
+        ket(sandwich)[dt; DNt DEt DSt DWt] *
+        removeunit(O, 4)[dstring db; dt]
+end
+
+function end_correlator_denominator(
+        j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1},
+        env::CTMRGEnv
+    ) where {T, S}
+    r, c = Tuple(j)
+    C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)]
+    E_east = env.edges[EAST, mod1(r, end), _next(c, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)]
+
+    return @autoopt @tensor V[χS DEt DEb; χN] *
+        C_northeast[χN; χNE] *
+        E_east[χNE DEt DEb; χSE] *
+        C_southeast[χSE; χS]
+end

src/algorithms/contractions/vumps_contractions.jl

@@ -32,19 +32,6 @@ function _transfer_left(
         A[χ_NW D_N_above D_N_below; χ_NE]
 end
 
-# transfer with excited GL
-function MPSKit.transfer_left(
-        GL::GenericMPSTensor{S, 4}, O::PEPSSandwich,
-        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3},
-    ) where {S}
-    return @autoopt @tensor GL′[χ_SE D_E_above d_string D_E_below; χ_NE] :=
-        GL[χ_SW D_W_above d_string D_W_below; χ_NW] *
-        conj(Ā[χ_SW D_S_above D_S_below; χ_SE]) *
-        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]) *
-        A[χ_NW D_N_above D_N_below; χ_NE]
-end
-
 function _transfer_right(
         GR::GenericMPSTensor{S, 3}, O::PEPSSandwich,
         A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3},