diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b3730aa46..0b7a817a2 100644 --- a/src/PEPSKit.jl +++ b/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") diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl new file mode 100644 index 000000000..a8e7f0490 --- /dev/null +++ b/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 diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl new file mode 100644 index 000000000..697ad2a72 --- /dev/null +++ b/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 diff --git a/src/algorithms/contractions/transfer.jl b/src/algorithms/contractions/transfer.jl new file mode 100644 index 000000000..3743396ed --- /dev/null +++ b/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 diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 97f231905..58584f7e6 100644 --- a/src/algorithms/contractions/vumps_contractions.jl +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -7,7 +7,7 @@ function MPSKit.transfer_left( A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N} Ā = twistdual(Ā, 2:N) - return _transfer_left(GL, O, A, Ā) + return mps_transfer_left(GL, O, A, Ā) end function MPSKit.transfer_right( @@ -15,12 +15,12 @@ function MPSKit.transfer_right( A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N} Ā = twistdual(Ā, 2:N) - return _transfer_right(GR, O, A, Ā) + return mps_transfer_right(GR, O, A, Ā) end ## PEPS -function _transfer_left( +function mps_transfer_left( GL::GenericMPSTensor{S, 3}, O::PEPSSandwich, A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, ) where {S} @@ -32,20 +32,7 @@ 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( +function mps_transfer_right( GR::GenericMPSTensor{S, 3}, O::PEPSSandwich, A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, ) where {S} @@ -59,7 +46,7 @@ end ## PEPO -@generated function _transfer_left( +@generated function mps_transfer_left( GL::GenericMPSTensor{S, N}, O::PEPOSandwich{H}, A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N, H} @@ -83,7 +70,7 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $GL´_e := $rhs)) end -@generated function _transfer_right( +@generated function mps_transfer_right( GR::GenericMPSTensor{S, N}, O::PEPOSandwich{H}, A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N, H} diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index f97d03b84..674644aab 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -1,3 +1,7 @@ +const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianIndices{N}} + +# Correlators in InfinitePEPS + function correlator_horizontal( bra::InfinitePEPS, operator, @@ -11,10 +15,8 @@ function correlator_horizontal( 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, @@ -22,148 +24,59 @@ function correlator_horizontal( 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 - Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] sandwich = ( ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)], ) - T = TransferMatrix(Atop, sandwich, _dag(Abot)) + T = edge_transfermatrix(Etop, sandwich, Ebot) Vo = Vo * T - twistdual!(T.below, 2:numout(T.below)) 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 - Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] sandwich = ( ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)], ) - T = TransferMatrix(Atop, sandwich, _dag(Abot)) - Vo = Vo * T - twistdual!(T.below, 2:numout(T.below)) + 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 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, - above::InfinitePEPS, - O::MPOTensor, - below::InfinitePEPS, - env::CTMRGEnv, - ) - 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, env::CTMRGEnv) - 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 - function correlator_vertical( bra::InfinitePEPS, - O, + 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_j = map(j -> siterotl90(j, size(bra)), js) - + rotated_js = map(j -> siterotl90(j, size(bra)), js) return correlator_horizontal( - rotated_bra, O, rotated_i, rotated_j, rotated_ket, rotl90(env) + rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env) ) end -const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianIndices{N}} - function MPSKit.correlator( bra::InfinitePEPS, O, @@ -195,3 +108,97 @@ end 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 = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + 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 = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + 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 +end + +function MPSKit.correlator( + ρ::InfinitePEPO, O, + i::CartesianIndex{2}, j::CartesianIndex{2}, + env::CTMRGEnv, + ) + return only(correlator(ρ, O, i, j:j, env)) +end + +# TODO: Correlators in InfinitePEPO (⟨ρ|O|ρ⟩) diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index 3f2c26a9d..5cfa17106 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -38,7 +38,7 @@ function gauge_fix(envprev::CTMRGEnv{C, T}, envfinal::CTMRGEnv{C, T}) where {C, # Find right fixed points of mixed transfer matrices ρinit = randn( - scalartype(T), MPSKit._lastspace(Tsfinal[end])' ← MPSKit._lastspace(M[end])' + scalartype(T), space(Tsfinal[end], numind(Tsfinal[end]))' ← space(M[end], numind(M[end]))' ) ρprev = transfermatrix_fixedpoint(Tsprev, M, ρinit) ρfinal = transfermatrix_fixedpoint(Tsfinal, M, ρinit) @@ -56,7 +56,10 @@ end # this is a bit of a hack to get the fixed point of the mixed transfer matrix # because MPSKit is not compatible with AD -@generated function _transfer_right( +# NOTE: the action of the transfer operator here is NOT the same as that of +# MPSKit.transfer_right; to match the MPSKit implementation, it would require a twist on any +# index where spaces(Abar, 2:end) is dual (which we can just ignore as long as we use a random Abar) +@generated function mps_transfer_right( v::AbstractTensorMap{<:Any, S, 1, N₁}, A::GenericMPSTensor{S, N₂}, Abar::GenericMPSTensor{S, N₂}, ) where {S, N₁, N₂} @@ -71,7 +74,7 @@ end function transfermatrix_fixedpoint(tops, bottoms, ρinit) _, vecs, info = eigsolve(ρinit, 1, :LM, Arnoldi()) do ρ return foldr(zip(tops, bottoms); init = ρ) do (top, bottom), ρ - return _transfer_right(ρ, top, bottom) + return mps_transfer_right(ρ, top, bottom) end end info.converged > 0 || @warn "eigsolve did not converge" diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index cdd2adb6d..13d4652ae 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -240,14 +240,36 @@ end end """ -Adjoint of an MPS tensor, but permutes the physical spaces back into the codomain. -Intuitively, this conjugates a tensor and then reinterprets its 'direction' as an MPS -tensor. + edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, + sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} + +Calculate the partial spectrum of the left edge transfer matrix corresponding to the given +`top` vector of edges and a `bot` vector of edge. The `sector` keyword argument can be used +to specify a non-trivial total charge for the transfer matrix eigenvectors. Specifically, an +auxiliary space `ℂ[typeof(sector)](sector => 1)'` will be added to the domain of each +eigenvector. The `tol` and `num_vals` keyword arguments are passed to `KrylovKit.eigolve`. """ -function _dag(A::MPSKit.GenericMPSTensor{S, N}) where {S, N} - return permute(A', ((1, (3:(N + 1))...), (2,))) +function edge_transfer_spectrum( + top::Vector{E}, bot::Vector{E}; tol = MPSKit.Defaults.tol, num_vals = 20, + sector = one(sectortype(E)) + ) where {E <: CTMRGEdgeTensor} + init = randn( + scalartype(E), + space(first(bot), numind(first(bot)))' ← ℂ[typeof(sector)](sector => 1)' ⊗ space(first(top), 1), + ) + + transferspace = fuse(space(first(top), 1) * space(first(bot), numind(first(bot)))') + num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values + eigenvals, eigenvecs, convhist = eigsolve( + flip(edge_transfermatrix(top, bot)), init, num_vals, :LM; tol = tol + ) + convhist.converged < num_vals && + @warn "correlation length failed to converge: normres = $(convhist.normres)" + + return eigenvals end + # TODO: decide on appropriate signature and returns for the more generic case """ correlation_length(state, env::CTMRGEnv; num_vals=2, kwargs...) @@ -271,15 +293,15 @@ function _correlation_length( # Horizontal λ_h = map(1:n_rows) do r - above = InfiniteMPS(env.edges[NORTH, r, :]) - below = InfiniteMPS(_dag.(env.edges[SOUTH, _next(r, n_rows), :])) - vals = MPSKit.transfer_spectrum(above; below, num_vals, sector, kwargs...) + top = env.edges[NORTH, r, :] + bot = env.edges[SOUTH, _next(r, n_rows), :] + vals = edge_transfer_spectrum(top, bot; num_vals, sector, kwargs...) # normalize using largest eigenvalue in trivial sector if isone(sector) N = first(vals) else - vals_triv = MPSKit.transfer_spectrum(above; below, num_vals = 1, kwargs...) + vals_triv = edge_transfer_spectrum(top, bot; num_vals = 1, kwargs...) N = first(vals_triv) end return vals ./ N # normalize largest eigenvalue @@ -287,15 +309,15 @@ function _correlation_length( # Vertical λ_v = map(1:n_cols) do c - above = InfiniteMPS(env.edges[EAST, :, c]) - below = InfiniteMPS(_dag.(env.edges[WEST, :, _next(c, n_cols)])) - vals = MPSKit.transfer_spectrum(above; below, num_vals, sector, kwargs...) + top = env.edges[EAST, :, c] + bot = env.edges[WEST, :, _next(c, n_cols)] + vals = edge_transfer_spectrum(top, bot; num_vals, sector, kwargs...) # normalize using largest eigenvalue in trivial sector if isone(sector) N = first(vals) else - vals_triv = MPSKit.transfer_spectrum(above; below, num_vals = 1, kwargs...) + vals_triv = edge_transfer_spectrum(top, bot; num_vals = 1, kwargs...) N = first(vals_triv) end return vals ./ N # normalize largest eigenvalue diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl new file mode 100644 index 000000000..354d3fe73 --- /dev/null +++ b/src/algorithms/transfermatrix.jl @@ -0,0 +1,38 @@ +# +# Edge transfer matrices +# + +# single site transfer +struct EdgeTransferMatrix{A <: CTMRGEdgeTensor, B, C <: CTMRGEdgeTensor} <: MPSKit.AbstractTransferMatrix + top::A + mid::B + bot::C + isflipped::Bool +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... +function TensorKit.flip(tm::EdgeTransferMatrix) + return EdgeTransferMatrix(tm.top, tm.mid, tm.bot, !tm.isflipped) +end + +# action on a vector using * is dispatched in MPSKit to regular function application +function (d::EdgeTransferMatrix)(vec) + return if d.isflipped + edge_transfer_left(vec, d.mid, d.top, d.bot) + else + edge_transfer_right(vec, d.mid, d.top, d.bot) + end +end + +# constructors +edge_transfermatrix(a) = edge_transfermatrix(a, nothing, a) +edge_transfermatrix(a, b) = edge_transfermatrix(a, nothing, b) +function edge_transfermatrix(a::CTMRGEdgeTensor, b, c::CTMRGEdgeTensor, isflipped = false) + return EdgeTransferMatrix(a, b, c, isflipped) +end +function edge_transfermatrix(a::AbstractVector, b, c::AbstractVector, isflipped = false) + tot = ProductTransferMatrix(convert(Vector, edge_transfermatrix.(a, b, c))) + return isflipped ? flip(tot) : tot +end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index e16b4cfca..e771cefb6 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -68,14 +68,17 @@ end function InfinitePEPO( Pspaces::A, Nspaces::A, Espaces::A = Nspaces ) where {A <: AbstractArray{<:ElementarySpace, 2}} + return InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) +end +function InfinitePEPO( + f, T, Pspaces::A, Nspaces::A, Espaces::A = Nspaces + ) where {A <: AbstractArray{<:ElementarySpace, 2}} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) - Pspaces = reshape(Pspaces, (size(Pspaces)..., 1)) Nspaces = reshape(Pspaces, (size(Nspaces)..., 1)) Espaces = reshape(Pspaces, (size(Espaces)..., 1)) - - return InfinitePEPO(Pspaces, Nspaces, Espaces) + return InfinitePEPO(f, T, Pspaces, Nspaces, Espaces) end """ diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl index 3ba9786ee..766457c23 100644 --- a/test/utility/correlator.jl +++ b/test/utility/correlator.jl @@ -2,27 +2,66 @@ 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) +const syms = (Z2Irrep, FermionParity) + +function get_spaces(sym::Type{<:Sector}) + @assert sym in syms + Nr, Nc = 2, 2 + Vphy = Vect[sym](0 => 1, 1 => 1) + V = Vect[sym](0 => 1, 1 => 2) + Venv = Vect[sym](0 => 2, 1 => 2) + Nspaces = [V' V; V V'] + Espaces = [V V'; V' V] + return Vphy, Venv, Nspaces, Espaces +end site0 = CartesianIndex(1, 1) -maxsep = 8 -site1s = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep) +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) -op = tJ.S_exchange(ComplexF64, Trivial, Trivial); +@testset "Correlator in InfinitePEPS ($(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)) + peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) + env = CTMRGEnv(randn, ComplexF64, peps, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(peps, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(peps, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 + end + end +end -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 +@testset "Correlator in 1-layer 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) + pf = InfinitePartitionFunction(pepo) + env = CTMRGEnv(randn, ComplexF64, pf, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(pepo, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(pepo, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 + end + end +end