Add correlator for InfinitePEPO (mixed state)

src/PEPSKit.jl

@@ -15,7 +15,7 @@ using ChainRulesCore, Zygote
 using LoggingExtras
 
 using MPSKit
-using MPSKit: MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix
+using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, ProductTransferMatrix
 import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 import MPSKit: infinite_temperature_density_matrix
 
@@ -52,12 +52,15 @@ include("environments/vumps_environments.jl")
 include("environments/suweight.jl")
 
 include("algorithms/contractions/ctmrg_contractions.jl")
+include("algorithms/contractions/transfer.jl")
 include("algorithms/contractions/localoperator.jl")
 include("algorithms/contractions/vumps_contractions.jl")
 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")
@@ -74,6 +77,7 @@ include("algorithms/time_evolution/evoltools.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
 
+include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
 include("algorithms/correlators.jl")
 

src/algorithms/contractions/correlator/pepo_1layer.jl

@@ -0,0 +1,56 @@
+function start_correlator(
+        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 = 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[χSE De dstring; χ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::CTMRGEdgeTensor{T, S, 3},
+        ρ::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 = 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[χSW DW dstring; χ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::CTMRGEdgeTensor{T, S, 2}, 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,80 @@
+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/transfer.jl

@@ -0,0 +1,178 @@
+#
+# Transfer function for (CTMRG) edges
+#
+
+edge_transfer_left(v, ::Nothing, A, B) = edge_transfer_left(v, A, B)
+edge_transfer_right(v, ::Nothing, A, B) = edge_transfer_right(v, A, B)
+
+"""
+    edge_transfer_left(v, Et, Eb)
+
+Apply an edge transfer matrix to the left.
+
+```
+ ┌─Et─
+-v │
+ └─qƎ─
+```
+"""
+@generated function edge_transfer_left(
+        v::AbstractTensorMap{<:Any, S, 1, N₁},
+        Etop::CTMRGEdgeTensor{<:Any, S, N₂},
+        Ebot::CTMRGEdgeTensor{<:Any, S, N₂}
+    ) where {S, N₁, N₂}
+    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))
+    return macroexpand(
+        @__MODULE__, :(return @tensor $t_out := $t_in * $t_top * $t_bot)
+    )
+end
+
+
+"""
+    edge_transfer_right(v, Et, Eb)
+
+Apply an edge transfer matrix to the right.
+
+```
+─Et─┐
+ │  v-
+─qƎ─┘
+```
+"""
+@generated function edge_transfer_right(
+        v::AbstractTensorMap{<:Any, S, 1, N₁},
+        Etop::CTMRGEdgeTensor{<:Any, S, N₂},
+        Ebot::CTMRGEdgeTensor{<:Any, S, N₂}
+    ) where {S, N₁, N₂}
+    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))
+    return macroexpand(
+        @__MODULE__, :(return @tensor $t_out := $t_top * $t_bot * $t_in)
+    )
+end
+
+"""
+    edge_transfer_left(v, O, Et, Eb)
+
+Apply an edge transfer matrix to the left.
+
+```
+ ┌──Et─
+ │  │
+ v──O──
+ │  │
+ └──qƎ─
+```
+"""
+function edge_transfer_left(
+        v::CTMRGEdgeTensor{<:Any, S, 3},
+        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] :=
+        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},
+    ) where {S}
+    @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
+
+"""
+    transfer_right(v, Et, Eb)
+
+Apply an edge transfer matrix to the right.
+
+```
+──Et─┐
+  │  │
+──O──v
+  │  │
+──qƎ─┘
+```
+"""
+function edge_transfer_right(
+        v::CTMRGEdgeTensor{<:Any, S, 3},
+        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] :=
+        v[χ_NE D_E_above D_E_below; χ_SE] *
+        Etop[χ_NW D_N_above D_N_below; χ_NE] *
+        Ebot[χ_SE D_S_below D_S_above; χ_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},
+    ) 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
+
+"""
+    edge_transfer_left(v, O, Et, Eb)
+
+Apply an edge transfer matrix to the left on an excited vector.
+
+```
+ ┌──Et─
+ │  │
+-v──O──
+ │  │
+ └──qƎ─
+```
+"""
+function edge_transfer_left(
+        v::CTMRGEdgeTensor{<:Any, S, 4}, 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] :=
+        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(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_left(
+        v::CTMRGEdgeTensor{<:Any, S, 3}, O::PFTensor,
+        Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2},
+    ) where {S}
+    return @autoopt @tensor v′[χ_SE D_E d_string; χ_NE] :=
+        v[χ_SW D_W d_string; χ_NW] *
+        Etop[χ_NW D_N; χ_NE] *
+        Ebot[χ_SE D_S; χ_SW] *
+        O[D_W D_S; D_N D_E]
+end