From 1ee96413f17f865cb1b882c6b44a46e014fb790a Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:35:11 +0200 Subject: [PATCH 01/32] Slight reorganization --- src/PEPSKit.jl | 2 +- .../ctmrg/characteristic_equations.jl | 24 ---- .../contractions/ctmrg/contract_site.jl | 30 ---- .../contractions/ctmrg/network_value.jl | 131 ++++++++++++++++++ src/algorithms/toolbox.jl | 94 +------------ src/utility/util.jl | 60 ++++++++ 6 files changed, 198 insertions(+), 143 deletions(-) delete mode 100644 src/algorithms/contractions/ctmrg/contract_site.jl create mode 100644 src/algorithms/contractions/ctmrg/network_value.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index e6cae54f1..b6e0e626b 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -94,7 +94,7 @@ include("algorithms/contractions/ctmrg/halfinf_env.jl") include("algorithms/contractions/ctmrg/fullinf_env.jl") include("algorithms/contractions/ctmrg/renormalize_corner.jl") include("algorithms/contractions/ctmrg/renormalize_edge.jl") -include("algorithms/contractions/ctmrg/contract_site.jl") +include("algorithms/contractions/ctmrg/network_value.jl") include("algorithms/contractions/ctmrg/gaugefix.jl") include("algorithms/contractions/ctmrg/characteristic_equations.jl") diff --git a/src/algorithms/contractions/ctmrg/characteristic_equations.jl b/src/algorithms/contractions/ctmrg/characteristic_equations.jl index 477a039f9..740c087f5 100644 --- a/src/algorithms/contractions/ctmrg/characteristic_equations.jl +++ b/src/algorithms/contractions/ctmrg/characteristic_equations.jl @@ -349,30 +349,6 @@ function ChainRulesCore.rrule( return tout, twistnondual_pullback end -function absorb_left( - E::AbstractTensorMap{T, S}, C::CornerTensor{S} - ) where {T, S} - pC = (codomainind(C), domainind(C)) - pE = ((codomainind(E)[1],), (codomainind(E)[2:end]..., domainind(E)...)) - pCE = (codomainind(E), domainind(E)) - return tensorcontract(C, pC, false, E, pE, false, pCE) -end -function absorb_right( - P::AbstractTensorMap{T, S}, C::CornerTensor{S} - ) where {T, S} - pP = ((codomainind(P)..., domainind(P)[2:end]...), (domainind(P)[1],)) - pC = (codomainind(C), domainind(C)) - pPC = (codomainind(P), (domainind(P)[end], domainind(P)[1:(end - 1)]...)) - return tensorcontract(P, pP, false, C, pC, false, pPC) -end -# specialized versions; TODO: probably remove, this is a terrible idea for fermionic tensors... -function absorb_right(E::EdgeTensor{S}, C::CornerTensor{S}) where {S} - return E * C -end -function absorb_left_right(T::AbstractTensorMap, CL::CornerTensor, CR::CornerTensor) - return absorb_right(absorb_left(T, CL), CR) -end - # Partial contractions # -------------------- diff --git a/src/algorithms/contractions/ctmrg/contract_site.jl b/src/algorithms/contractions/ctmrg/contract_site.jl deleted file mode 100644 index aa8487909..000000000 --- a/src/algorithms/contractions/ctmrg/contract_site.jl +++ /dev/null @@ -1,30 +0,0 @@ -## Site contraction -@generated function _contract_site( - C_northwest, C_northeast, C_southeast, C_southwest, - E_north::TE, E_east::TE, E_south::TE, E_west::TE, - O::PEPOSandwich{H}, - ) where {TE <: CTMRGEdgeTensor, H} - @assert numout(TE) == H + 3 - - C_northwest_e = _corner_expr(:C_northwest, :WNW, :NNW) - C_northeast_e = _corner_expr(:C_northeast, :NNE, :ENE) - C_southeast_e = _corner_expr(:C_southeast, :ESE, :SSE) - C_southwest_e = _corner_expr(:C_southwest, :SSW, :WSW) - - E_north_e = _pepo_edge_expr(:E_north, :NNW, :NNE, :N, H) - E_east_e = _pepo_edge_expr(:E_east, :ENE, :ESE, :E, H) - E_south_e = _pepo_edge_expr(:E_south, :SSE, :SSW, :S, H) - E_west_e = _pepo_edge_expr(:E_west, :WSW, :WNW, :W, H) - - ket_e, bra_e, pepo_es = _pepo_sandwich_expr(:O, H) - - rhs = Expr( - :call, :*, - C_northwest_e, C_northeast_e, C_southeast_e, C_southwest_e, - E_north_e, E_east_e, E_south_e, E_west_e, - ket_e, Expr(:call, :conj, bra_e), - pepo_es..., - ) - - return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) -end diff --git a/src/algorithms/contractions/ctmrg/network_value.jl b/src/algorithms/contractions/ctmrg/network_value.jl new file mode 100644 index 000000000..9430a0df2 --- /dev/null +++ b/src/algorithms/contractions/ctmrg/network_value.jl @@ -0,0 +1,131 @@ +## Network value contractions +# +# The contractions making up `network_value`: the single site contraction, which +# dispatches on the local sandwich making up the network, and the corner and edge +# contractions which normalize it. The latter only involve environment tensors and are +# therefore independent of the sandwich type. + +## Site contraction +function _contract_site( + C_northwest, C_northeast, C_southeast, C_southwest, + E_north::CTMRG_PEPS_EdgeTensor, E_east::CTMRG_PEPS_EdgeTensor, + E_south::CTMRG_PEPS_EdgeTensor, E_west::CTMRG_PEPS_EdgeTensor, + O::PEPSSandwich, + ) + return @autoopt @tensor E_west[χ_WSW D_W_above D_W_below; χ_WNW] * + C_northwest[χ_WNW; χ_NNW] * + E_north[χ_NNW D_N_above D_N_below; χ_NNE] * + C_northeast[χ_NNE; χ_ENE] * + E_east[χ_ENE D_E_above D_E_below; χ_ESE] * + C_southeast[χ_ESE; χ_SSE] * + E_south[χ_SSE D_S_above D_S_below; χ_SSW] * + C_southwest[χ_SSW; χ_WSW] * + 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 _contract_site( + C_northwest, C_northeast, C_southeast, C_southwest, + E_north::CTMRG_PF_EdgeTensor, E_east::CTMRG_PF_EdgeTensor, + E_south::CTMRG_PF_EdgeTensor, E_west::CTMRG_PF_EdgeTensor, + O::PFTensor, + ) + return @autoopt @tensor E_west[χ_WSW D_W; χ_WNW] * + C_northwest[χ_WNW; χ_NNW] * + E_north[χ_NNW D_N; χ_NNE] * + C_northeast[χ_NNE; χ_ENE] * + E_east[χ_ENE D_E; χ_ESE] * + C_southeast[χ_ESE; χ_SSE] * + E_south[χ_SSE D_S; χ_SSW] * + C_southwest[χ_SSW; χ_WSW] * + O[D_W D_S; D_N D_E] +end + +@generated function _contract_site( + C_northwest, C_northeast, C_southeast, C_southwest, + E_north::TE, E_east::TE, E_south::TE, E_west::TE, + O::PEPOSandwich{H}, + ) where {TE <: CTMRGEdgeTensor, H} + @assert numout(TE) == H + 3 + + C_northwest_e = _corner_expr(:C_northwest, :WNW, :NNW) + C_northeast_e = _corner_expr(:C_northeast, :NNE, :ENE) + C_southeast_e = _corner_expr(:C_southeast, :ESE, :SSE) + C_southwest_e = _corner_expr(:C_southwest, :SSW, :WSW) + + E_north_e = _pepo_edge_expr(:E_north, :NNW, :NNE, :N, H) + E_east_e = _pepo_edge_expr(:E_east, :ENE, :ESE, :E, H) + E_south_e = _pepo_edge_expr(:E_south, :SSE, :SSW, :S, H) + E_west_e = _pepo_edge_expr(:E_west, :WSW, :WNW, :W, H) + + ket_e, bra_e, pepo_es = _pepo_sandwich_expr(:O, H) + + rhs = Expr( + :call, :*, + C_northwest_e, C_northeast_e, C_southeast_e, C_southwest_e, + E_north_e, E_east_e, E_south_e, E_west_e, + ket_e, Expr(:call, :conj, bra_e), + pepo_es..., + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) +end + +## Normalization contractions +function _contract_corners( + C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, + C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, + ) + return @tensor C_northwest[1; 2] * C_northeast[2; 3] * + C_southeast[3; 4] * C_southwest[4; 1] +end + +@generated function _contract_vertical_edges( + C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, + C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, + E_east::CTMRGEdgeTensor{T, S, N}, + E_west::CTMRGEdgeTensor{T, S, N}, + ) where {T, S, N} + C_northwest_e = tensorexpr(:C_northwest, (envlabel(:NW),), (envlabel(:N),)) + C_northeast_e = tensorexpr(:C_northeast, (envlabel(:N),), (envlabel(:NE),)) + C_southeast_e = tensorexpr(:C_southeast, (envlabel(:SE),), (envlabel(:S),)) + C_southwest_e = tensorexpr(:C_southwest, (envlabel(:S),), (envlabel(:SW),)) + + E_east_e = tensorexpr( + :E_east, (envlabel(:NE), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:SE),) + ) + E_west_e = tensorexpr( + :E_west, (envlabel(:SW), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:NW),) + ) + + rhs = Expr( + :call, :*, + E_west_e, C_northwest_e, C_northeast_e, E_east_e, C_southeast_e, C_southwest_e, + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) +end + +@generated function _contract_horizontal_edges( + C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, + C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, + E_north::CTMRGEdgeTensor{T, S, N}, E_south::CTMRGEdgeTensor{T, S, N}, + ) where {T, S, N} + C_northwest_e = tensorexpr(:C_northwest, (envlabel(:W),), (envlabel(:NW),)) + C_northeast_e = tensorexpr(:C_northeast, (envlabel(:NE),), (envlabel(:E),)) + C_southeast_e = tensorexpr(:C_southeast, (envlabel(:E),), (envlabel(:SE),)) + C_southwest_e = tensorexpr(:C_southwest, (envlabel(:SW),), (envlabel(:W),)) + + E_north_e = tensorexpr( + :E_north, (envlabel(:NW), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:NE),) + ) + E_south_e = tensorexpr( + :E_south, (envlabel(:SE), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:SW),) + ) + + rhs = Expr( + :call, :*, + C_northwest_e, E_north_e, C_northeast_e, C_southeast_e, E_south_e, C_southwest_e, + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) +end diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index b53e0f443..e30af47fa 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -104,39 +104,6 @@ function _contract_site(ind::Tuple{Int, Int}, network, env::CTMRGEnv) network[r, c], ) end -function _contract_site( - C_northwest, C_northeast, C_southeast, C_southwest, - E_north::CTMRG_PEPS_EdgeTensor, E_east::CTMRG_PEPS_EdgeTensor, - E_south::CTMRG_PEPS_EdgeTensor, E_west::CTMRG_PEPS_EdgeTensor, - O::PEPSSandwich, - ) - return @autoopt @tensor E_west[χ_WSW D_W_above D_W_below; χ_WNW] * - C_northwest[χ_WNW; χ_NNW] * - E_north[χ_NNW D_N_above D_N_below; χ_NNE] * - C_northeast[χ_NNE; χ_ENE] * - E_east[χ_ENE D_E_above D_E_below; χ_ESE] * - C_southeast[χ_ESE; χ_SSE] * - E_south[χ_SSE D_S_above D_S_below; χ_SSW] * - C_southwest[χ_SSW; χ_WSW] * - 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 _contract_site( - C_northwest, C_northeast, C_southeast, C_southwest, - E_north::CTMRG_PF_EdgeTensor, E_east::CTMRG_PF_EdgeTensor, - E_south::CTMRG_PF_EdgeTensor, E_west::CTMRG_PF_EdgeTensor, - O::PFTensor, - ) - return @autoopt @tensor E_west[χ_WSW D_W; χ_WNW] * - C_northwest[χ_WNW; χ_NNW] * - E_north[χ_NNW D_N; χ_NNE] * - C_northeast[χ_NNE; χ_ENE] * - E_east[χ_ENE D_E; χ_ESE] * - C_southeast[χ_ESE; χ_SSE] * - E_south[χ_SSE D_S; χ_SSW] * - C_southwest[χ_SSW; χ_WSW] * - O[D_W D_S; D_N D_E] -end """ _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv) @@ -146,11 +113,12 @@ environment `env`. """ function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) r, c = ind - C_NW = corner(env, NORTHWEST, r - 1, c - 1) - C_NE = corner(env, NORTHEAST, r - 1, c) - C_SE = corner(env, SOUTHEAST, r, c) - C_SW = corner(env, SOUTHWEST, r, c - 1) - return @tensor C_NW[1; 2] * C_NE[2; 3] * C_SE[3; 4] * C_SW[4; 1] + return _contract_corners( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r, c), + corner(env, SOUTHWEST, r, c - 1), + ) end """ @@ -170,31 +138,6 @@ function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) edge(env, WEST, r, c - 1), ) end -@generated function _contract_vertical_edges( - C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, - C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, - E_east::CTMRGEdgeTensor{T, S, N}, - E_west::CTMRGEdgeTensor{T, S, N}, - ) where {T, S, N} - C_northwest_e = tensorexpr(:C_northwest, (envlabel(:NW),), (envlabel(:N),)) - C_northeast_e = tensorexpr(:C_northeast, (envlabel(:N),), (envlabel(:NE),)) - C_southeast_e = tensorexpr(:C_southeast, (envlabel(:SE),), (envlabel(:S),)) - C_southwest_e = tensorexpr(:C_southwest, (envlabel(:S),), (envlabel(:SW),)) - - E_east_e = tensorexpr( - :E_east, (envlabel(:NE), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:SE),) - ) - E_west_e = tensorexpr( - :E_west, (envlabel(:SW), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:NW),) - ) - - rhs = Expr( - :call, :*, - E_west_e, C_northwest_e, C_northeast_e, E_east_e, C_southeast_e, C_southwest_e, - ) - - return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) -end """ _contract_horizontal_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) @@ -213,31 +156,6 @@ function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) edge(env, SOUTH, r, c), ) end -@generated function _contract_horizontal_edges( - C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, - C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, - E_north::CTMRGEdgeTensor{T, S, N}, E_south::CTMRGEdgeTensor{T, S, N}, - ) where {T, S, N} - C_northwest_e = tensorexpr(:C_northwest, (envlabel(:W),), (envlabel(:NW),)) - C_northeast_e = tensorexpr(:C_northeast, (envlabel(:NE),), (envlabel(:E),)) - C_southeast_e = tensorexpr(:C_southeast, (envlabel(:E),), (envlabel(:SE),)) - C_southwest_e = tensorexpr(:C_southwest, (envlabel(:SW),), (envlabel(:W),)) - - E_north_e = tensorexpr( - :E_north, (envlabel(:NW), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:NE),) - ) - E_south_e = tensorexpr( - :E_south, (envlabel(:SE), ntuple(i -> virtuallabel(i), N - 1)...), (envlabel(:SW),) - ) - - rhs = Expr( - :call, :*, - C_northwest_e, E_north_e, C_northeast_e, C_southeast_e, E_south_e, C_southwest_e, - ) - - return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) -end - """ edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} diff --git a/src/utility/util.jl b/src/utility/util.jl index 87e3a65f4..5b10c6ecc 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -62,6 +62,66 @@ function absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorM return U * sqrt_S, sqrt_S * V end +""" + absorb_left( + E::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb a matrix `C` into the left of a tensor `E` by contracting the first leg of `E` +with the second leg of `C`. +""" +function absorb_left( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + pC = (codomainind(C), domainind(C)) + pA = ((codomainind(A)[1],), (codomainind(A)[2:end]..., domainind(A)...)) + pCA = (codomainind(A), domainind(A)) + return tensorcontract(C, pC, false, A, pA, false, pCA) +end +function absorb_left( + P::AbstractTensorMap{<:Any, S, 1, N}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S, N} + return twistnondual(C, 2) * P +end + +""" + absorb_right( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb a matrix `C` into the right of a tensor `A` by contracting the last leg of `A` +with the first leg of `C`. +""" +function absorb_right( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + pA = ((codomainind(A)..., domainind(A)[2:end]...), (domainind(A)[1],)) + pC = (codomainind(C), domainind(C)) + pAC = (codomainind(A), (domainind(A)[end], domainind(A)[1:(end - 1)]...)) + return tensorcontract(A, pA, false, C, pC, false, pAC) +end +function absorb_right( + E::AbstractTensorMap{<:Any, S, N, 1}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S, N} + return E * twistdual(C, 1) +end + +""" + absorb_left_right( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb a matrix `C` into the right of a tensor `A` by contracting the last leg of `A` +with the first leg of `C`. +""" +function absorb_left_right( + T::AbstractTensorMap{<:Any, S}, + CL::AbstractTensorMap{<:Any, S, 1, 1}, + CR::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + return absorb_right(absorb_left(T, CL), CR) +end + _fliptwist_s(s::DiagonalTensorMap) = twist!(DiagonalTensorMap(flip(s, 1:2)), 1) """ From 01f20b910a826c9f6df465f2ce1923497c403638 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:40:11 +0200 Subject: [PATCH 02/32] Remove environment type restriction in expectation value --- src/algorithms/toolbox.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index e30af47fa..4bdc08c44 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -1,6 +1,6 @@ """ - expectation_value(state, O::LocalOperator, env::CTMRGEnv) - expectation_value(bra, O::LocalOperator, ket, env::CTMRGEnv) + expectation_value(state, O::LocalOperator, env) + expectation_value(bra, O::LocalOperator, ket, env) Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ of a [`LocalOperator`](@ref) `O`. This can be done either for a PEPS, or alternatively for a density matrix PEPO. @@ -8,7 +8,7 @@ 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::S, O::LocalOperator, ket::S, env::CTMRGEnv + bra::S, O::LocalOperator, ket::S, env ) where {S <: InfiniteState} checklattice(bra, O, ket) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly @@ -17,10 +17,8 @@ function MPSKit.expectation_value( end return sum(term_vals) end -MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::CTMRGEnv) = expectation_value(peps, O, peps, env) -function MPSKit.expectation_value( - state::InfinitePEPO, O::LocalOperator, env::CTMRGEnv - ) +MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) +function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) checklattice(state, O) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly ρ = reduced_densitymatrix(inds, state, env) From f0da3aff86bb82a0d1a93ba4152e2b7b6912a2bf Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:42:35 +0200 Subject: [PATCH 03/32] Reuse corner absorption --- src/algorithms/contractions/localoperator.jl | 72 ++++++++++---------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 56364cdde..d07a4f9e0 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -399,18 +399,18 @@ function reduced_densitymatrix1x1( A = ket[row, col] Ā = bra[row, col] - E_north = - edge(env, NORTH, row - 1, col) * - twistdual(corner(env, NORTHEAST, row - 1, col + 1), 1) - E_east = - edge(env, EAST, row, col + 1) * - twistdual(corner(env, SOUTHEAST, row + 1, col + 1), 1) - E_south = - edge(env, SOUTH, row + 1, col) * - twistdual(corner(env, SOUTHWEST, row + 1, col - 1), 1) - E_west = - edge(env, WEST, row, col - 1) * - twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) + E_north = absorb_right( + edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) + ) + E_east = absorb_right( + edge(env, EAST, row, col + 1), corner(env, SOUTHEAST, row + 1, col + 1) + ) + E_south = absorb_right( + edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) + ) + E_west = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) @tensor EE_SW[χSE χNW DSb DWb; DSt DWt] := E_south[χSE DSt DSb; χSW] * E_west[χSW DWt DWb; χNW] @@ -455,20 +455,20 @@ function reduced_densitymatrix2x1( A_south = ket[row + 1, col] Ā_south = bra[row + 1, col] - E_north = - edge(env, NORTH, row - 1, col) * - twistdual(corner(env, NORTHEAST, row - 1, col + 1), 1) + E_north = absorb_right( + edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) + ) E_northeast = edge(env, EAST, row, col + 1) - E_southeast = - edge(env, EAST, row + 1, col + 1) * - twistdual(corner(env, SOUTHEAST, row + 2, col + 1), 1) - E_south = - edge(env, SOUTH, row + 2, col) * - twistdual(corner(env, SOUTHWEST, row + 2, col - 1), 1) + E_southeast = absorb_right( + edge(env, EAST, row + 1, col + 1), corner(env, SOUTHEAST, row + 2, col + 1) + ) + E_south = absorb_right( + edge(env, SOUTH, row + 2, col), corner(env, SOUTHWEST, row + 2, col - 1) + ) E_southwest = edge(env, WEST, row + 1, col - 1) - E_northwest = - edge(env, WEST, row, col - 1) * - twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) + E_northwest = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) @tensor EE_NW[χW χNE DNWt DNt; DNWb DNb] := E_northwest[χW DNWt DNWb; χNW] * E_north[χNW DNt DNb; χNE] @@ -506,19 +506,19 @@ function reduced_densitymatrix1x2( Ā_east = bra[row, col + 1] E_northwest = edge(env, NORTH, row - 1, col) - E_northeast = - edge(env, NORTH, row - 1, col + 1) * - twistdual(corner(env, NORTHEAST, row - 1, col + 2), 1) - E_east = - edge(env, EAST, row, col + 2) * - twistdual(corner(env, SOUTHEAST, row + 1, col + 2), 1) + E_northeast = absorb_right( + edge(env, NORTH, row - 1, col + 1), corner(env, NORTHEAST, row - 1, col + 2) + ) + E_east = absorb_right( + edge(env, EAST, row, col + 2), corner(env, SOUTHEAST, row + 1, col + 2) + ) E_southeast = edge(env, SOUTH, row + 1, col + 1) - E_southwest = - edge(env, SOUTH, row + 1, col) * - twistdual(corner(env, SOUTHWEST, row + 1, col - 1), 1) - E_west = - edge(env, WEST, row, col - 1) * - twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) + E_southwest = absorb_right( + edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) + ) + E_west = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) @tensor EE_SW[χS χNW DSWt DWt; DSWb DWb] := E_southwest[χS DSWt DSWb; χSW] * E_west[χSW DWt DWb; χNW] From 5ba32feb76c06ed47b88540fc27c32136bb55de6 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:43:43 +0200 Subject: [PATCH 04/32] Split off local expectation value hook, and route dense TensorMap version through reduced density matrix --- src/algorithms/toolbox.jl | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 4bdc08c44..275ac881d 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -12,8 +12,7 @@ function MPSKit.expectation_value( ) 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) - return trmul(operator, ρ) + return local_expectation_value(inds, bra, operator, ket, env) end return sum(term_vals) end @@ -21,12 +20,33 @@ MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectatio function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) checklattice(state, O) term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - ρ = reduced_densitymatrix(inds, state, env) - return trmul(operator, ρ) + return local_expectation_value(inds, state, operator, env) end return sum(term_vals) end +""" + local_expectation_value(inds, bra, operator, ket, env) + local_expectation_value(inds, state, operator, env) + +Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation +value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state), where `operator` is the local +term acting on the sites `inds`. + +The implementation is overloaded based on the type of operator to be evaluated +""" +function local_expectation_value end + +# AbstractTensorMap evaluation goes through reduced density matrix +function local_expectation_value(inds, bra, operator::AbstractTensorMap, ket, env) + ρ = reduced_densitymatrix(inds, ket, bra, env) + return trmul(operator, ρ) +end +function local_expectation_value(inds, state, operator::AbstractTensorMap, env) + ρ = reduced_densitymatrix(inds, state, env) + return trmul(operator, ρ) +end + """ expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) From 35788f4164d71531ea892c018e0093392804d9de Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:44:57 +0200 Subject: [PATCH 05/32] Evaluate `BPEnv` expectation values in exactly the same way, and add reduced density operators --- .../contractions/bp_contractions.jl | 174 ++++++++---------- 1 file changed, 81 insertions(+), 93 deletions(-) diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index b751cbd75..dd1279909 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -48,61 +48,34 @@ absorb_south_message(A::PEPSTensor, M::PEPSMessage) = absorb_west_message(A::PEPSTensor, M::PEPSMessage) = @tensor A′[d; N E S W'] := A[d; N E S W] * M[W'; W] -# Belief Propagation Expectation values -# ------------------------------------- -function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::BPEnv) - checklattice(peps, O) - term_vals = dtmap([O.terms...]) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - contract_local_operator(inds, operator, peps, peps, env) / - contract_local_norm(inds, peps, peps, env) - end - return sum(term_vals) -end - -function contract_local_operator( +# Belief Propagation reduced density matrices +# ------------------------------------------- +# BP messages live on the bonds of the network, so there is no analogue of the generated +# rectangular-patch contraction used for a `CTMRGEnv`: only single sites and nearest +# neighbor pairs can be contracted without introducing loop corrections. +function reduced_densitymatrix( inds::Vector{CartesianIndex{2}}, - O::AbstractTensorMap, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv, ) - length(inds) == 1 && return contract_local_operator1x1(only(inds), O, ket, bra, env) + length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) if length(inds) == 2 ind_relative = inds[2] - inds[1] if ind_relative == CartesianIndex(1, 0) - return contract_local_operator2x1(inds[1], O, ket, bra, env) + return reduced_densitymatrix2x1(inds[1], ket, bra, env) elseif ind_relative == CartesianIndex(0, 1) - return contract_local_operator1x2(inds[1], O, ket, bra, env) + return reduced_densitymatrix1x2(inds[1], ket, bra, env) end end error("No implementation for contractions for BP environments with $inds") end -function contract_local_norm( - inds::Vector{CartesianIndex{2}}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, - ) - length(inds) == 1 && return contract_local_norm1x1(only(inds), ket, bra, env) +reduced_densitymatrix(inds, ket::InfinitePEPS, env::BPEnv) = + reduced_densitymatrix(inds, ket, ket, env) - if length(inds) == 2 - ind_relative = inds[2] - inds[1] - if ind_relative == CartesianIndex(1, 0) - return contract_local_norm2x1(inds[1], ket, bra, env) - elseif ind_relative == CartesianIndex(0, 1) - return contract_local_norm1x2(inds[1], ket, bra, env) - end - end - error("No implementation for contractions for BP environments with $inds") -end - -function contract_local_operator1x1( - ind::CartesianIndex{2}, - O::AbstractTensorMap{<:Any, <:Any, 1, 1}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, +function reduced_densitymatrix1x1( + ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) row, col = Tuple(ind) M_north = env[NORTH, row - 1, col] @@ -110,42 +83,19 @@ function contract_local_operator1x1( M_south = env[SOUTH, row + 1, col] M_west = env[WEST, row, col - 1] - return @autoopt @tensor begin + @autoopt @tensor ρ[dt; db] := ket[row, col][dt; DNt DEt DSt DWt] * - conj(bra[row, col][db; DNb DEb DSb DWb]) * - O[db; dt] * - M_north[DNt; DNb] * - M_east[DEt; DEb] * - M_south[DSb; DSt] * - M_west[DWb; DWt] - end -end - -function contract_local_norm1x1( - ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv - ) - row, col = Tuple(ind) - M_north = env[NORTH, row - 1, col] - M_east = env[EAST, row, col + 1] - M_south = env[SOUTH, row + 1, col] - M_west = env[WEST, row, col - 1] + conj(bra[row, col][db; DNb DEb DSb DWb]) * + M_north[DNt; DNb] * + M_east[DEt; DEb] * + M_south[DSb; DSt] * + M_west[DWb; DWt] - return @autoopt @tensor begin - ket[row, col][d; DNt DEt DSt DWt] * - conj(bra[row, col][d; DNb DEb DSb DWb]) * - M_north[DNt; DNb] * - M_east[DEt; DEb] * - M_south[DSb; DSt] * - M_west[DWb; DWt] - end + return ρ / str(ρ) end -function contract_local_operator2x1( - coord::CartesianIndex{2}, - O::AbstractTensorMap{<:Any, <:Any, 2, 2}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, +function reduced_densitymatrix2x1( + coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) row, col = Tuple(coord) M_north = env[NORTH, row - 1, col] @@ -155,7 +105,8 @@ function contract_local_operator2x1( M_southwest = env[WEST, row + 1, col - 1] M_northwest = env[WEST, row, col - 1] - return @autoopt @tensor ket[row, col][dNt; DNt DNEt DMt DNWt] * + @autoopt @tensor ρ[dNt dSt; dNb dSb] := + ket[row, col][dNt; DNt DNEt DMt DNWt] * ket[row + 1, col][dSt; DMt DSEt DSt DSWt] * conj(bra[row, col][dNb; DNb DNEb DMb DNWb]) * conj(bra[row + 1, col][dSb; DMb DSEb DSb DSWb]) * @@ -164,16 +115,13 @@ function contract_local_operator2x1( M_southeast[DSEt; DSEb] * M_south[DSb; DSt] * M_southwest[DSWb; DSWt] * - M_northwest[DNWb; DNWt] * - O[dNb dSb; dNt dSt] + M_northwest[DNWb; DNWt] + + return ρ / str(ρ) end -function contract_local_operator1x2( - coord::CartesianIndex{2}, - O::AbstractTensorMap{<:Any, <:Any, 2, 2}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, +function reduced_densitymatrix1x2( + coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) row, col = Tuple(coord) M_west = env[WEST, row, col - 1] @@ -183,22 +131,62 @@ function contract_local_operator1x2( M_southeast = env[SOUTH, row + 1, col + 1] M_southwest = env[SOUTH, row + 1, col] A_west = ket[row, col] - Ā_west = bra[row, col] + Ā_west = bra[row, col] A_east = ket[row, col + 1] - Ā_east = bra[row, col + 1] + Ā_east = bra[row, col + 1] - return @autoopt @tensor begin + @autoopt @tensor ρ[dWt dEt; dWb dEb] := A_west[dWt; DNWt DMt DSWt DWt] * - A_east[dEt; DNEt DEt DSEt DMt] * - conj(Ā_west[dWb; DNWb DMb DSWb DWb]) * - conj(Ā_east[dEb; DNEb DEb DSEb DMb]) * - M_west[DWb; DWt] * - M_northwest[DNWt; DNWb] * - M_northeast[DNEt; DNEb] * + A_east[dEt; DNEt DEt DSEt DMt] * + conj(Ā_west[dWb; DNWb DMb DSWb DWb]) * + conj(Ā_east[dEb; DNEb DEb DSEb DMb]) * + M_west[DWb; DWt] * + M_northwest[DNWt; DNWb] * + M_northeast[DNEt; DNEb] * + M_east[DEt; DEb] * + M_southeast[DSEb; DSEt] * + M_southwest[DSWb; DSWt] + + return ρ / str(ρ) +end + +# Local norm contractions +# ---------------------- +function contract_local_norm( + inds::Vector{CartesianIndex{2}}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::BPEnv, + ) + length(inds) == 1 && return contract_local_norm1x1(only(inds), ket, bra, env) + + if length(inds) == 2 + ind_relative = inds[2] - inds[1] + if ind_relative == CartesianIndex(1, 0) + return contract_local_norm2x1(inds[1], ket, bra, env) + elseif ind_relative == CartesianIndex(0, 1) + return contract_local_norm1x2(inds[1], ket, bra, env) + end + end + error("No implementation for contractions for BP environments with $inds") +end + +function contract_local_norm1x1( + ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv + ) + row, col = Tuple(ind) + M_north = env[NORTH, row - 1, col] + M_east = env[EAST, row, col + 1] + M_south = env[SOUTH, row + 1, col] + M_west = env[WEST, row, col - 1] + + return @autoopt @tensor begin + ket[row, col][d; DNt DEt DSt DWt] * + conj(bra[row, col][d; DNb DEb DSb DWb]) * + M_north[DNt; DNb] * M_east[DEt; DEb] * - M_southeast[DSEb; DSEt] * - M_southwest[DSWb; DSWt] * - O[dWb dEb; dWt dEt] + M_south[DSb; DSt] * + M_west[DWb; DWt] end end From e5facb89ebea90067d27b9511ca039ea434299a4 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 27 Aug 2026 14:46:29 +0200 Subject: [PATCH 06/32] Remove unused `contract_local_operator` --- src/algorithms/contractions/localoperator.jl | 69 +------------------- test/toolbox/densitymatrices.jl | 8 +-- 2 files changed, 5 insertions(+), 72 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index d07a4f9e0..5b432c159 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -7,35 +7,7 @@ envlabel(args...) = tensorlabel(:χ, args...) virtuallabel(args...) = tensorlabel(:D, args...) physicallabel(args...) = tensorlabel(:d, args...) -""" -$(SIGNATURES) - -Contract a local operator `O` on the PEPS `peps` at the indices `inds` using the environment `env`. - -This works by generating the appropriate contraction on a rectangular patch with its corners -specified by `inds`. The `peps` is contracted with `O` from above and below, and the PEPS-operator -sandwich is surrounded with the appropriate environment tensors. -""" -function contract_local_operator( - inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, - ) - static_inds = Tuple(Val.(inds)) - return _contract_local_operator(static_inds, O, (ket, bra), env) -end -function contract_local_operator( - inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, - ) - return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env) -end - -Base.@deprecate( - contract_local_operator(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv), - contract_local_operator(collect(inds), ket, bra, env) -) - -# This implements the contraction of an operator acting on sites `inds`. +# This implements the contractions on a rectangular patch spanned by the sites `inds`. # The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance). function _contract_corner_expr(rowrange, colrange) @@ -220,48 +192,13 @@ function _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds = end end -@generated function _contract_local_operator( - inds::NTuple{N, Val}, - O::AbstractTensorMap{T, S, N, N}, - state::Tuple{InfinitePEPS, InfinitePEPS}, - env::CTMRGEnv, - ) where {T, S, N} - cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val - allunique(cartesian_inds) || - throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) - rowrange = getindex.(cartesian_inds, 1) - colrange = getindex.(cartesian_inds, 2) - - corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) - edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, 2) - operator = tensorexpr( - :O, - ntuple(i -> physicallabel(:O, 2, i), N), - ntuple(i -> physicallabel(:O, 1, i), N), - ) - ket, bra = _contract_state_expr(rowrange, colrange, 2, cartesian_inds) - - multiplication_ex = Expr( - :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ket..., map(x -> Expr(:call, :conj, x), bra)..., - operator, - ) - - returnex = quote - @autoopt @tensor $multiplication_ex - end - return macroexpand(@__MODULE__, returnex) -end - """ $(SIGNATURES) Contract a local norm of the PEPS `peps` around indices `inds`. -This works analogously to [`contract_local_operator`](@ref) by generating the contraction -on a rectangular patch based on `inds` but replacing the operator with an identity such +This works by generating the appropriate contraction on a rectangular patch with its +corners specified by `inds`, contracting the `peps` with itself from above and below such that the PEPS norm is computed. (Note that this is not the physical norm of the state.) """ function contract_local_norm( diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl index 30067d3d1..9e2c29986 100644 --- a/test/toolbox/densitymatrices.jl +++ b/test/toolbox/densitymatrices.jl @@ -1,6 +1,6 @@ using TensorKit using PEPSKit -using PEPSKit: contract_local_operator, contract_local_norm +using PEPSKit: contract_local_norm, _contract_site, InfiniteSquareNetwork using Test using TestExtras @@ -57,9 +57,8 @@ end O_doubled_singlesite = LocalOperator(physicalspace(ρ_peps), (site,) => O_doubled) E2 = expectation_value(ρ_peps, O_doubled_singlesite, ρ_peps, env) @test E1 ≈ E2 - val = contract_local_operator([site], O_doubled, ρ_peps, ρ_peps, env) nrm = contract_local_norm([site], ρ_peps, ρ_peps, env) - @test E1 ≈ val / nrm + @test nrm ≈ _contract_site(site, InfiniteSquareNetwork(ρ_peps), env) # two sites for inds in zip( @@ -71,8 +70,5 @@ end O_doubled_twosite = LocalOperator(physicalspace(ρ_peps), inds => O_doubled ⊗ O_doubled) E2 = expectation_value(ρ_peps, O_doubled_twosite, ρ_peps, env) @test E1 ≈ E2 - val = contract_local_operator(collect(inds), O_doubled ⊗ O_doubled, ρ_peps, ρ_peps, env) - nrm = contract_local_norm(collect(inds), ρ_peps, ρ_peps, env) - @test E1 ≈ val / nrm end end From 496f063f6863dfa10a354299c22256f16e2e6dc9 Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 31 Aug 2026 13:22:18 +0200 Subject: [PATCH 07/32] Add vanity line break --- src/algorithms/toolbox.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 275ac881d..4ef530183 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -174,6 +174,7 @@ function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) edge(env, SOUTH, r, c), ) end + """ edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} From 463f48d6583b1bd7ec31195825aa3756adcafe2a Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 31 Aug 2026 15:34:39 +0200 Subject: [PATCH 08/32] Drop more environment type specifiers --- src/algorithms/toolbox.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 4ef530183..b44ca4f0f 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -61,24 +61,24 @@ convention. function MPSKit.expectation_value( pf::InfinitePartitionFunction, op::Pair{CartesianIndex{2}, <:AbstractTensorMap{T, S, 2, 2}}, - env::CTMRGEnv, + env, ) where {T, S} return contract_local_tensor(op[1], op[2], env) / contract_local_tensor(op[1], pf[op[1]], env) end function MPSKit.expectation_value( - pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env::CTMRGEnv + pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env ) return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) end """ - cost_function(peps::InfinitePEPS, env::CTMRGEnv, O::LocalOperator) + cost_function(peps::InfinitePEPS, env, O::LocalOperator) Real part of expectation value of `O`. Prints a warning if the expectation value yields a finite imaginary part (up to a tolerance). """ -function cost_function(peps::InfinitePEPS, env::CTMRGEnv, O::LocalOperator) +function cost_function(peps::InfinitePEPS, env, O::LocalOperator) E = MPSKit.expectation_value(peps, O, env) ignore_derivatives() do isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || From 00c33b6823100567aad8388907905c8eb4e75e2d Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 31 Aug 2026 16:28:30 +0200 Subject: [PATCH 09/32] Remove unused method (accidentally missed during deprecation) --- src/algorithms/contractions/localoperator.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 5b432c159..45615d140 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -366,19 +366,6 @@ function reduced_densitymatrix1x1( return ρ / str(ρ) end -function reduced_densitymatrix( - inds::NTuple{2, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) - if inds[2] - inds[1] == CartesianIndex(1, 0) - return reduced_densitymatrix2x1(inds[1], ket, bra, env) - elseif inds[2] - inds[1] == CartesianIndex(0, 1) - return reduced_densitymatrix1x2(inds[1], ket, bra, env) - else - static_inds = Val.(inds) - return _contract_densitymatrix(static_inds, (ket, bra), env) - end -end - # Special case 2x1 density matrix: # Keep contraction order but try to optimize intermediate permutations: function reduced_densitymatrix2x1( From e1cf86e0bb68a0568872829f466fca2684538bfb Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 31 Aug 2026 16:29:46 +0200 Subject: [PATCH 10/32] Drop more type specs --- src/algorithms/contractions/localoperator.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 45615d140..e329f7c51 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -253,7 +253,7 @@ specified by `inds`. The result is normalized such that `tr(ρ) = 1`. """ reduced_densitymatrix function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env ) length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) @@ -269,57 +269,57 @@ function reduced_densitymatrix( return _contract_densitymatrix(static_inds, (ket, bra), env) end function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv + inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env ) size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) static_inds = Tuple(Val.(inds)) return _contract_densitymatrix(static_inds, (state,), env) end function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env ) size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) static_inds = Tuple(Val.(inds)) return _contract_densitymatrix(static_inds, (ket, bra), env) end -reduced_densitymatrix(inds, ket::InfinitePEPS, env::CTMRGEnv) = +reduced_densitymatrix(inds, ket::InfinitePEPS, env) = reduced_densitymatrix(inds, ket, ket, env) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env ) where {N}, reduced_densitymatrix(collect(inds), ket, bra, env) ) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv + inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env ) where {N}, reduced_densitymatrix(collect(inds), state, env) ) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv + inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env ) where {N}, reduced_densitymatrix(collect(inds), ket, bra, env) ) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env ) where {N}, reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env) ) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv + inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPO, bra::InfinitePEPO, env ) where {N}, reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env) ) Base.@deprecate( reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env::CTMRGEnv + inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env ) where {N}, reduced_densitymatrix(collect(CartesianIndex.(inds)), state, env) ) From 15d2330297c1b18b6a2edc4edbc423b3f4e1dc1e Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 2 Sep 2026 09:06:11 +0200 Subject: [PATCH 11/32] Restore contract_local_operator for future use --- .../contractions/bp_contractions.jl | 94 ------------------- src/algorithms/contractions/localoperator.jl | 75 ++++++++++++++- test/toolbox/densitymatrices.jl | 7 +- 3 files changed, 77 insertions(+), 99 deletions(-) diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index dd1279909..557a4e232 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -149,97 +149,3 @@ function reduced_densitymatrix1x2( return ρ / str(ρ) end - -# Local norm contractions -# ---------------------- -function contract_local_norm( - inds::Vector{CartesianIndex{2}}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, - ) - length(inds) == 1 && return contract_local_norm1x1(only(inds), ket, bra, env) - - if length(inds) == 2 - ind_relative = inds[2] - inds[1] - if ind_relative == CartesianIndex(1, 0) - return contract_local_norm2x1(inds[1], ket, bra, env) - elseif ind_relative == CartesianIndex(0, 1) - return contract_local_norm1x2(inds[1], ket, bra, env) - end - end - error("No implementation for contractions for BP environments with $inds") -end - -function contract_local_norm1x1( - ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv - ) - row, col = Tuple(ind) - M_north = env[NORTH, row - 1, col] - M_east = env[EAST, row, col + 1] - M_south = env[SOUTH, row + 1, col] - M_west = env[WEST, row, col - 1] - - return @autoopt @tensor begin - ket[row, col][d; DNt DEt DSt DWt] * - conj(bra[row, col][d; DNb DEb DSb DWb]) * - M_north[DNt; DNb] * - M_east[DEt; DEb] * - M_south[DSb; DSt] * - M_west[DWb; DWt] - end -end - -function contract_local_norm2x1( - coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv - ) - row, col = Tuple(coord) - M_north = env[NORTH, row - 1, col] - M_northeast = env[EAST, row, col + 1] - M_southeast = env[EAST, row + 1, col + 1] - M_south = env[SOUTH, row + 2, col] - M_southwest = env[WEST, row + 1, col - 1] - M_northwest = env[WEST, row, col - 1] - - return @autoopt @tensor ket[row, col][dN; DNt DNEt DMt DNWt] * - ket[row + 1, col][dS; DMt DSEt DSt DSWt] * - conj(bra[row, col][dN; DNb DNEb DMb DNWb]) * - conj(bra[row + 1, col][dS; DMb DSEb DSb DSWb]) * - M_north[DNt; DNb] * - M_northeast[DNEt; DNEb] * - M_southeast[DSEt; DSEb] * - M_south[DSb; DSt] * - M_southwest[DSWb; DSWt] * - M_northwest[DNWb; DNWt] -end - -function contract_local_norm1x2( - coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv - ) - row, col = Tuple(coord) - - M_west = env[WEST, row, col - 1] - M_northwest = env[NORTH, row - 1, col] - M_northeast = env[NORTH, row - 1, col + 1] - M_east = env[EAST, row, col + 2] - M_southeast = env[SOUTH, row + 1, col + 1] - M_southwest = env[SOUTH, row + 1, col] - - A_west = ket[row, col] - Ā_west = bra[row, col] - A_east = ket[row, col + 1] - Ā_east = bra[row, col + 1] - - return @autoopt @tensor begin - A_west[dW; DNWt DMt DSWt DWt] * - A_east[dE; DNEt DEt DSEt DMt] * - conj(Ā_west[dW; DNWb DMb DSWb DWb]) * - conj(Ā_east[dE; DNEb DEb DSEb DMb]) * - M_west[DWb; DWt] * - M_northwest[DNWt; DNWb] * - M_northeast[DNEt; DNEb] * - M_east[DEt; DEb] * - M_southeast[DSEb; DSEt] * - M_southwest[DSWb; DSWt] - end -end diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index e329f7c51..7ea89e06d 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -7,7 +7,39 @@ envlabel(args...) = tensorlabel(:χ, args...) virtuallabel(args...) = tensorlabel(:D, args...) physicallabel(args...) = tensorlabel(:d, args...) -# This implements the contractions on a rectangular patch spanned by the sites `inds`. +""" +$(SIGNATURES) + +Contract a local operator `O` on the PEPS `ket`, `bra` at the indices `inds` using the +environment `env`. + +This works by generating the appropriate contraction on a rectangular patch with its corners +specified by `inds`. The `ket` and `bra` are contracted with `O` in between, and the +PEPS-operator sandwich is surrounded with the appropriate environment tensors. +""" +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap, + ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator( + inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap, + ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ) + return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env) +end + +Base.@deprecate( + contract_local_operator( + inds::NTuple, O::AbstractTensorMap, + ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ), + contract_local_operator(collect(inds), O, ket, bra, env) +) + +# This implements the contraction of an operator acting on sites `inds`. # The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance). function _contract_corner_expr(rowrange, colrange) @@ -192,13 +224,48 @@ function _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds = end end +@generated function _contract_local_operator( + inds::NTuple{N, Val}, + O::AbstractTensorMap{T, S, N, N}, + state::Tuple{InfinitePEPS, InfinitePEPS}, + env::CTMRGEnv, + ) where {T, S, N} + cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val + allunique(cartesian_inds) || + throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) + rowrange = getindex.(cartesian_inds, 1) + colrange = getindex.(cartesian_inds, 2) + + corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) + edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, 2) + operator = tensorexpr( + :O, + ntuple(i -> physicallabel(:O, 2, i), N), + ntuple(i -> physicallabel(:O, 1, i), N), + ) + ket, bra = _contract_state_expr(rowrange, colrange, 2, cartesian_inds) + + multiplication_ex = Expr( + :call, :*, + corner_NW, corner_NE, corner_SE, corner_SW, + edges_N..., edges_E..., edges_S..., edges_W..., + ket..., map(x -> Expr(:call, :conj, x), bra)..., + operator, + ) + + returnex = quote + @autoopt @tensor $multiplication_ex + end + return macroexpand(@__MODULE__, returnex) +end + """ $(SIGNATURES) -Contract a local norm of the PEPS `peps` around indices `inds`. +Contract a local norm of the PEPS `ket`, `bra` around indices `inds`. -This works by generating the appropriate contraction on a rectangular patch with its -corners specified by `inds`, contracting the `peps` with itself from above and below such +This works analogously to [`contract_local_operator`](@ref) by generating the contraction +on a rectangular patch based on `inds` but replacing the operator with an identity such that the PEPS norm is computed. (Note that this is not the physical norm of the state.) """ function contract_local_norm( diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl index 9e2c29986..803ab2283 100644 --- a/test/toolbox/densitymatrices.jl +++ b/test/toolbox/densitymatrices.jl @@ -1,6 +1,6 @@ using TensorKit using PEPSKit -using PEPSKit: contract_local_norm, _contract_site, InfiniteSquareNetwork +using PEPSKit: contract_local_operator, contract_local_norm, _contract_site, InfiniteSquareNetwork using Test using TestExtras @@ -57,8 +57,10 @@ end O_doubled_singlesite = LocalOperator(physicalspace(ρ_peps), (site,) => O_doubled) E2 = expectation_value(ρ_peps, O_doubled_singlesite, ρ_peps, env) @test E1 ≈ E2 + val = contract_local_operator([site], O_doubled, ρ_peps, ρ_peps, env) nrm = contract_local_norm([site], ρ_peps, ρ_peps, env) @test nrm ≈ _contract_site(site, InfiniteSquareNetwork(ρ_peps), env) + @test E1 ≈ val / nrm # two sites for inds in zip( @@ -70,5 +72,8 @@ end O_doubled_twosite = LocalOperator(physicalspace(ρ_peps), inds => O_doubled ⊗ O_doubled) E2 = expectation_value(ρ_peps, O_doubled_twosite, ρ_peps, env) @test E1 ≈ E2 + val = contract_local_operator(collect(inds), O_doubled ⊗ O_doubled, ρ_peps, ρ_peps, env) + nrm = contract_local_norm(collect(inds), ρ_peps, ρ_peps, env) + @test E1 ≈ val / nrm end end From f9b10f2ef683614f8f8b82e76c295a4aea45b03b Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 2 Sep 2026 10:04:13 +0200 Subject: [PATCH 12/32] More merging --- .../contractions/bp_contractions.jl | 35 +++++++++---------- src/algorithms/contractions/localoperator.jl | 12 +++---- 2 files changed, 22 insertions(+), 25 deletions(-) diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index 557a4e232..a9ae6494f 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -51,28 +51,25 @@ absorb_west_message(A::PEPSTensor, M::PEPSMessage) = # Belief Propagation reduced density matrices # ------------------------------------------- # BP messages live on the bonds of the network, so there is no analogue of the generated -# rectangular-patch contraction used for a `CTMRGEnv`: only single sites and nearest -# neighbor pairs can be contracted without introducing loop corrections. -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, - ket::InfinitePEPS, - bra::InfinitePEPS, - env::BPEnv, +# rectangular-patch contraction used for a `CTMRGEnv`: only single sites and nearest neighbor +# pairs can be contracted without introducing loop corrections. + +function _contract_densitymatrix(inds::NTuple{N, Val}, state::Tuple, env::BPEnv) where {N} + sites = map(v -> typeof(v).parameters[1], inds) + return throw( + ArgumentError( + "Cannot contract a $(_patch_shape_string(sites)) patch using a `BPEnv`; + only 1x1, 2x1, 1x2 patches are supported." + ) ) - length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) +end - if length(inds) == 2 - ind_relative = inds[2] - inds[1] - if ind_relative == CartesianIndex(1, 0) - return reduced_densitymatrix2x1(inds[1], ket, bra, env) - elseif ind_relative == CartesianIndex(0, 1) - return reduced_densitymatrix1x2(inds[1], ket, bra, env) - end - end - error("No implementation for contractions for BP environments with $inds") +function _patch_shape_string(sites) + rows, cols = getindex.(sites, 1), getindex.(sites, 2) + nrows = maximum(rows) - minimum(rows) + 1 + ncols = maximum(cols) - minimum(cols) + 1 + return "$(nrows)x$(ncols)" end -reduced_densitymatrix(inds, ket::InfinitePEPS, env::BPEnv) = - reduced_densitymatrix(inds, ket, ket, env) function reduced_densitymatrix1x1( ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 7ea89e06d..2ec2269fe 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -19,14 +19,14 @@ PEPS-operator sandwich is surrounded with the appropriate environment tensors. """ function contract_local_operator( inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ket::InfinitePEPS, bra::InfinitePEPS, env, ) static_inds = Tuple(Val.(inds)) return _contract_local_operator(static_inds, O, (ket, bra), env) end function contract_local_operator( inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ket::InfinitePEPS, bra::InfinitePEPS, env, ) return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env) end @@ -34,7 +34,7 @@ end Base.@deprecate( contract_local_operator( inds::NTuple, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, + ket::InfinitePEPS, bra::InfinitePEPS, env, ), contract_local_operator(collect(inds), O, ket, bra, env) ) @@ -269,13 +269,13 @@ on a rectangular patch based on `inds` but replacing the operator with an identi that the PEPS norm is computed. (Note that this is not the physical norm of the state.) """ function contract_local_norm( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env ) static_inds = Tuple(Val.(inds)) return _contract_local_norm(static_inds, (ket, bra), env) end function contract_local_norm( - inds::Vector{Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + inds::Vector{Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env ) return contract_local_norm(CartesianIndex.(inds), ket, bra, env) end @@ -305,7 +305,7 @@ end end Base.@deprecate( - contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv), + contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env), contract_local_norm(collect(inds), ket, bra, env) ) From 9b24a74ee2d978afc0c5401c04c4a6eb7e262438 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 2 Sep 2026 13:22:53 +0200 Subject: [PATCH 13/32] Reorganize `expectation_value` code - Reorganize file structure - Refactor in terms of generic expression generators that can be easily overloaded --- src/PEPSKit.jl | 1 + .../contractions/bp_contractions.jl | 11 +- src/algorithms/contractions/localoperator.jl | 655 +++++++++--------- src/algorithms/expectation_value.jl | 319 +++++++++ .../optimization/peps_optimization.jl | 15 + src/algorithms/toolbox.jl | 173 ----- test/toolbox/densitymatrices.jl | 26 +- 7 files changed, 689 insertions(+), 511 deletions(-) create mode 100644 src/algorithms/expectation_value.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b6e0e626b..8524e4946 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -140,6 +140,7 @@ include("algorithms/bp/beliefpropagation.jl") include("algorithms/bp/gaugefix.jl") include("algorithms/transfermatrix.jl") +include("algorithms/expectation_value.jl") include("algorithms/toolbox.jl") include("algorithms/correlator_adapters.jl") include("algorithms/correlators.jl") diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index a9ae6494f..f2bb7eb5d 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -54,8 +54,8 @@ absorb_west_message(A::PEPSTensor, M::PEPSMessage) = # rectangular-patch contraction used for a `CTMRGEnv`: only single sites and nearest neighbor # pairs can be contracted without introducing loop corrections. -function _contract_densitymatrix(inds::NTuple{N, Val}, state::Tuple, env::BPEnv) where {N} - sites = map(v -> typeof(v).parameters[1], inds) +function _contract_densitymatrix(inds::NTuple{N, Val}, state, env::BPEnv) where {N} + sites = _patch_inds(inds) return throw( ArgumentError( "Cannot contract a $(_patch_shape_string(sites)) patch using a `BPEnv`; @@ -64,13 +64,6 @@ function _contract_densitymatrix(inds::NTuple{N, Val}, state::Tuple, env::BPEnv) ) end -function _patch_shape_string(sites) - rows, cols = getindex.(sites, 1), getindex.(sites, 2) - nrows = maximum(rows) - minimum(rows) + 1 - ncols = maximum(cols) - minimum(cols) + 1 - return "$(nrows)x$(ncols)" -end - function reduced_densitymatrix1x1( ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 2ec2269fe..bd8a35fb3 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -1,5 +1,10 @@ # Contraction of local operators on arbitrary lattice locations # ------------------------------------------------------------- + + +# Contraction label helpers +# ------------------------- + function tensorlabel(args...) return Symbol(ntuple(i -> iseven(i) ? :_ : args[(i + 1) >> 1], 2 * length(args) - 1)...) end @@ -10,39 +15,141 @@ physicallabel(args...) = tensorlabel(:d, args...) """ $(SIGNATURES) -Contract a local operator `O` on the PEPS `ket`, `bra` at the indices `inds` using the -environment `env`. +Returns which slot of `open` the patch position `(r, c)` occupies, or `nothing` if that site +carries no open physical leg. `open` itself may be `nothing`, meaning no site does. +""" +_open_slot(open, r, c) = isnothing(open) ? nothing : findfirst(==(CartesianIndex(r, c)), open) -This works by generating the appropriate contraction on a rectangular patch with its corners -specified by `inds`. The `ket` and `bra` are contracted with `O` in between, and the -PEPS-operator sandwich is surrounded with the appropriate environment tensors. """ -function contract_local_operator( - inds::Vector{CartesianIndex{2}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env, +$(SIGNATURES) + +Returns the virtual leg labels of the bulk factor at patch position `(i, j)` of `layer`, in +domain order `(N, E, S, W)`. + +This is the bulk half of the label protocol, and is the same for every state type: legs facing +the perimeter of a `gridsize` patch take the `virtuallabel(SIDE, layer, ...)` labels that +[`boundary_contraction_expr`](@ref) consumes, while legs facing another site take an interior +`:horizontal` or `:vertical` label shared with that neighbour. +""" +function _bulk_virtuallabels(i, j, layer, gridsize) + return ( + i == 1 ? virtuallabel(NORTH, layer, j) : virtuallabel(:vertical, layer, i - 1, j), + j == gridsize[2] ? virtuallabel(EAST, layer, i) : + virtuallabel(:horizontal, layer, i, j), + i == gridsize[1] ? virtuallabel(SOUTH, layer, j) : virtuallabel(:vertical, layer, i, j), + j == 1 ? virtuallabel(WEST, layer, i) : virtuallabel(:horizontal, layer, i, j - 1), ) - static_inds = Tuple(Val.(inds)) - return _contract_local_operator(static_inds, O, (ket, bra), env) end -function contract_local_operator( - inds::Vector{Tuple{Int, Int}}, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env, - ) - return contract_local_operator(CartesianIndex.(inds), O, ket, bra, env) + + +# Patch geometry +# -------------- + +""" +$(SIGNATURES) + +Recover the patch coordinates from `Val`-encoded indices, checking that they do not overlap. +Uses an implementation in the type domain, so that the patch geometry is available to the +generated contraction expressions. +""" +function _patch_inds(inds::Type) + sites = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) + allunique(sites) || throw(ArgumentError("Indices should not overlap: $sites.")) + return sites end +_patch_inds(inds::Tuple{Vararg{Val}}) = _patch_inds(typeof(inds)) -Base.@deprecate( - contract_local_operator( - inds::NTuple, O::AbstractTensorMap, - ket::InfinitePEPS, bra::InfinitePEPS, env, - ), - contract_local_operator(collect(inds), O, ket, bra, env) -) +""" +$(SIGNATURES) + +Row and column ranges of the rectangular patch spanned by `sites`. +""" +function _patch_ranges(sites) + rows, cols = getindex.(sites, 1), getindex.(sites, 2) + return UnitRange(extrema(rows)...), UnitRange(extrema(cols)...) +end + +""" +$(SIGNATURES) + +Number of rows and columns of the rectangular patch spanned by `sites`. +""" +function _patch_gridsize(sites) + rowrange, colrange = _patch_ranges(sites) + return length(rowrange), length(colrange) +end + +""" +$(SIGNATURES) + +Shape of the patch spanned by `sites` as a string, e.g. `"2x2"`. For error messages and shape +guards of environments which only support a restricted patch geometry. +""" +function _patch_shape_string(sites) + nrows, ncols = _patch_gridsize(sites) + return "$(nrows)x$(ncols)" +end + + +# Assembly helpers +# ---------------- + +""" +$(SIGNATURES) + +Check whether the `@tensor` contraction for a given state type should be emitted with +`contractcheck = true`. Enabled for PEPO sandwiches, disabled for other state types, +in order to preserve the legacy behavior for different state types. +""" +contractcheck(state::Type) = false +contractcheck(::Type{<:InfinitePEPO}) = true +contractcheck(::Type{<:Tuple{Vararg{InfinitePEPO}}}) = true + +""" +$(SIGNATURES) + +Wrap an assembled product in `@autoopt @tensor`, with `lhs` on the left if given and as a +scalar otherwise, emitting `contractcheck = true` only when `check` is set. +""" +function _tensor_expr(prod, lhs, check::Bool) + return if isnothing(lhs) + check ? :(@autoopt @tensor contractcheck = true $prod) : :(@autoopt @tensor $prod) + else + check ? :(@autoopt @tensor contractcheck = true $lhs := $prod) : + :(@autoopt @tensor $lhs := $prod) + end +end + + +# Contraction expression generators +# --------------------------------- -# This implements the contraction of an operator acting on sites `inds`. -# The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance). +## Boundary: environment -function _contract_corner_expr(rowrange, colrange) +""" + boundary_contraction_expr(::Type{Env}, rowrange, colrange) + +Build the contraction expressions for the environment factors surrounding the patch spanned +by `rowrange` and `colrange`, dispatching on the type of the environment. + +The returned factors must consume exactly the perimeter labels the bulk exposes - +`virtuallabel(NORTH, layer, j)`, `virtuallabel(EAST, layer, i)`, +`virtuallabel(SOUTH, layer, j)` and `virtuallabel(WEST, layer, i)`, one per layer - and must +close every label they introduce themselves among their own factors. + +!!! note + The expression generator assumes the environment variable name is `env`. +""" +boundary_contraction_expr(env::Type, rowrange, colrange) = throw( + ArgumentError("No patch boundary contraction defined for environments of type $env.") +) + +function boundary_contraction_expr( + ::Type{<:CTMRGEnv{C, T}}, rowrange, colrange + ) where {C, T} + # the edges carry one virtual leg per layer of the sandwich, on top of the two environment + # indices threading the ring, so the height follows from the edge tensor type + height = numout(T) - 1 rmin, rmax = extrema(rowrange) cmin, cmax = extrema(colrange) gridsize = (rmax - rmin + 1, cmax - cmin + 1) @@ -59,14 +166,6 @@ function _contract_corner_expr(rowrange, colrange) C_SW = :(corner(env, SOUTHWEST, $(rmax + 1), $(cmin - 1))) corner_SW = tensorexpr(C_SW, envlabel(SOUTH, 0), envlabel(WEST, gridsize[1])) - return corner_NW, corner_NE, corner_SE, corner_SW -end - -function _contract_edge_expr(rowrange, colrange, height) - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - edges_N = map(1:gridsize[2]) do i E_N = :(edge(env, NORTH, $(rmin - 1), $(cmin + i - 1))) return tensorexpr( @@ -103,21 +202,46 @@ function _contract_edge_expr(rowrange, colrange, height) ) end - return edges_N, edges_E, edges_S, edges_W + return [ + corner_NW, corner_NE, corner_SE, corner_SW, + edges_N..., edges_E..., edges_S..., edges_W..., + ] end -function _contract_state_expr(rowrange, colrange, height, cartesian_inds = nothing) + +## Bulk: state + +""" + bulk_contraction_expr(::Type{State}, rowrange, colrange, open) + +Build the contraction expressions for the factors inside the patch spanned by `rowrange` +and `colrange`, dispatching on the type of the state. + +`open` is the vector of sites whose physical legs are left open, or `nothing` when every +physical leg is contracted within the bulk. Open sites are labelled +`physicallabel(:O, layer, slot)`, where `slot` is the position of the site in `open`; closed +sites share a label between the layers. Physical labels are this generator's business alone. + +Interior bonds must be closed among the returned factors, which must expose their +perimeter-facing virtual legs under the labels [`boundary_contraction_expr`](@ref) consumes. + +!!! note + The expression generator assumes the state variable name is `state`. +""" +bulk_contraction_expr(state::Type, rowrange, colrange, open) = throw( + ArgumentError("No patch bulk contraction defined for states of type $state.") +) + +function bulk_contraction_expr( + ::Type{<:Tuple{InfinitePEPS, InfinitePEPS}}, rowrange, colrange, open + ) rmin, rmax = extrema(rowrange) cmin, cmax = extrema(colrange) gridsize = (rmax - rmin + 1, cmax - cmin + 1) - return map(1:height) do side + layers = map(1:2) do side return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) - inds_id = if isnothing(cartesian_inds) - nothing - else - findfirst(==(CartesianIndex(rmin + i - 1, cmin + j - 1)), cartesian_inds) - end + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) physical_label = if isnothing(inds_id) physicallabel(i, j) else @@ -126,270 +250,221 @@ function _contract_state_expr(rowrange, colrange, height, cartesian_inds = nothi return tensorexpr( :(state[$(side)][$(rmin + i - 1), $(cmin + j - 1)]), (physical_label,), - ( - if i == 1 - virtuallabel(NORTH, side, j) - else - virtuallabel(:vertical, side, i - 1, j) - end, - if j == gridsize[2] - virtuallabel(EAST, side, i) - else - virtuallabel(:horizontal, side, i, j) - end, - if i == gridsize[1] - virtuallabel(SOUTH, side, j) - else - virtuallabel(:vertical, side, i, j) - end, - if j == 1 - virtuallabel(WEST, side, i) - else - virtuallabel(:horizontal, side, i, j - 1) - end, - ), + _bulk_virtuallabels(i, j, side, gridsize), ) end end + + ket, bra = layers + return [ket..., map(x -> Expr(:call, :conj, x), bra)...] end -function _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds = nothing) - @assert height == 1 || height == 2 "Contraction with more than 2 layers of PEPO is unintended." +function bulk_contraction_expr(::Type{<:InfinitePEPO}, rowrange, colrange, open) rmin, rmax = extrema(rowrange) cmin, cmax = extrema(colrange) gridsize = (rmax - rmin + 1, cmax - cmin + 1) - return map(1:height) do side + + # a single layer is not wrapped in a tuple, so it is indexed directly + layer = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) + physical_label_out = if isnothing(inds_id) + physicallabel(i, j) # traced over the layer + else + physicallabel(:O, 1, inds_id) + end + physical_label_in = if isnothing(inds_id) + physicallabel(i, j) + else + physicallabel(:O, 2, inds_id) + end + return tensorexpr( + :(twistdual(state[$(rmin + i - 1), $(cmin + j - 1)], 2)), + (physical_label_out, physical_label_in), + _bulk_virtuallabels(i, j, 1, gridsize), + ) + end + + return vec(layer) +end + +function bulk_contraction_expr( + ::Type{<:Tuple{InfinitePEPO, InfinitePEPO}}, rowrange, colrange, open + ) + rmin, rmax = extrema(rowrange) + cmin, cmax = extrema(colrange) + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + layers = map(1:2) do side return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) - inds_id = if isnothing(cartesian_inds) - nothing + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) + physical_label_out = if isnothing(inds_id) + physicallabel(:out, i, j) else - findfirst(==(CartesianIndex(rmin + i - 1, cmin + j - 1)), cartesian_inds) + physicallabel(:O, side, inds_id) end - if height == 1 - physical_label_in = if isnothing(inds_id) - physicallabel(i, j) - else - physicallabel(:O, 2, inds_id) - end - physical_label_out = if isnothing(inds_id) - physicallabel(i, j) - else - physicallabel(:O, 1, inds_id) - end - tensor_name = :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], 2)) + # the two layers are linked through a shared label on the open sites + physical_label_in = if isnothing(inds_id) + physicallabel(:in, i, j) else - physical_label_in = if isnothing(inds_id) - physicallabel(:in, i, j) - else - physicallabel(:Oopen, inds_id) - end - physical_label_out = if isnothing(inds_id) - physicallabel(:out, i, j) - else - physicallabel(:O, side, inds_id) - end - tensor_name = if side == 2 - :(state[2][$(rmin + i - 1), $(cmin + j - 1)]) - else - :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], (1, 2))) - end + physicallabel(:Oopen, inds_id) + end + tensor_name = if side == 2 + :(state[2][$(rmin + i - 1), $(cmin + j - 1)]) + else + :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], (1, 2))) end - return tensorexpr( tensor_name, (physical_label_out, physical_label_in), - ( - if i == 1 - virtuallabel(NORTH, side, j) - else - virtuallabel(:vertical, side, i - 1, j) - end, - if j == gridsize[2] - virtuallabel(EAST, side, i) - else - virtuallabel(:horizontal, side, i, j) - end, - if i == gridsize[1] - virtuallabel(SOUTH, side, j) - else - virtuallabel(:vertical, side, i, j) - end, - if j == 1 - virtuallabel(WEST, side, i) - else - virtuallabel(:horizontal, side, i, j - 1) - end, - ), + _bulk_virtuallabels(i, j, side, gridsize), ) end end + + ket, bra = layers + return [ket..., map(x -> Expr(:call, :conj, x), bra)...] end -@generated function _contract_local_operator( - inds::NTuple{N, Val}, - O::AbstractTensorMap{T, S, N, N}, - state::Tuple{InfinitePEPS, InfinitePEPS}, - env::CTMRGEnv, - ) where {T, S, N} - cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val - allunique(cartesian_inds) || - throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) - rowrange = getindex.(cartesian_inds, 1) - colrange = getindex.(cartesian_inds, 2) - - corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) - edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, 2) - operator = tensorexpr( - :O, - ntuple(i -> physicallabel(:O, 2, i), N), - ntuple(i -> physicallabel(:O, 1, i), N), +function bulk_contraction_expr( + state::Type{<:Tuple{Vararg{InfinitePEPO}}}, rowrange, colrange, open ) - ket, bra = _contract_state_expr(rowrange, colrange, 2, cartesian_inds) + return throw( + ArgumentError( + "Cannot contract a patch of $(length(state.parameters)) PEPO layers; only a \ + single layer or a two-layer sandwich are supported." + ) + ) +end + + +## Operator + +""" + operator_contraction_expr(::Type{Operator}, nsites) + +Build the contraction expressions for the operator factors inserted into a patch acting on +`nsites` sites, dispatching on the type of the operator. + +The factors refer to the enclosing generated function's `operator` argument by name, and +attach to the open physical legs of the bulk: `physicallabel(:O, 2, k)` is the bra index of +site `k` and `physicallabel(:O, 1, k)` its ket index. Every state type labels its open legs +with that same pair, so this generator does not depend on the state. Any label it introduces +itself must be closed among its own factors. + +!!! note + The expression generator assumes the operator variable name is `operator`. +""" +operator_contraction_expr(operator::Type, nsites) = throw( + ArgumentError("No patch operator contraction defined for operators of type $operator.") +) + +function operator_contraction_expr(::Type{<:AbstractTensorMap}, nsites) + return [ + tensorexpr( + :operator, + ntuple(i -> physicallabel(:O, 2, i), nsites), + ntuple(i -> physicallabel(:O, 1, i), nsites), + ), + ] +end + + +# Low-level patch contractions using generated contraction expressions +# -------------------------------------------------------------------- + +""" +$(SIGNATURES) + +Contract the rectangular patch spanned by `inds` with `operator` inserted on those sites, +leaving no physical leg open, and return the resulting scalar. + +The sites are carried as `Val` parameters so that the patch geometry is available while the +contraction is generated, and `state` is the tuple of layers making up the sandwich. The +expression is assembled from [`boundary_contraction_expr`](@ref), +[`bulk_contraction_expr`](@ref) and [`operator_contraction_expr`](@ref), which dispatch on the +types of `env`, `state` and `operator` respectively, so this single method covers every +combination those three have methods for. +""" +@generated function _contract_local_operator( + inds::NTuple{N, Val}, operator, state, env + ) where {N} + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) multiplication_ex = Expr( :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ket..., map(x -> Expr(:call, :conj, x), bra)..., - operator, + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, sites)..., + operator_contraction_expr(operator, N)..., ) - returnex = quote - @autoopt @tensor $multiplication_ex - end + returnex = _tensor_expr(multiplication_ex, nothing, contractcheck(state)) return macroexpand(@__MODULE__, returnex) end """ $(SIGNATURES) -Contract a local norm of the PEPS `ket`, `bra` around indices `inds`. +Contract the rectangular patch spanned by `inds` with no operator inserted, pairing the +physical legs of the layers on every site, and return the resulting scalar. -This works analogously to [`contract_local_operator`](@ref) by generating the contraction -on a rectangular patch based on `inds` but replacing the operator with an identity such -that the PEPS norm is computed. (Note that this is not the physical norm of the state.) +Assembled as [`_contract_local_operator`](@ref), except that +[`bulk_contraction_expr`](@ref) is passed `nothing` in place of the open sites, so no physical +leg is left open and no operator factor is needed. Note this is the norm of the patch within +the given environment, not the physical norm of the state. """ -function contract_local_norm( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) - static_inds = Tuple(Val.(inds)) - return _contract_local_norm(static_inds, (ket, bra), env) -end -function contract_local_norm( - inds::Vector{Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) - return contract_local_norm(CartesianIndex.(inds), ket, bra, env) -end @generated function _contract_local_norm( - inds::NTuple{N, Val}, state::Tuple{InfinitePEPS, InfinitePEPS}, env::CTMRGEnv + inds::NTuple{N, Val}, state, env ) where {N} - cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val - allunique(cartesian_inds) || throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) - rowrange = getindex.(cartesian_inds, 1) - colrange = getindex.(cartesian_inds, 2) - - corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) - edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, 2) - ket, bra = _contract_state_expr(rowrange, colrange, 2) + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) multiplication_ex = Expr( :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ket..., map(x -> Expr(:call, :conj, x), bra)..., + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, nothing)..., # legs paired, not open ) - returnex = quote - @autoopt @tensor $multiplication_ex - end + returnex = _tensor_expr(multiplication_ex, nothing, contractcheck(state)) return macroexpand(@__MODULE__, returnex) end -Base.@deprecate( - contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env), - contract_local_norm(collect(inds), ket, bra, env) -) - -@doc """ +""" $(SIGNATURES) -Construct the reduced density matrix `ρ` of PEPS or PEPO (representing PEPS with ancilla legs) `ket`, `bra` with open indices `inds` using the environment `env`. -Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the density matrix PEPO `state` with open indices `inds` using the environment `env`. +Contract the rectangular patch spanned by `inds` leaving the physical legs of those sites +open, and return the resulting reduced density matrix, normalized by its supertrace. -This works by generating the appropriate contraction on a rectangular patch with its corners -specified by `inds`. The result is normalized such that `tr(ρ) = 1`. -""" reduced_densitymatrix +Assembled as [`_contract_local_operator`](@ref) but without an operator factor, so the open +legs become the indices of `ρ`: `physicallabel(:O, 1, k)` in its codomain and +`physicallabel(:O, 2, k)` in its domain, for each site `k`. Normalization uses `str` rather +than `tr`, since the supertrace carries the fermionic signs. +""" +@generated function _contract_densitymatrix( + inds::NTuple{N, Val}, state, env + ) where {N} + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env + multiplication_ex = Expr( + :call, :*, + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, sites)..., + ) + result = tensorexpr( + :ρ, + ntuple(i -> physicallabel(:O, 1, i), N), + ntuple(i -> physicallabel(:O, 2, i), N), ) - length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) - if length(inds) == 2 - if inds[2] - inds[1] == CartesianIndex(1, 0) - return reduced_densitymatrix2x1(inds[1], ket, bra, env) - elseif inds[2] - inds[1] == CartesianIndex(0, 1) - return reduced_densitymatrix1x2(inds[1], ket, bra, env) - end + multex = _tensor_expr(multiplication_ex, result, contractcheck(state)) + return quote + $(macroexpand(@__MODULE__, multex)) + return ρ / str(ρ) end - - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, (ket, bra), env) -end -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env - ) - size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, (state,), env) end -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env - ) - size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) - size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, (ket, bra), env) -end -reduced_densitymatrix(inds, ket::InfinitePEPS, env) = - reduced_densitymatrix(inds, ket, ket, env) - -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) where {N}, - reduced_densitymatrix(collect(inds), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env - ) where {N}, - reduced_densitymatrix(collect(inds), state, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env - ) where {N}, - reduced_densitymatrix(collect(inds), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, ket::InfinitePEPO, bra::InfinitePEPO, env - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), state, env) -) +# Fast path specializations +# ------------------------- # Special case 1x1 density matrix: # Keep contraction order but try to optimize intermediate permutations: @@ -534,79 +609,3 @@ function reduced_densitymatrix1x2( return ρ / str(ρ) end - -@generated function _contract_densitymatrix( - inds::NTuple{N, Val}, state::Tuple{InfinitePEPS, InfinitePEPS}, env::CTMRGEnv - ) where {N} - cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val - allunique(cartesian_inds) || - throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) - rowrange = getindex.(cartesian_inds, 1) - colrange = getindex.(cartesian_inds, 2) - - corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) - edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, 2) - result = tensorexpr( - :ρ, - ntuple(i -> physicallabel(:O, 1, i), N), - ntuple(i -> physicallabel(:O, 2, i), N), - ) - ket, bra = _contract_state_expr(rowrange, colrange, 2, cartesian_inds) - - multiplication_ex = Expr( - :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ket..., map(x -> Expr(:call, :conj, x), bra)..., - ) - multex = :(@autoopt @tensor $result := $multiplication_ex) - return quote - $(macroexpand(@__MODULE__, multex)) - return ρ / str(ρ) - end -end - -@generated function _contract_densitymatrix( - inds::NTuple{N, Val}, state::Tuple{InfinitePEPO, Vararg{InfinitePEPO, M}}, env::CTMRGEnv - ) where {N, M} - height = M + 1 - @assert height == 1 || height == 2 "Contraction with more than 2 layers of PEPO is unintended." - cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val - allunique(cartesian_inds) || - throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) - rowrange = getindex.(cartesian_inds, 1) - colrange = getindex.(cartesian_inds, 2) - - corner_NW, corner_NE, corner_SE, corner_SW = _contract_corner_expr(rowrange, colrange) - edges_N, edges_E, edges_S, edges_W = _contract_edge_expr(rowrange, colrange, height) - result = tensorexpr( - :ρ, - ntuple(i -> physicallabel(:O, 1, i), N), - ntuple(i -> physicallabel(:O, 2, i), N), - ) - - pepos = _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds) - - if height == 2 - ket, bra = pepos - multiplication_ex = Expr( - :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ket..., map(x -> Expr(:call, :conj, x), bra)..., - ) - else - multiplication_ex = Expr( - :call, :*, - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - pepos[1]... - ) - end - - multex = :(@autoopt @tensor contractcheck = true $result := $multiplication_ex) - return quote - $(macroexpand(@__MODULE__, multex)) - return ρ / str(ρ) - end -end diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl new file mode 100644 index 000000000..5ed52da6d --- /dev/null +++ b/src/algorithms/expectation_value.jl @@ -0,0 +1,319 @@ +# Expectation value of a LocalOperator +# ------------------------------------ + +""" + expectation_value(state, O::LocalOperator, env) + expectation_value(bra, O::LocalOperator, ket, env) + +Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state) of a [`LocalOperator`](@ref) `O`. +This can be done either for a PEPS, or alternatively for a density matrix PEPO. +In the latter case the first signature corresponds to a single layer PEPO contraction, while +the second signature yields a bilayer contraction instead. +""" +function MPSKit.expectation_value( + bra::S, O::LocalOperator, ket::S, env + ) where {S <: InfiniteState} + checklattice(bra, O, ket) + term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly + return local_expectation_value(inds, bra, operator, ket, env) + end + return sum(term_vals) +end +MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) +function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) + checklattice(state, O) + term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly + return local_expectation_value(inds, state, operator, env) + end + return sum(term_vals) +end + + +# Expectation value of an individual local term +# --------------------------------------------- + +""" + local_expectation_value(inds, (bra, ket), operator, env) + local_expectation_value(inds, state, operator, env) + +Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation +value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state), where `operator` is the local +term acting on the sites `inds`. + +The implementation is overloaded based on the type of operator to be evaluated +""" +function local_expectation_value end + +# AbstractTensorMap evaluation goes through reduced density matrix +function local_expectation_value(inds, bra, operator::AbstractTensorMap, ket, env) + ρ = reduced_densitymatrix(inds, ket, bra, env) + return trmul(operator, ρ) +end +function local_expectation_value(inds, state, operator::AbstractTensorMap, env) + ρ = reduced_densitymatrix(inds, state, env) + return trmul(operator, ρ) +end + + +# Local patch contractions +# ------------------------ + +@doc """ + reduced_densitymatrix(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + reduced_densitymatrix(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) + reduced_densitymatrix(inds, state::InfinitePEPO, env) + +Construct the reduced density matrix `ρ` of `|ket⟩⟨bra|`, where both `ket` and `bra` +correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. +Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the +density matrix PEPO `state`. The reduced density matrix is contracted around the open +indices `inds` using the environment `env`, and is normalized such that `tr(ρ) = 1`. +""" reduced_densitymatrix + +# PEPS case has fast-path specializations +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) + + if length(inds) == 2 + if inds[2] - inds[1] == CartesianIndex(1, 0) + return reduced_densitymatrix2x1(inds[1], ket, bra, env) + elseif inds[2] - inds[1] == CartesianIndex(0, 1) + return reduced_densitymatrix1x2(inds[1], ket, bra, env) + end + end + + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, (ket, bra), env) +end +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env + ) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, state, env) +end +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, (ket, bra), env) +end +reduced_densitymatrix(inds, ket::InfinitePEPS, env) = + reduced_densitymatrix(inds, ket, ket, env) + +# handle deprecations of Tuple inds specifications +Base.@deprecate( + reduced_densitymatrix( + inds::NTuple{N, CartesianIndex{2}}, args... + ) where {N}, + reduced_densitymatrix(collect(inds), args...) +) +Base.@deprecate( + reduced_densitymatrix( + inds::NTuple{N, Tuple{Int, Int}}, args... + ) where {N}, + reduced_densitymatrix(collect(CartesianIndex.(inds)), args...) +) + +""" + contract_local_operator(inds, O, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_operator(inds, O, ket::InfinitePEPO, bra::InfinitePEPO, env) + contract_local_operator(inds, O, state::InfinitePEPO, env) + +Contract a local operator `O` between `ket` and `bra` states, computing `⟨bra|O|ket⟩`, where +`ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary +legs. Alternatively, contract a local operator `O` with a density matrix PEPO `state`, +computing `tr(O * state)`. `O` is applied to the open indices `inds`, and the result is +contracted using the environment `env`. +""" +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, + ket::InfinitePEPS, bra::InfinitePEPS, env, + ) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, state::InfinitePEPO, env + ) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, state, env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator(inds::Vector{Tuple{Int, Int}}, O, args...) + return contract_local_operator(CartesianIndex.(inds), O, args...) +end + +Base.@deprecate( + contract_local_operator( + inds::NTuple, args... + ), + contract_local_operator(collect(inds), args...) +) + +""" + contract_local_norm(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_norm(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) + contract_local_norm(inds, state::InfinitePEPO, env) + +Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch +of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a +PEPS with ancillary legs. Alternatively, contract a local norm patch of a density matrix +PEPO `state`, computing a patch of `tr(state)`. The norm patch is contracted around the open +indices `inds` using the environment `env`. +""" +function contract_local_norm( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env + ) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, (ket, bra), env) +end +function contract_local_norm(inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, state, env) +end +function contract_local_norm( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, (ket, bra), env) +end +function contract_local_norm(inds::Vector{Tuple{Int, Int}}, args...) + return contract_local_norm(CartesianIndex.(inds), args...) +end + +Base.@deprecate( + contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env), + contract_local_norm(collect(inds), ket, bra, env) +) + + +# Expectation value of a local partition function tensor +# ------------------------------------------------------ + +""" + expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) + +Compute the expectation value corresponding to inserting a local tensor(s) `O` at +position `inds` in the partition function `pf` and contracting the chole using a given CTMRG +environment `env`. + +Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O` +should be a rank-4 tensor conforming to the [`PartitionFunctionTensor`](@ref) indexing +convention. +""" +function MPSKit.expectation_value( + pf::InfinitePartitionFunction, + op::Pair{CartesianIndex{2}, <:AbstractTensorMap{T, S, 2, 2}}, + env, + ) where {T, S} + return contract_local_tensor(op[1], op[2], env) / + contract_local_tensor(op[1], pf[op[1]], env) +end +function MPSKit.expectation_value( + pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env + ) + return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) +end + +# Network values +# -------------- + +""" + network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) + +Return the value (per unit cell) of a given contractible network contracted using a given +CTMRG environment. +""" +function network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) + return prod(Iterators.product(axes(network)...)) do (r, c) + return _contract_site((r, c), network, env) * _contract_corners((r, c), env) / + _contract_vertical_edges((r, c), env) / _contract_horizontal_edges((r, c), env) + end +end +network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state), env) + +""" + _contract_site(ind::Tuple{Int,Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) + +Contract around a single site `ind` of a square network using a given CTMRG environment. +""" +function _contract_site(ind::Tuple{Int, Int}, network, env::CTMRGEnv) + r, c = ind + return _contract_site( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), + network[r, c], + ) +end + +""" + _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract all corners around the south-east at position `ind` of the CTMRG +environment `env`. +""" +function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_corners( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r, c), + corner(env, SOUTHWEST, r, c - 1), + ) +end + +""" + _contract_vertical_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract the vertical edges and corners around the east edge at position `ind` of the +CTMRG environment `env`. +""" +function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_vertical_edges( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r + 1, c), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, EAST, r, c), + edge(env, WEST, r, c - 1), + ) +end + +""" + _contract_horizontal_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract the horizontal edges and corners around the south edge at position `ind` of the +CTMRG environment `env`. +""" +function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_horizontal_edges( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r, c + 1), + corner(env, SOUTHWEST, r, c - 1), + edge(env, NORTH, r - 1, c), + edge(env, SOUTH, r, c), + ) +end diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl index 2151ca957..634292d3e 100644 --- a/src/algorithms/optimization/peps_optimization.jl +++ b/src/algorithms/optimization/peps_optimization.jl @@ -1,3 +1,18 @@ +""" + cost_function(peps::InfinitePEPS, env, O::LocalOperator) + +Real part of expectation value of `O`, used as the cost function in variational PEPS optimization. +Prints a warning if the expectation value yields a finite imaginary part (up to a tolerance). +""" +function cost_function(peps::InfinitePEPS, env, O::LocalOperator) + E = MPSKit.expectation_value(peps, O, env) + ignore_derivatives() do + isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || + @warn "Expectation value is not real: $E." + end + return real(E) +end + """ $(TYPEDEF) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index b44ca4f0f..2b93d85c0 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -1,180 +1,7 @@ -""" - expectation_value(state, O::LocalOperator, env) - expectation_value(bra, O::LocalOperator, ket, env) - -Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ of a [`LocalOperator`](@ref) `O`. -This can be done either for a PEPS, or alternatively for a density matrix PEPO. -In the latter case the first signature corresponds to a single layer PEPO contraction, while -the second signature yields a bilayer contraction instead. -""" -function MPSKit.expectation_value( - bra::S, O::LocalOperator, ket::S, env - ) where {S <: InfiniteState} - checklattice(bra, O, ket) - term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - return local_expectation_value(inds, bra, operator, ket, env) - end - return sum(term_vals) -end -MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) -function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) - checklattice(state, O) - term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - return local_expectation_value(inds, state, operator, env) - end - return sum(term_vals) -end - -""" - local_expectation_value(inds, bra, operator, ket, env) - local_expectation_value(inds, state, operator, env) - -Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation -value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state), where `operator` is the local -term acting on the sites `inds`. - -The implementation is overloaded based on the type of operator to be evaluated -""" -function local_expectation_value end - -# AbstractTensorMap evaluation goes through reduced density matrix -function local_expectation_value(inds, bra, operator::AbstractTensorMap, ket, env) - ρ = reduced_densitymatrix(inds, ket, bra, env) - return trmul(operator, ρ) -end -function local_expectation_value(inds, state, operator::AbstractTensorMap, env) - ρ = reduced_densitymatrix(inds, state, env) - return trmul(operator, ρ) -end - -""" - expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) - -Compute the expectation value corresponding to inserting a local tensor(s) `O` at -position `inds` in the partition function `pf` and contracting the chole using a given CTMRG -environment `env`. - -Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O` -should be a rank-4 tensor conforming to the [`PartitionFunctionTensor`](@ref) indexing -convention. -""" -function MPSKit.expectation_value( - pf::InfinitePartitionFunction, - op::Pair{CartesianIndex{2}, <:AbstractTensorMap{T, S, 2, 2}}, - env, - ) where {T, S} - return contract_local_tensor(op[1], op[2], env) / - contract_local_tensor(op[1], pf[op[1]], env) -end -function MPSKit.expectation_value( - pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env - ) - return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) -end - -""" - cost_function(peps::InfinitePEPS, env, O::LocalOperator) - -Real part of expectation value of `O`. Prints a warning if the expectation value -yields a finite imaginary part (up to a tolerance). -""" -function cost_function(peps::InfinitePEPS, env, O::LocalOperator) - E = MPSKit.expectation_value(peps, O, env) - ignore_derivatives() do - isapprox(imag(E), 0; atol = sqrt(eps(real(E)))) || - @warn "Expectation value is not real: $E." - end - return real(E) -end - function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) return network_value(InfiniteSquareNetwork(peps), env) end -""" - network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) - -Return the value (per unit cell) of a given contractible network contracted using a given -CTMRG environment. -""" -function network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) - return prod(Iterators.product(axes(network)...)) do (r, c) - return _contract_site((r, c), network, env) * _contract_corners((r, c), env) / - _contract_vertical_edges((r, c), env) / _contract_horizontal_edges((r, c), env) - end -end -network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state), env) - -""" - _contract_site(ind::Tuple{Int,Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) - -Contract around a single site `ind` of a square network using a given CTMRG environment. -""" -function _contract_site(ind::Tuple{Int, Int}, network, env::CTMRGEnv) - r, c = ind - return _contract_site( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r + 1, c + 1), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), - edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), - network[r, c], - ) -end - -""" - _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract all corners around the south-east at position `ind` of the CTMRG -environment `env`. -""" -function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_corners( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c), - corner(env, SOUTHEAST, r, c), - corner(env, SOUTHWEST, r, c - 1), - ) -end - -""" - _contract_vertical_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract the vertical edges and corners around the east edge at position `ind` of the -CTMRG environment `env`. -""" -function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_vertical_edges( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c), - corner(env, SOUTHEAST, r + 1, c), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, EAST, r, c), - edge(env, WEST, r, c - 1), - ) -end - -""" - _contract_horizontal_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract the horizontal edges and corners around the south edge at position `ind` of the -CTMRG environment `env`. -""" -function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_horizontal_edges( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r, c + 1), - corner(env, SOUTHWEST, r, c - 1), - edge(env, NORTH, r - 1, c), - edge(env, SOUTH, r, c), - ) -end - """ edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl index 803ab2283..8c3307a88 100644 --- a/test/toolbox/densitymatrices.jl +++ b/test/toolbox/densitymatrices.jl @@ -21,10 +21,16 @@ Ds = Dict(Trivial => ℂ^3, U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (1 @plansor O_pf[W S; N E] := O[p'; p] * ρ[1, 1, 1][p p'; N E S W] # Single site - O_singlesite = LocalOperator(physicalspace(ρ), ((1, 1),) => O) + site = (1, 1) + O_singlesite = LocalOperator(physicalspace(ρ), (site,) => O) E1 = expectation_value(ρ, O_singlesite, env) E2 = expectation_value(ρ_pf, CartesianIndex(1, 1) => O_pf, env) @test E1 ≈ E2 + # the operator/norm route agrees with the density matrix one; a single layer is passed + # to the contractions bare rather than wrapped in a tuple + val = contract_local_operator([site], O, ρ, env) + nrm = contract_local_norm([site], ρ, env) + @test E1 ≈ val / nrm # two sites for inds in zip( @@ -34,6 +40,9 @@ Ds = Dict(Trivial => ℂ^3, U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (1 O_twosite = LocalOperator(physicalspace(ρ), inds => O ⊗ O) E3 = expectation_value(ρ, O_twosite, env) # TODO: not defined for partition functions... + val = contract_local_operator(collect(inds), O ⊗ O, ρ, env) + nrm = contract_local_norm(collect(inds), ρ, env) + @test E3 ≈ val / nrm end end @@ -61,6 +70,10 @@ end nrm = contract_local_norm([site], ρ_peps, ρ_peps, env) @test nrm ≈ _contract_site(site, InfiniteSquareNetwork(ρ_peps), env) @test E1 ≈ val / nrm + # same through the two-layer PEPO sandwich rather than its fused PEPS view + val_pepo = contract_local_operator([site], O, ρ, ρ, env) + nrm_pepo = contract_local_norm([site], ρ, ρ, env) + @test E1 ≈ val_pepo / nrm_pepo # two sites for inds in zip( @@ -75,5 +88,16 @@ end val = contract_local_operator(collect(inds), O_doubled ⊗ O_doubled, ρ_peps, ρ_peps, env) nrm = contract_local_norm(collect(inds), ρ_peps, ρ_peps, env) @test E1 ≈ val / nrm + val_pepo = contract_local_operator(collect(inds), O ⊗ O, ρ, ρ, env) + nrm_pepo = contract_local_norm(collect(inds), ρ, ρ, env) + @test E1 ≈ val_pepo / nrm_pepo end end + +@testset "Three or more PEPO layers are rejected" begin + d, D, χ = ds[Trivial], Ds[Trivial], χs[Trivial] + ρ = InfinitePEPO(d, D; unitcell = (2, 2, 1)) + env = CTMRGEnv(InfinitePEPS(ρ), χ) + inds = Tuple(Val.([CartesianIndex(1, 1)])) + @test_throws ArgumentError PEPSKit._contract_densitymatrix(inds, (ρ, ρ, ρ), env) +end From 16bdf9dde1b575ec949ec14444bbc73712869508 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:21:16 +0200 Subject: [PATCH 14/32] Fix absorb docstrings --- src/utility/util.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/utility/util.jl b/src/utility/util.jl index 5b10c6ecc..8d2ba7ab8 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -64,10 +64,10 @@ end """ absorb_left( - E::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the left of a tensor `E` by contracting the first leg of `E` +Absorb a matrix `C` into the left of a tensor map `A` by contracting the first leg of `A` with the second leg of `C`. """ function absorb_left( @@ -89,7 +89,7 @@ end A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the right of a tensor `A` by contracting the last leg of `A` +Absorb a matrix `C` into the right of a tensor map `A` by contracting the last leg of `A` with the first leg of `C`. """ function absorb_right( @@ -108,18 +108,21 @@ end """ absorb_left_right( - A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + A::AbstractTensorMap{<:Any, S}, + CL::AbstractTensorMap{<:Any, S, 1, 1}, + CR::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the right of a tensor `A` by contracting the last leg of `A` -with the first leg of `C`. +Absorb matrices `CL` and `CR` into the left and right of a tensor map `A` by contracting the +first leg of `A` with the second leg of `CL` and the last leg of `A` with the first leg of +`CR`. """ function absorb_left_right( - T::AbstractTensorMap{<:Any, S}, + A::AbstractTensorMap{<:Any, S}, CL::AbstractTensorMap{<:Any, S, 1, 1}, CR::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} - return absorb_right(absorb_left(T, CL), CR) + return absorb_right(absorb_left(A, CL), CR) end _fliptwist_s(s::DiagonalTensorMap) = twist!(DiagonalTensorMap(flip(s, 1:2)), 1) From 3a0af9c524139ea812ee003931e026676314c1fa Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:22:00 +0200 Subject: [PATCH 15/32] Fix typo --- src/algorithms/expectation_value.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index 5ed52da6d..0c403c5f0 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -210,7 +210,7 @@ Base.@deprecate( expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) Compute the expectation value corresponding to inserting a local tensor(s) `O` at -position `inds` in the partition function `pf` and contracting the chole using a given CTMRG +position `inds` in the partition function `pf` and contracting the whole using a given CTMRG environment `env`. Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O` From f1b38d694accb212e18cfd7a3e1084dbefa76730 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:29:33 +0200 Subject: [PATCH 16/32] Fix up more docstrings --- src/algorithms/contractions/localoperator.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index bd8a35fb3..08a137262 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -374,12 +374,12 @@ end """ $(SIGNATURES) -Contract the rectangular patch spanned by `inds` with `operator` inserted on those sites, -leaving no physical leg open, and return the resulting scalar. +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with +`operator` inserted on those sites, leaving no physical leg open, and return the resulting +scalar. The sites are carried as `Val` parameters so that the patch geometry is available while the -contraction is generated, and `state` is the tuple of layers making up the sandwich. The -expression is assembled from [`boundary_contraction_expr`](@ref), +contraction is generated. The expression is assembled from [`boundary_contraction_expr`](@ref), [`bulk_contraction_expr`](@ref) and [`operator_contraction_expr`](@ref), which dispatch on the types of `env`, `state` and `operator` respectively, so this single method covers every combination those three have methods for. @@ -404,8 +404,9 @@ end """ $(SIGNATURES) -Contract the rectangular patch spanned by `inds` with no operator inserted, pairing the -physical legs of the layers on every site, and return the resulting scalar. +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with no +operator inserted, pairing the physical legs of the layers on every site, and return the +resulting scalar. Assembled as [`_contract_local_operator`](@ref), except that [`bulk_contraction_expr`](@ref) is passed `nothing` in place of the open sites, so no physical @@ -431,8 +432,8 @@ end """ $(SIGNATURES) -Contract the rectangular patch spanned by `inds` leaving the physical legs of those sites -open, and return the resulting reduced density matrix, normalized by its supertrace. +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, leaving +the physical legs of those sites open, and return the resulting reduced density matrix, normalized by its supertrace. Assembled as [`_contract_local_operator`](@ref) but without an operator factor, so the open legs become the indices of `ρ`: `physicallabel(:O, 1, k)` in its codomain and From 4d021cde097b990cd4b7bbcf4a19b78772bf526f Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:30:36 +0200 Subject: [PATCH 17/32] Add network type annotation --- src/algorithms/expectation_value.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index 0c403c5f0..08fb3e314 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -253,7 +253,7 @@ network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state) Contract around a single site `ind` of a square network using a given CTMRG environment. """ -function _contract_site(ind::Tuple{Int, Int}, network, env::CTMRGEnv) +function _contract_site(ind::Tuple{Int, Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) r, c = ind return _contract_site( corner(env, NORTHWEST, r - 1, c - 1), From c00fe0847ea3582f42ff9a37e5457defcf68bd64 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:40:33 +0200 Subject: [PATCH 18/32] Revert docstring to match implementation --- src/algorithms/expectation_value.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index 08fb3e314..1dc617407 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -33,7 +33,7 @@ end # --------------------------------------------- """ - local_expectation_value(inds, (bra, ket), operator, env) + local_expectation_value(inds, bra, operator, ket, env) local_expectation_value(inds, state, operator, env) Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation From d92d6c2bc5f05decdf2d2ef4ba5e9c729eb37693 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 11:41:00 +0200 Subject: [PATCH 19/32] Clarify use of supertrace and add explicit reference. --- src/algorithms/expectation_value.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index 1dc617407..ce381375e 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -67,7 +67,9 @@ Construct the reduced density matrix `ρ` of `|ket⟩⟨bra|`, where both `ket` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the density matrix PEPO `state`. The reduced density matrix is contracted around the open -indices `inds` using the environment `env`, and is normalized such that `tr(ρ) = 1`. +indices `inds` using the environment `env`, and is normalized such that `str(ρ) = 1`. + +See also [`str`](@ref). """ reduced_densitymatrix # PEPS case has fast-path specializations From 68b6514dd568df54f962c3982ed1df702079605c Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 4 Sep 2026 13:05:13 +0200 Subject: [PATCH 20/32] Try to make docstrings more specific --- src/algorithms/expectation_value.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index ce381375e..4e2e93ace 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -66,8 +66,9 @@ end Construct the reduced density matrix `ρ` of `|ket⟩⟨bra|`, where both `ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the -density matrix PEPO `state`. The reduced density matrix is contracted around the open -indices `inds` using the environment `env`, and is normalized such that `str(ρ) = 1`. +density matrix PEPO `state`. The reduced density matrix is contracted over the virtual +indices surrounding the open physical indices at sites `inds` using the environment `env`, +and is normalized such that `str(ρ) = 1`. See also [`str`](@ref). """ reduced_densitymatrix @@ -129,8 +130,8 @@ Base.@deprecate( Contract a local operator `O` between `ket` and `bra` states, computing `⟨bra|O|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. Alternatively, contract a local operator `O` with a density matrix PEPO `state`, -computing `tr(O * state)`. `O` is applied to the open indices `inds`, and the result is -contracted using the environment `env`. +computing `tr(O * state)`. `O` is applied to the open phyisical indices at sites `inds`, and +the result is contracted over the surrounding virtual indices using the environment `env`. """ function contract_local_operator( inds::Vector{CartesianIndex{2}}, O, @@ -173,8 +174,10 @@ Base.@deprecate( Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. Alternatively, contract a local norm patch of a density matrix -PEPO `state`, computing a patch of `tr(state)`. The norm patch is contracted around the open -indices `inds` using the environment `env`. +PEPO `state`, computing a patch of `tr(state)`. The contracted rectangular norm patch is +determined by the open physical indices `inds` and is contracted over the surrounding virtual +indices using the environment `env`. In particular, the patch location is precisely the same +as that of the patch used in [`contract_local_operator`](@ref). """ function contract_local_norm( inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env From cd81b82356ad96a118f90a65a1b8e8826fa143f9 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 11 Sep 2026 14:49:17 +0200 Subject: [PATCH 21/32] Update BP patch contraction comment --- src/algorithms/contractions/bp_contractions.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index f2bb7eb5d..8fedd3f05 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -50,9 +50,9 @@ absorb_west_message(A::PEPSTensor, M::PEPSMessage) = # Belief Propagation reduced density matrices # ------------------------------------------- -# BP messages live on the bonds of the network, so there is no analogue of the generated -# rectangular-patch contraction used for a `CTMRGEnv`: only single sites and nearest neighbor -# pairs can be contracted without introducing loop corrections. + +# NOTE: currently restricted to 1x1, 2x1, and 1x2 patches, since evaluating larger patches +# without including loop corrections tends to be a crude and not very useful appoximation. function _contract_densitymatrix(inds::NTuple{N, Val}, state, env::BPEnv) where {N} sites = _patch_inds(inds) From 1685a5704cb0fe0e11acfb162ec5002665d7c209 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 11 Sep 2026 14:49:49 +0200 Subject: [PATCH 22/32] Remove `contractcheck` altogether --- src/algorithms/contractions/localoperator.jl | 28 +++++--------------- 1 file changed, 6 insertions(+), 22 deletions(-) diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index 08a137262..76f521340 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -97,27 +97,11 @@ end """ $(SIGNATURES) -Check whether the `@tensor` contraction for a given state type should be emitted with -`contractcheck = true`. Enabled for PEPO sandwiches, disabled for other state types, -in order to preserve the legacy behavior for different state types. -""" -contractcheck(state::Type) = false -contractcheck(::Type{<:InfinitePEPO}) = true -contractcheck(::Type{<:Tuple{Vararg{InfinitePEPO}}}) = true - -""" -$(SIGNATURES) - Wrap an assembled product in `@autoopt @tensor`, with `lhs` on the left if given and as a -scalar otherwise, emitting `contractcheck = true` only when `check` is set. +scalar otherwise. """ -function _tensor_expr(prod, lhs, check::Bool) - return if isnothing(lhs) - check ? :(@autoopt @tensor contractcheck = true $prod) : :(@autoopt @tensor $prod) - else - check ? :(@autoopt @tensor contractcheck = true $lhs := $prod) : - :(@autoopt @tensor $lhs := $prod) - end +function _tensor_expr(prod, lhs = nothing) + return isnothing(lhs) ? :(@autoopt @tensor $prod) : :(@autoopt @tensor $lhs := $prod) end @@ -397,7 +381,7 @@ combination those three have methods for. operator_contraction_expr(operator, N)..., ) - returnex = _tensor_expr(multiplication_ex, nothing, contractcheck(state)) + returnex = _tensor_expr(multiplication_ex) return macroexpand(@__MODULE__, returnex) end @@ -425,7 +409,7 @@ the given environment, not the physical norm of the state. bulk_contraction_expr(state, rowrange, colrange, nothing)..., # legs paired, not open ) - returnex = _tensor_expr(multiplication_ex, nothing, contractcheck(state)) + returnex = _tensor_expr(multiplication_ex) return macroexpand(@__MODULE__, returnex) end @@ -457,7 +441,7 @@ than `tr`, since the supertrace carries the fermionic signs. ntuple(i -> physicallabel(:O, 2, i), N), ) - multex = _tensor_expr(multiplication_ex, result, contractcheck(state)) + multex = _tensor_expr(multiplication_ex, result) return quote $(macroexpand(@__MODULE__, multex)) return ρ / str(ρ) From d8ef1b6d3dd7e199f65ea85e61ea16f8ede42dbf Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 14 Sep 2026 10:21:56 +0200 Subject: [PATCH 23/32] Fix typo --- src/algorithms/expectation_value.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl index 4e2e93ace..a8f5048b8 100644 --- a/src/algorithms/expectation_value.jl +++ b/src/algorithms/expectation_value.jl @@ -130,7 +130,7 @@ Base.@deprecate( Contract a local operator `O` between `ket` and `bra` states, computing `⟨bra|O|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. Alternatively, contract a local operator `O` with a density matrix PEPO `state`, -computing `tr(O * state)`. `O` is applied to the open phyisical indices at sites `inds`, and +computing `tr(O * state)`. `O` is applied to the open physical indices at sites `inds`, and the result is contracted over the surrounding virtual indices using the environment `env`. """ function contract_local_operator( From 7f4b40b728f017c21cef2b3ba07a61634ce7ebc9 Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 14 Sep 2026 10:38:51 +0200 Subject: [PATCH 24/32] Reorganize files --- src/PEPSKit.jl | 20 +- src/algorithms/contractions/absorb.jl | 65 ++ src/algorithms/contractions/bp_messages.jl | 49 ++ .../contractions/ctmrg/network_value.jl | 70 ++ .../densitymatrix/bp.jl} | 50 -- .../local_patch/densitymatrix/ctmrg.jl | 145 +++++ .../local_patch/densitymatrix/generic.jl | 37 ++ .../contractions/local_patch/expr_utils.jl | 90 +++ .../contractions/local_patch/network_expr.jl | 262 ++++++++ .../local_patch/patch_contractions.jl | 60 ++ src/algorithms/contractions/localoperator.jl | 596 ------------------ src/algorithms/expectation_value.jl | 324 ---------- .../correlator_adapters.jl | 0 .../{ => expectation_value}/correlators.jl | 0 .../expectation_value/expectation_value.jl | 83 +++ .../expectation_value/network_value.jl | 73 +++ .../expectation_value/patch_contractions.jl | 87 +++ .../reduced_densitymatrix.jl | 109 ++++ src/algorithms/toolbox.jl | 56 -- src/utility/tensor_traces.jl | 27 + src/utility/util.jl | 88 --- 21 files changed, 1172 insertions(+), 1119 deletions(-) create mode 100644 src/algorithms/contractions/absorb.jl create mode 100644 src/algorithms/contractions/bp_messages.jl rename src/algorithms/contractions/{bp_contractions.jl => local_patch/densitymatrix/bp.jl} (60%) create mode 100644 src/algorithms/contractions/local_patch/densitymatrix/ctmrg.jl create mode 100644 src/algorithms/contractions/local_patch/densitymatrix/generic.jl create mode 100644 src/algorithms/contractions/local_patch/expr_utils.jl create mode 100644 src/algorithms/contractions/local_patch/network_expr.jl create mode 100644 src/algorithms/contractions/local_patch/patch_contractions.jl delete mode 100644 src/algorithms/contractions/localoperator.jl delete mode 100644 src/algorithms/expectation_value.jl rename src/algorithms/{ => expectation_value}/correlator_adapters.jl (100%) rename src/algorithms/{ => expectation_value}/correlators.jl (100%) create mode 100644 src/algorithms/expectation_value/expectation_value.jl create mode 100644 src/algorithms/expectation_value/network_value.jl create mode 100644 src/algorithms/expectation_value/patch_contractions.jl create mode 100644 src/algorithms/expectation_value/reduced_densitymatrix.jl create mode 100644 src/utility/tensor_traces.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index f4760367f..2cce6b142 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -52,6 +52,7 @@ using DocStringExtensions include("Defaults.jl") # Include first to allow for docstring interpolation with Defaults values include("utility/util.jl") +include("utility/tensor_traces.jl") include("utility/indexing.jl") include("utility/diffable_threads.jl") include("utility/twistdual.jl") @@ -99,11 +100,17 @@ include("algorithms/contractions/ctmrg/network_value.jl") include("algorithms/contractions/ctmrg/gaugefix.jl") include("algorithms/contractions/ctmrg/characteristic_equations.jl") +include("algorithms/contractions/absorb.jl") include("algorithms/contractions/absorb_weight.jl") include("algorithms/contractions/transfer.jl") -include("algorithms/contractions/localoperator.jl") include("algorithms/contractions/vumps_contractions.jl") -include("algorithms/contractions/bp_contractions.jl") +include("algorithms/contractions/bp_messages.jl") +include("algorithms/contractions/local_patch/expr_utils.jl") +include("algorithms/contractions/local_patch/network_expr.jl") +include("algorithms/contractions/local_patch/patch_contractions.jl") +include("algorithms/contractions/local_patch/densitymatrix/generic.jl") +include("algorithms/contractions/local_patch/densitymatrix/ctmrg.jl") +include("algorithms/contractions/local_patch/densitymatrix/bp.jl") include("algorithms/contractions/bondenv/benv_tools.jl") include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") @@ -141,10 +148,13 @@ include("algorithms/bp/beliefpropagation.jl") include("algorithms/bp/gaugefix.jl") include("algorithms/transfermatrix.jl") -include("algorithms/expectation_value.jl") +include("algorithms/expectation_value/expectation_value.jl") +include("algorithms/expectation_value/reduced_densitymatrix.jl") +include("algorithms/expectation_value/patch_contractions.jl") +include("algorithms/expectation_value/network_value.jl") +include("algorithms/expectation_value/correlator_adapters.jl") +include("algorithms/expectation_value/correlators.jl") include("algorithms/toolbox.jl") -include("algorithms/correlator_adapters.jl") -include("algorithms/correlators.jl") include("algorithms/optimization/implicit_differentiation.jl") include("algorithms/optimization/preconditioning.jl") diff --git a/src/algorithms/contractions/absorb.jl b/src/algorithms/contractions/absorb.jl new file mode 100644 index 000000000..b0387ef4c --- /dev/null +++ b/src/algorithms/contractions/absorb.jl @@ -0,0 +1,65 @@ +# Absorption of matrices into tensor legs +# --------------------------------------- + +""" + absorb_left( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb a matrix `C` into the left of a tensor map `A` by contracting the first leg of `A` +with the second leg of `C`. +""" +function absorb_left( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + pC = (codomainind(C), domainind(C)) + pA = ((codomainind(A)[1],), (codomainind(A)[2:end]..., domainind(A)...)) + pCA = (codomainind(A), domainind(A)) + return tensorcontract(C, pC, false, A, pA, false, pCA) +end +function absorb_left( + P::AbstractTensorMap{<:Any, S, 1, N}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S, N} + return twistnondual(C, 2) * P +end + +""" + absorb_right( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb a matrix `C` into the right of a tensor map `A` by contracting the last leg of `A` +with the first leg of `C`. +""" +function absorb_right( + A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + pA = ((codomainind(A)..., domainind(A)[2:end]...), (domainind(A)[1],)) + pC = (codomainind(C), domainind(C)) + pAC = (codomainind(A), (domainind(A)[end], domainind(A)[1:(end - 1)]...)) + return tensorcontract(A, pA, false, C, pC, false, pAC) +end +function absorb_right( + E::AbstractTensorMap{<:Any, S, N, 1}, C::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S, N} + return E * twistdual(C, 1) +end + +""" + absorb_left_right( + A::AbstractTensorMap{<:Any, S}, + CL::AbstractTensorMap{<:Any, S, 1, 1}, + CR::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + +Absorb matrices `CL` and `CR` into the left and right of a tensor map `A` by contracting the +first leg of `A` with the second leg of `CL` and the last leg of `A` with the first leg of +`CR`. +""" +function absorb_left_right( + A::AbstractTensorMap{<:Any, S}, + CL::AbstractTensorMap{<:Any, S, 1, 1}, + CR::AbstractTensorMap{<:Any, S, 1, 1} + ) where {S} + return absorb_right(absorb_left(A, CL), CR) +end diff --git a/src/algorithms/contractions/bp_messages.jl b/src/algorithms/contractions/bp_messages.jl new file mode 100644 index 000000000..55fdc51b9 --- /dev/null +++ b/src/algorithms/contractions/bp_messages.jl @@ -0,0 +1,49 @@ +const PEPSMessage = AbstractTensorMap{<:Any, <:Any, 1, 1} + +# Belief Propagation Updates +# -------------------------- +function contract_north_message( + A::PEPSSandwich, M_west::PEPSMessage, M_north::PEPSMessage, M_east::PEPSMessage + ) + return @autoopt @tensor begin + M_north′[DSt; DSb] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_west[DWb; DWt] * M_north[DNt; DNb] * M_east[DEt; DEb] + end +end +function contract_east_message( + A::PEPSSandwich, M_north::PEPSMessage, M_east::PEPSMessage, M_south::PEPSMessage + ) + return @autoopt @tensor begin + M_east′[DWt; DWb] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_north[DNt; DNb] * M_east[DEt; DEb] * M_south[DSb; DSt] + end +end +function contract_south_message( + A::PEPSSandwich, M_east::PEPSMessage, M_south::PEPSMessage, M_west::PEPSMessage + ) + return @autoopt @tensor begin + M_south′[DNb; DNt] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_east[DEt; DEb] * M_south[DSb; DSt] * M_west[DWb; DWt] + end +end +function contract_west_message( + A::PEPSSandwich, M_south::PEPSMessage, M_west::PEPSMessage, M_north::PEPSMessage + ) + return @autoopt @tensor begin + M_west′[DEb; DEt] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_south[DSb; DSt] * M_west[DWb; DWt] * M_north[DNt; DNb] + end +end + +absorb_north_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N' E S W] := A[d; N E S W] * M[N; N'] +absorb_east_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E' S W] := A[d; N E S W] * M[E; E'] +absorb_south_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E S' W] := A[d; N E S W] * M[S'; S] +absorb_west_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E S W'] := A[d; N E S W] * M[W'; W] diff --git a/src/algorithms/contractions/ctmrg/network_value.jl b/src/algorithms/contractions/ctmrg/network_value.jl index 9430a0df2..3961cc061 100644 --- a/src/algorithms/contractions/ctmrg/network_value.jl +++ b/src/algorithms/contractions/ctmrg/network_value.jl @@ -70,6 +70,24 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) end +""" + _contract_site(ind::Tuple{Int,Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) + +Contract around a single site `ind` of a square network using a given CTMRG environment. +""" +function _contract_site(ind::Tuple{Int, Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) + r, c = ind + return _contract_site( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), + network[r, c], + ) +end + ## Normalization contractions function _contract_corners( C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, @@ -79,6 +97,22 @@ function _contract_corners( C_southeast[3; 4] * C_southwest[4; 1] end +""" + _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract all corners around the south-east at position `ind` of the CTMRG +environment `env`. +""" +function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_corners( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r, c), + corner(env, SOUTHWEST, r, c - 1), + ) +end + @generated function _contract_vertical_edges( C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, @@ -105,6 +139,24 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) end +""" + _contract_vertical_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract the vertical edges and corners around the east edge at position `ind` of the +CTMRG environment `env`. +""" +function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_vertical_edges( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r + 1, c), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, EAST, r, c), + edge(env, WEST, r, c - 1), + ) +end + @generated function _contract_horizontal_edges( C_northwest::CTMRGCornerTensor, C_northeast::CTMRGCornerTensor, C_southeast::CTMRGCornerTensor, C_southwest::CTMRGCornerTensor, @@ -129,3 +181,21 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) end + +""" + _contract_horizontal_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) + +Contract the horizontal edges and corners around the south edge at position `ind` of the +CTMRG environment `env`. +""" +function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) + r, c = ind + return _contract_horizontal_edges( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r, c + 1), + corner(env, SOUTHWEST, r, c - 1), + edge(env, NORTH, r - 1, c), + edge(env, SOUTH, r, c), + ) +end diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/local_patch/densitymatrix/bp.jl similarity index 60% rename from src/algorithms/contractions/bp_contractions.jl rename to src/algorithms/contractions/local_patch/densitymatrix/bp.jl index 8fedd3f05..16c087b52 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/local_patch/densitymatrix/bp.jl @@ -1,53 +1,3 @@ -const PEPSMessage = AbstractTensorMap{<:Any, <:Any, 1, 1} - -# Belief Propagation Updates -# -------------------------- -function contract_north_message( - A::PEPSSandwich, M_west::PEPSMessage, M_north::PEPSMessage, M_east::PEPSMessage - ) - return @autoopt @tensor begin - M_north′[DSt; DSb] := - ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * - M_west[DWb; DWt] * M_north[DNt; DNb] * M_east[DEt; DEb] - end -end -function contract_east_message( - A::PEPSSandwich, M_north::PEPSMessage, M_east::PEPSMessage, M_south::PEPSMessage - ) - return @autoopt @tensor begin - M_east′[DWt; DWb] := - ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * - M_north[DNt; DNb] * M_east[DEt; DEb] * M_south[DSb; DSt] - end -end -function contract_south_message( - A::PEPSSandwich, M_east::PEPSMessage, M_south::PEPSMessage, M_west::PEPSMessage - ) - return @autoopt @tensor begin - M_south′[DNb; DNt] := - ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * - M_east[DEt; DEb] * M_south[DSb; DSt] * M_west[DWb; DWt] - end -end -function contract_west_message( - A::PEPSSandwich, M_south::PEPSMessage, M_west::PEPSMessage, M_north::PEPSMessage - ) - return @autoopt @tensor begin - M_west′[DEb; DEt] := - ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * - M_south[DSb; DSt] * M_west[DWb; DWt] * M_north[DNt; DNb] - end -end - -absorb_north_message(A::PEPSTensor, M::PEPSMessage) = - @tensor A′[d; N' E S W] := A[d; N E S W] * M[N; N'] -absorb_east_message(A::PEPSTensor, M::PEPSMessage) = - @tensor A′[d; N E' S W] := A[d; N E S W] * M[E; E'] -absorb_south_message(A::PEPSTensor, M::PEPSMessage) = - @tensor A′[d; N E S' W] := A[d; N E S W] * M[S'; S] -absorb_west_message(A::PEPSTensor, M::PEPSMessage) = - @tensor A′[d; N E S W'] := A[d; N E S W] * M[W'; W] - # Belief Propagation reduced density matrices # ------------------------------------------- diff --git a/src/algorithms/contractions/local_patch/densitymatrix/ctmrg.jl b/src/algorithms/contractions/local_patch/densitymatrix/ctmrg.jl new file mode 100644 index 000000000..823cb579f --- /dev/null +++ b/src/algorithms/contractions/local_patch/densitymatrix/ctmrg.jl @@ -0,0 +1,145 @@ +# Fast path reduced density matrices using CTMRG environments +# ----------------------------------------------------------- + +# Keep contraction order but try to optimize intermediate permutations: +# EE_SWA is largest object so keep largest legs to the front there +function reduced_densitymatrix1x1( + inds::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + ) + row, col = Tuple(inds) + + # Unpack variables and absorb corners + A = ket[row, col] + Ā = bra[row, col] + + E_north = absorb_right( + edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) + ) + E_east = absorb_right( + edge(env, EAST, row, col + 1), corner(env, SOUTHEAST, row + 1, col + 1) + ) + E_south = absorb_right( + edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) + ) + E_west = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) + + @tensor EE_SW[χSE χNW DSb DWb; DSt DWt] := + E_south[χSE DSt DSb; χSW] * E_west[χSW DWt DWb; χNW] + + @tensor EE_SWA[χSE χNW DNt DEt; dt DSb DWb] := + EE_SW[χSE χNW DSb DWb; DSt DWt] * A[dt; DNt DEt DSt DWt] + + @tensor EE_NE[DNb DEb; χSE χNW DNt DEt] := + E_north[χNW DNt DNb; χNE] * E_east[χNE DEt DEb; χSE] + + @tensor EEAEE[dt; DNb DEb DSb DWb] := + EE_NE[DNb DEb; χSE χNW DNt DEt] * EE_SWA[χSE χNW DNt DEt; dt DSb DWb] + + @tensor ρ[dt; db] := EEAEE[dt; DNb DEb DSb DWb] * conj(Ā[db; DNb DEb DSb DWb]) + + return ρ / str(ρ) +end + +# Special case 2x1 density matrix: +# Keep contraction order but try to optimize intermediate permutations: +function reduced_densitymatrix2x1( + ind::CartesianIndex, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + ) + row, col = Tuple(ind) + + # Unpack variables and absorb corners + A_north = ket[row, col] + Ā_north = bra[row, col] + A_south = ket[row + 1, col] + Ā_south = bra[row + 1, col] + + E_north = absorb_right( + edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) + ) + E_northeast = edge(env, EAST, row, col + 1) + E_southeast = absorb_right( + edge(env, EAST, row + 1, col + 1), corner(env, SOUTHEAST, row + 2, col + 1) + ) + E_south = absorb_right( + edge(env, SOUTH, row + 2, col), corner(env, SOUTHWEST, row + 2, col - 1) + ) + E_southwest = edge(env, WEST, row + 1, col - 1) + E_northwest = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) + + @tensor EE_NW[χW χNE DNWt DNt; DNWb DNb] := + E_northwest[χW DNWt DNWb; χNW] * E_north[χNW DNt DNb; χNE] + @tensor EEA_NW[χW DMb dNb χNE DNEb; DNWt DNt] := + EE_NW[χW χNE DNWt DNt; DNWb DNb] * conj(Ā_north[dNb; DNb DNEb DMb DNWb]) + @tensor EEAA_NW[χW DMb dNb dNt DMt; χNE DNEt DNEb] := + EEA_NW[χW DMb dNb χNE DNEb; DNWt DNt] * A_north[dNt; DNt DNEt DMt DNWt] + @tensor EEEAA_N[dNt dNb; χW DMt DMb χE] := + EEAA_NW[χW DMb dNb dNt DMt; χNE DNEt DNEb] * E_northeast[χNE DNEt DNEb; χE] + + @tensor EE_SE[χE χSW DSEt DSt; DSEb DSb] := + E_southeast[χE DSEt DSEb; χSE] * E_south[χSE DSt DSb; χSW] + @tensor EEA_SE[χE DMb dSb χSW DSWb; DSEt DSt] := + EE_SE[χE χSW DSEt DSt; DSEb DSb] * conj(Ā_south[dSb; DMb DSEb DSb DSWb]) + @tensor EEAA_SE[χE DMb dSb dSt DMt; χSW DSWt DSWb] := + EEA_SE[χE DMb dSb χSW DSWb; DSEt DSt] * A_south[dSt; DMt DSEt DSt DSWt] + @tensor EEEAA_S[χW DMt DMb χE; dSt dSb] := + EEAA_SE[χE DMb dSb dSt DMt; χSW DSWt DSWb] * E_southwest[χSW DSWt DSWb; χW] + + @tensor ρ[dNt dSt; dNb dSb] := + EEEAA_N[dNt dNb; χW DMt DMb χE] * EEEAA_S[χW DMt DMb χE; dSt dSb] + + return ρ / str(ρ) +end + +function reduced_densitymatrix1x2( + ind::CartesianIndex, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv + ) + row, col = Tuple(ind) + + # Unpack variables and absorb corners + A_west = ket[row, col] + Ā_west = bra[row, col] + A_east = ket[row, col + 1] + Ā_east = bra[row, col + 1] + + E_northwest = edge(env, NORTH, row - 1, col) + E_northeast = absorb_right( + edge(env, NORTH, row - 1, col + 1), corner(env, NORTHEAST, row - 1, col + 2) + ) + E_east = absorb_right( + edge(env, EAST, row, col + 2), corner(env, SOUTHEAST, row + 1, col + 2) + ) + E_southeast = edge(env, SOUTH, row + 1, col + 1) + E_southwest = absorb_right( + edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) + ) + E_west = absorb_right( + edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) + ) + + @tensor EE_SW[χS χNW DSWt DWt; DSWb DWb] := + E_southwest[χS DSWt DSWb; χSW] * E_west[χSW DWt DWb; χNW] + @tensor EEA_SW[χS DMb dWb χNW DNWb; DSWt DWt] := + EE_SW[χS χNW DSWt DWt; DSWb DWb] * conj(Ā_west[dWb; DNWb DMb DSWb DWb]) + @tensor EEAA_SW[χS DMb dWb dWt DMt; χNW DNWt DNWb] := + EEA_SW[χS DMb dWb χNW DNWb; DSWt DWt] * A_west[dWt; DNWt DMt DSWt DWt] + @tensor EEEAA_W[dWt dWb; χS DMt DMb χN] := + EEAA_SW[χS DMb dWb dWt DMt; χNW DNWt DNWb] * E_northwest[χNW DNWt DNWb; χN] + + @tensor EE_NE[χN χSE DNEt DEt; DNEb DEb] := + E_northeast[χN DNEt DNEb; χNE] * E_east[χNE DEt DEb; χSE] + @tensor EEA_NE[χN DMb dEb χSE DSEb; DNEt DEt] := + EE_NE[χN χSE DNEt DEt; DNEb DEb] * conj(Ā_east[dEb; DNEb DEb DSEb DMb]) + @tensor EEAA_NE[χN DMb dEb dEt DMt; χSE DSEt DSEb] := + EEA_NE[χN DMb dEb χSE DSEb; DNEt DEt] * A_east[dEt; DNEt DEt DSEt DMt] + @tensor EEEAA_E[χS DMt DMb χN; dEt dEb] := + EEAA_NE[χN DMb dEb dEt DMt; χSE DSEt DSEb] * E_southeast[χSE DSEt DSEb; χS] + + @tensor ρ[dWt dEt; dWb dEb] := + EEEAA_W[dWt dWb; χS DMt DMb χN] * EEEAA_E[χS DMt DMb χN; dEt dEb] + + return ρ / str(ρ) +end diff --git a/src/algorithms/contractions/local_patch/densitymatrix/generic.jl b/src/algorithms/contractions/local_patch/densitymatrix/generic.jl new file mode 100644 index 000000000..28f2c65ef --- /dev/null +++ b/src/algorithms/contractions/local_patch/densitymatrix/generic.jl @@ -0,0 +1,37 @@ +# Generic reduced density matrix patch contraction +# ------------------------------------------------ + +""" +$(SIGNATURES) + +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, leaving +the physical legs of those sites open, and return the resulting reduced density matrix, normalized by its supertrace. + +Assembled as [`_contract_local_operator`](@ref) but without an operator factor, so the open +legs become the indices of `ρ`: `physicallabel(:O, 1, k)` in its codomain and +`physicallabel(:O, 2, k)` in its domain, for each site `k`. Normalization uses `str` rather +than `tr`, since the supertrace carries the fermionic signs. +""" +@generated function _contract_densitymatrix( + inds::NTuple{N, Val}, state, env + ) where {N} + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) + + multiplication_ex = Expr( + :call, :*, + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, sites)..., + ) + result = tensorexpr( + :ρ, + ntuple(i -> physicallabel(:O, 1, i), N), + ntuple(i -> physicallabel(:O, 2, i), N), + ) + + multex = _tensor_expr(multiplication_ex, result) + return quote + $(macroexpand(@__MODULE__, multex)) + return ρ / str(ρ) + end +end diff --git a/src/algorithms/contractions/local_patch/expr_utils.jl b/src/algorithms/contractions/local_patch/expr_utils.jl new file mode 100644 index 000000000..6c0125614 --- /dev/null +++ b/src/algorithms/contractions/local_patch/expr_utils.jl @@ -0,0 +1,90 @@ +# Contraction expression utilities for local patches +# -------------------------------------------------- + +# Contraction label helpers +# ------------------------- + +function tensorlabel(args...) + return Symbol(ntuple(i -> iseven(i) ? :_ : args[(i + 1) >> 1], 2 * length(args) - 1)...) +end +envlabel(args...) = tensorlabel(:χ, args...) +virtuallabel(args...) = tensorlabel(:D, args...) +physicallabel(args...) = tensorlabel(:d, args...) + +""" +$(SIGNATURES) + +Returns which slot of `open` the patch position `(r, c)` occupies, or `nothing` if that site +carries no open physical leg. `open` itself may be `nothing`, meaning no site does. +""" +_open_slot(open, r, c) = isnothing(open) ? nothing : findfirst(==(CartesianIndex(r, c)), open) + +""" +$(SIGNATURES) + +Returns the virtual leg labels of the bulk factor at patch position `(i, j)` of `layer`, in +domain order `(N, E, S, W)`. + +This is the bulk half of the label protocol, and is the same for every state type: legs facing +the perimeter of a `gridsize` patch take the `virtuallabel(SIDE, layer, ...)` labels that +[`boundary_contraction_expr`](@ref) consumes, while legs facing another site take an interior +`:horizontal` or `:vertical` label shared with that neighbour. +""" +function _bulk_virtuallabels(i, j, layer, gridsize) + return ( + i == 1 ? virtuallabel(NORTH, layer, j) : virtuallabel(:vertical, layer, i - 1, j), + j == gridsize[2] ? virtuallabel(EAST, layer, i) : + virtuallabel(:horizontal, layer, i, j), + i == gridsize[1] ? virtuallabel(SOUTH, layer, j) : virtuallabel(:vertical, layer, i, j), + j == 1 ? virtuallabel(WEST, layer, i) : virtuallabel(:horizontal, layer, i, j - 1), + ) +end + + +# Patch geometry +# -------------- + +""" +$(SIGNATURES) + +Recover the patch coordinates from `Val`-encoded indices, checking that they do not overlap. +Uses an implementation in the type domain, so that the patch geometry is available to the +generated contraction expressions. +""" +function _patch_inds(inds::Type) + sites = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) + allunique(sites) || throw(ArgumentError("Indices should not overlap: $sites.")) + return sites +end +_patch_inds(inds::Tuple{Vararg{Val}}) = _patch_inds(typeof(inds)) + +""" +$(SIGNATURES) + +Row and column ranges of the rectangular patch spanned by `sites`. +""" +function _patch_ranges(sites) + rows, cols = getindex.(sites, 1), getindex.(sites, 2) + return UnitRange(extrema(rows)...), UnitRange(extrema(cols)...) +end + +""" +$(SIGNATURES) + +Number of rows and columns of the rectangular patch spanned by `sites`. +""" +function _patch_gridsize(sites) + rowrange, colrange = _patch_ranges(sites) + return length(rowrange), length(colrange) +end + +""" +$(SIGNATURES) + +Shape of the patch spanned by `sites` as a string, e.g. `"2x2"`. For error messages and shape +guards of environments which only support a restricted patch geometry. +""" +function _patch_shape_string(sites) + nrows, ncols = _patch_gridsize(sites) + return "$(nrows)x$(ncols)" +end diff --git a/src/algorithms/contractions/local_patch/network_expr.jl b/src/algorithms/contractions/local_patch/network_expr.jl new file mode 100644 index 000000000..6120fa3b3 --- /dev/null +++ b/src/algorithms/contractions/local_patch/network_expr.jl @@ -0,0 +1,262 @@ +# Local patch contraction expression generators +# --------------------------------------------- + +# Assembly helpers +# ---------------- + +""" +$(SIGNATURES) + +Wrap an assembled product in `@autoopt @tensor`, with `lhs` on the left if given and as a +scalar otherwise. +""" +function _tensor_expr(prod, lhs = nothing) + return isnothing(lhs) ? :(@autoopt @tensor $prod) : :(@autoopt @tensor $lhs := $prod) +end + + +# Contraction expression generators +# --------------------------------- + +## Boundary: environment + +""" + boundary_contraction_expr(::Type{Env}, rowrange, colrange) + +Build the contraction expressions for the environment factors surrounding the patch spanned +by `rowrange` and `colrange`, dispatching on the type of the environment. + +The returned factors must consume exactly the perimeter labels the bulk exposes - +`virtuallabel(NORTH, layer, j)`, `virtuallabel(EAST, layer, i)`, +`virtuallabel(SOUTH, layer, j)` and `virtuallabel(WEST, layer, i)`, one per layer - and must +close every label they introduce themselves among their own factors. + +!!! note + The expression generator assumes the environment variable name is `env`. +""" +boundary_contraction_expr(env::Type, rowrange, colrange) = throw( + ArgumentError("No patch boundary contraction defined for environments of type $env.") +) + +function boundary_contraction_expr( + ::Type{<:CTMRGEnv{C, T}}, rowrange, colrange + ) where {C, T} + # the edges carry one virtual leg per layer of the sandwich, on top of the two environment + # indices threading the ring, so the height follows from the edge tensor type + height = numout(T) - 1 + rmin, rmax = extrema(rowrange) + cmin, cmax = extrema(colrange) + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + C_NW = :(corner(env, NORTHWEST, $(rmin - 1), $(cmin - 1))) + corner_NW = tensorexpr(C_NW, envlabel(WEST, 0), envlabel(NORTH, 0)) + + C_NE = :(corner(env, NORTHEAST, $(rmin - 1), $(cmax + 1))) + corner_NE = tensorexpr(C_NE, envlabel(NORTH, gridsize[2]), envlabel(EAST, 0)) + + C_SE = :(corner(env, SOUTHEAST, $(rmax + 1), $(cmax + 1))) + corner_SE = tensorexpr(C_SE, envlabel(EAST, gridsize[1]), envlabel(SOUTH, gridsize[2])) + + C_SW = :(corner(env, SOUTHWEST, $(rmax + 1), $(cmin - 1))) + corner_SW = tensorexpr(C_SW, envlabel(SOUTH, 0), envlabel(WEST, gridsize[1])) + + edges_N = map(1:gridsize[2]) do i + E_N = :(edge(env, NORTH, $(rmin - 1), $(cmin + i - 1))) + return tensorexpr( + E_N, + (envlabel(NORTH, i - 1), virtuallabel.(NORTH, ntuple(identity, height), i)...), + envlabel(NORTH, i), + ) + end + + edges_E = map(1:gridsize[1]) do i + E_E = :(edge(env, EAST, $(rmin + i - 1), $(cmax + 1))) + return tensorexpr( + E_E, + (envlabel(EAST, i - 1), virtuallabel.(EAST, ntuple(identity, height), i)...), + envlabel(EAST, i), + ) + end + + edges_S = map(1:gridsize[2]) do i + E_S = :(edge(env, SOUTH, $(rmax + 1), $(cmin + i - 1))) + return tensorexpr( + E_S, + (envlabel(SOUTH, i), virtuallabel.(SOUTH, ntuple(identity, height), i)...), + envlabel(SOUTH, i - 1), + ) + end + + edges_W = map(1:gridsize[1]) do i + E_W = :(edge(env, WEST, $(rmin + i - 1), $(cmin - 1))) + return tensorexpr( + E_W, + (envlabel(WEST, i), virtuallabel.(WEST, ntuple(identity, height), i)...), + envlabel(WEST, i - 1), + ) + end + + return [ + corner_NW, corner_NE, corner_SE, corner_SW, + edges_N..., edges_E..., edges_S..., edges_W..., + ] +end + + +## Bulk: state + +""" + bulk_contraction_expr(::Type{State}, rowrange, colrange, open) + +Build the contraction expressions for the factors inside the patch spanned by `rowrange` +and `colrange`, dispatching on the type of the state. + +`open` is the vector of sites whose physical legs are left open, or `nothing` when every +physical leg is contracted within the bulk. Open sites are labelled +`physicallabel(:O, layer, slot)`, where `slot` is the position of the site in `open`; closed +sites share a label between the layers. Physical labels are this generator's business alone. + +Interior bonds must be closed among the returned factors, which must expose their +perimeter-facing virtual legs under the labels [`boundary_contraction_expr`](@ref) consumes. + +!!! note + The expression generator assumes the state variable name is `state`. +""" +bulk_contraction_expr(state::Type, rowrange, colrange, open) = throw( + ArgumentError("No patch bulk contraction defined for states of type $state.") +) + +function bulk_contraction_expr( + ::Type{<:Tuple{InfinitePEPS, InfinitePEPS}}, rowrange, colrange, open + ) + rmin, rmax = extrema(rowrange) + cmin, cmax = extrema(colrange) + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + layers = map(1:2) do side + return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) + physical_label = if isnothing(inds_id) + physicallabel(i, j) + else + physicallabel(:O, side, inds_id) + end + return tensorexpr( + :(state[$(side)][$(rmin + i - 1), $(cmin + j - 1)]), + (physical_label,), + _bulk_virtuallabels(i, j, side, gridsize), + ) + end + end + + ket, bra = layers + return [ket..., map(x -> Expr(:call, :conj, x), bra)...] +end + +function bulk_contraction_expr(::Type{<:InfinitePEPO}, rowrange, colrange, open) + rmin, rmax = extrema(rowrange) + cmin, cmax = extrema(colrange) + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + # a single layer is not wrapped in a tuple, so it is indexed directly + layer = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) + physical_label_out = if isnothing(inds_id) + physicallabel(i, j) # traced over the layer + else + physicallabel(:O, 1, inds_id) + end + physical_label_in = if isnothing(inds_id) + physicallabel(i, j) + else + physicallabel(:O, 2, inds_id) + end + return tensorexpr( + :(twistdual(state[$(rmin + i - 1), $(cmin + j - 1)], 2)), + (physical_label_out, physical_label_in), + _bulk_virtuallabels(i, j, 1, gridsize), + ) + end + + return vec(layer) +end + +function bulk_contraction_expr( + ::Type{<:Tuple{InfinitePEPO, InfinitePEPO}}, rowrange, colrange, open + ) + rmin, rmax = extrema(rowrange) + cmin, cmax = extrema(colrange) + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + layers = map(1:2) do side + return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) + physical_label_out = if isnothing(inds_id) + physicallabel(:out, i, j) + else + physicallabel(:O, side, inds_id) + end + # the two layers are linked through a shared label on the open sites + physical_label_in = if isnothing(inds_id) + physicallabel(:in, i, j) + else + physicallabel(:Oopen, inds_id) + end + tensor_name = if side == 2 + :(state[2][$(rmin + i - 1), $(cmin + j - 1)]) + else + :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], (1, 2))) + end + return tensorexpr( + tensor_name, (physical_label_out, physical_label_in), + _bulk_virtuallabels(i, j, side, gridsize), + ) + end + end + + ket, bra = layers + return [ket..., map(x -> Expr(:call, :conj, x), bra)...] +end + +function bulk_contraction_expr( + state::Type{<:Tuple{Vararg{InfinitePEPO}}}, rowrange, colrange, open + ) + return throw( + ArgumentError( + "Cannot contract a patch of $(length(state.parameters)) PEPO layers; only a \ + single layer or a two-layer sandwich are supported." + ) + ) +end + + +## Operator + +""" + operator_contraction_expr(::Type{Operator}, nsites) + +Build the contraction expressions for the operator factors inserted into a patch acting on +`nsites` sites, dispatching on the type of the operator. + +The factors refer to the enclosing generated function's `operator` argument by name, and +attach to the open physical legs of the bulk: `physicallabel(:O, 2, k)` is the bra index of +site `k` and `physicallabel(:O, 1, k)` its ket index. Every state type labels its open legs +with that same pair, so this generator does not depend on the state. Any label it introduces +itself must be closed among its own factors. + +!!! note + The expression generator assumes the operator variable name is `operator`. +""" +operator_contraction_expr(operator::Type, nsites) = throw( + ArgumentError("No patch operator contraction defined for operators of type $operator.") +) + +function operator_contraction_expr(::Type{<:AbstractTensorMap}, nsites) + return [ + tensorexpr( + :operator, + ntuple(i -> physicallabel(:O, 2, i), nsites), + ntuple(i -> physicallabel(:O, 1, i), nsites), + ), + ] +end diff --git a/src/algorithms/contractions/local_patch/patch_contractions.jl b/src/algorithms/contractions/local_patch/patch_contractions.jl new file mode 100644 index 000000000..caa830d73 --- /dev/null +++ b/src/algorithms/contractions/local_patch/patch_contractions.jl @@ -0,0 +1,60 @@ +# Low-level patch contractions using generated contraction expressions +# -------------------------------------------------------------------- + +""" +$(SIGNATURES) + +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with +`operator` inserted on those sites, leaving no physical leg open, and return the resulting +scalar. + +The sites are carried as `Val` parameters so that the patch geometry is available while the +contraction is generated. The expression is assembled from [`boundary_contraction_expr`](@ref), +[`bulk_contraction_expr`](@ref) and [`operator_contraction_expr`](@ref), which dispatch on the +types of `env`, `state` and `operator` respectively, so this single method covers every +combination those three have methods for. +""" +@generated function _contract_local_operator( + inds::NTuple{N, Val}, operator, state, env + ) where {N} + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) + + multiplication_ex = Expr( + :call, :*, + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, sites)..., + operator_contraction_expr(operator, N)..., + ) + + returnex = _tensor_expr(multiplication_ex) + return macroexpand(@__MODULE__, returnex) +end + +""" +$(SIGNATURES) + +Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with no +operator inserted, pairing the physical legs of the layers on every site, and return the +resulting scalar. + +Assembled as [`_contract_local_operator`](@ref), except that +[`bulk_contraction_expr`](@ref) is passed `nothing` in place of the open sites, so no physical +leg is left open and no operator factor is needed. Note this is the norm of the patch within +the given environment, not the physical norm of the state. +""" +@generated function _contract_local_norm( + inds::NTuple{N, Val}, state, env + ) where {N} + sites = _patch_inds(inds) + rowrange, colrange = _patch_ranges(sites) + + multiplication_ex = Expr( + :call, :*, + boundary_contraction_expr(env, rowrange, colrange)..., + bulk_contraction_expr(state, rowrange, colrange, nothing)..., # legs paired, not open + ) + + returnex = _tensor_expr(multiplication_ex) + return macroexpand(@__MODULE__, returnex) +end diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl deleted file mode 100644 index 76f521340..000000000 --- a/src/algorithms/contractions/localoperator.jl +++ /dev/null @@ -1,596 +0,0 @@ -# Contraction of local operators on arbitrary lattice locations -# ------------------------------------------------------------- - - -# Contraction label helpers -# ------------------------- - -function tensorlabel(args...) - return Symbol(ntuple(i -> iseven(i) ? :_ : args[(i + 1) >> 1], 2 * length(args) - 1)...) -end -envlabel(args...) = tensorlabel(:χ, args...) -virtuallabel(args...) = tensorlabel(:D, args...) -physicallabel(args...) = tensorlabel(:d, args...) - -""" -$(SIGNATURES) - -Returns which slot of `open` the patch position `(r, c)` occupies, or `nothing` if that site -carries no open physical leg. `open` itself may be `nothing`, meaning no site does. -""" -_open_slot(open, r, c) = isnothing(open) ? nothing : findfirst(==(CartesianIndex(r, c)), open) - -""" -$(SIGNATURES) - -Returns the virtual leg labels of the bulk factor at patch position `(i, j)` of `layer`, in -domain order `(N, E, S, W)`. - -This is the bulk half of the label protocol, and is the same for every state type: legs facing -the perimeter of a `gridsize` patch take the `virtuallabel(SIDE, layer, ...)` labels that -[`boundary_contraction_expr`](@ref) consumes, while legs facing another site take an interior -`:horizontal` or `:vertical` label shared with that neighbour. -""" -function _bulk_virtuallabels(i, j, layer, gridsize) - return ( - i == 1 ? virtuallabel(NORTH, layer, j) : virtuallabel(:vertical, layer, i - 1, j), - j == gridsize[2] ? virtuallabel(EAST, layer, i) : - virtuallabel(:horizontal, layer, i, j), - i == gridsize[1] ? virtuallabel(SOUTH, layer, j) : virtuallabel(:vertical, layer, i, j), - j == 1 ? virtuallabel(WEST, layer, i) : virtuallabel(:horizontal, layer, i, j - 1), - ) -end - - -# Patch geometry -# -------------- - -""" -$(SIGNATURES) - -Recover the patch coordinates from `Val`-encoded indices, checking that they do not overlap. -Uses an implementation in the type domain, so that the patch geometry is available to the -generated contraction expressions. -""" -function _patch_inds(inds::Type) - sites = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) - allunique(sites) || throw(ArgumentError("Indices should not overlap: $sites.")) - return sites -end -_patch_inds(inds::Tuple{Vararg{Val}}) = _patch_inds(typeof(inds)) - -""" -$(SIGNATURES) - -Row and column ranges of the rectangular patch spanned by `sites`. -""" -function _patch_ranges(sites) - rows, cols = getindex.(sites, 1), getindex.(sites, 2) - return UnitRange(extrema(rows)...), UnitRange(extrema(cols)...) -end - -""" -$(SIGNATURES) - -Number of rows and columns of the rectangular patch spanned by `sites`. -""" -function _patch_gridsize(sites) - rowrange, colrange = _patch_ranges(sites) - return length(rowrange), length(colrange) -end - -""" -$(SIGNATURES) - -Shape of the patch spanned by `sites` as a string, e.g. `"2x2"`. For error messages and shape -guards of environments which only support a restricted patch geometry. -""" -function _patch_shape_string(sites) - nrows, ncols = _patch_gridsize(sites) - return "$(nrows)x$(ncols)" -end - - -# Assembly helpers -# ---------------- - -""" -$(SIGNATURES) - -Wrap an assembled product in `@autoopt @tensor`, with `lhs` on the left if given and as a -scalar otherwise. -""" -function _tensor_expr(prod, lhs = nothing) - return isnothing(lhs) ? :(@autoopt @tensor $prod) : :(@autoopt @tensor $lhs := $prod) -end - - -# Contraction expression generators -# --------------------------------- - -## Boundary: environment - -""" - boundary_contraction_expr(::Type{Env}, rowrange, colrange) - -Build the contraction expressions for the environment factors surrounding the patch spanned -by `rowrange` and `colrange`, dispatching on the type of the environment. - -The returned factors must consume exactly the perimeter labels the bulk exposes - -`virtuallabel(NORTH, layer, j)`, `virtuallabel(EAST, layer, i)`, -`virtuallabel(SOUTH, layer, j)` and `virtuallabel(WEST, layer, i)`, one per layer - and must -close every label they introduce themselves among their own factors. - -!!! note - The expression generator assumes the environment variable name is `env`. -""" -boundary_contraction_expr(env::Type, rowrange, colrange) = throw( - ArgumentError("No patch boundary contraction defined for environments of type $env.") -) - -function boundary_contraction_expr( - ::Type{<:CTMRGEnv{C, T}}, rowrange, colrange - ) where {C, T} - # the edges carry one virtual leg per layer of the sandwich, on top of the two environment - # indices threading the ring, so the height follows from the edge tensor type - height = numout(T) - 1 - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - - C_NW = :(corner(env, NORTHWEST, $(rmin - 1), $(cmin - 1))) - corner_NW = tensorexpr(C_NW, envlabel(WEST, 0), envlabel(NORTH, 0)) - - C_NE = :(corner(env, NORTHEAST, $(rmin - 1), $(cmax + 1))) - corner_NE = tensorexpr(C_NE, envlabel(NORTH, gridsize[2]), envlabel(EAST, 0)) - - C_SE = :(corner(env, SOUTHEAST, $(rmax + 1), $(cmax + 1))) - corner_SE = tensorexpr(C_SE, envlabel(EAST, gridsize[1]), envlabel(SOUTH, gridsize[2])) - - C_SW = :(corner(env, SOUTHWEST, $(rmax + 1), $(cmin - 1))) - corner_SW = tensorexpr(C_SW, envlabel(SOUTH, 0), envlabel(WEST, gridsize[1])) - - edges_N = map(1:gridsize[2]) do i - E_N = :(edge(env, NORTH, $(rmin - 1), $(cmin + i - 1))) - return tensorexpr( - E_N, - (envlabel(NORTH, i - 1), virtuallabel.(NORTH, ntuple(identity, height), i)...), - envlabel(NORTH, i), - ) - end - - edges_E = map(1:gridsize[1]) do i - E_E = :(edge(env, EAST, $(rmin + i - 1), $(cmax + 1))) - return tensorexpr( - E_E, - (envlabel(EAST, i - 1), virtuallabel.(EAST, ntuple(identity, height), i)...), - envlabel(EAST, i), - ) - end - - edges_S = map(1:gridsize[2]) do i - E_S = :(edge(env, SOUTH, $(rmax + 1), $(cmin + i - 1))) - return tensorexpr( - E_S, - (envlabel(SOUTH, i), virtuallabel.(SOUTH, ntuple(identity, height), i)...), - envlabel(SOUTH, i - 1), - ) - end - - edges_W = map(1:gridsize[1]) do i - E_W = :(edge(env, WEST, $(rmin + i - 1), $(cmin - 1))) - return tensorexpr( - E_W, - (envlabel(WEST, i), virtuallabel.(WEST, ntuple(identity, height), i)...), - envlabel(WEST, i - 1), - ) - end - - return [ - corner_NW, corner_NE, corner_SE, corner_SW, - edges_N..., edges_E..., edges_S..., edges_W..., - ] -end - - -## Bulk: state - -""" - bulk_contraction_expr(::Type{State}, rowrange, colrange, open) - -Build the contraction expressions for the factors inside the patch spanned by `rowrange` -and `colrange`, dispatching on the type of the state. - -`open` is the vector of sites whose physical legs are left open, or `nothing` when every -physical leg is contracted within the bulk. Open sites are labelled -`physicallabel(:O, layer, slot)`, where `slot` is the position of the site in `open`; closed -sites share a label between the layers. Physical labels are this generator's business alone. - -Interior bonds must be closed among the returned factors, which must expose their -perimeter-facing virtual legs under the labels [`boundary_contraction_expr`](@ref) consumes. - -!!! note - The expression generator assumes the state variable name is `state`. -""" -bulk_contraction_expr(state::Type, rowrange, colrange, open) = throw( - ArgumentError("No patch bulk contraction defined for states of type $state.") -) - -function bulk_contraction_expr( - ::Type{<:Tuple{InfinitePEPS, InfinitePEPS}}, rowrange, colrange, open - ) - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - - layers = map(1:2) do side - return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) - inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) - physical_label = if isnothing(inds_id) - physicallabel(i, j) - else - physicallabel(:O, side, inds_id) - end - return tensorexpr( - :(state[$(side)][$(rmin + i - 1), $(cmin + j - 1)]), - (physical_label,), - _bulk_virtuallabels(i, j, side, gridsize), - ) - end - end - - ket, bra = layers - return [ket..., map(x -> Expr(:call, :conj, x), bra)...] -end - -function bulk_contraction_expr(::Type{<:InfinitePEPO}, rowrange, colrange, open) - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - - # a single layer is not wrapped in a tuple, so it is indexed directly - layer = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) - inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) - physical_label_out = if isnothing(inds_id) - physicallabel(i, j) # traced over the layer - else - physicallabel(:O, 1, inds_id) - end - physical_label_in = if isnothing(inds_id) - physicallabel(i, j) - else - physicallabel(:O, 2, inds_id) - end - return tensorexpr( - :(twistdual(state[$(rmin + i - 1), $(cmin + j - 1)], 2)), - (physical_label_out, physical_label_in), - _bulk_virtuallabels(i, j, 1, gridsize), - ) - end - - return vec(layer) -end - -function bulk_contraction_expr( - ::Type{<:Tuple{InfinitePEPO, InfinitePEPO}}, rowrange, colrange, open - ) - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - - layers = map(1:2) do side - return map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) - inds_id = _open_slot(open, rmin + i - 1, cmin + j - 1) - physical_label_out = if isnothing(inds_id) - physicallabel(:out, i, j) - else - physicallabel(:O, side, inds_id) - end - # the two layers are linked through a shared label on the open sites - physical_label_in = if isnothing(inds_id) - physicallabel(:in, i, j) - else - physicallabel(:Oopen, inds_id) - end - tensor_name = if side == 2 - :(state[2][$(rmin + i - 1), $(cmin + j - 1)]) - else - :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], (1, 2))) - end - return tensorexpr( - tensor_name, (physical_label_out, physical_label_in), - _bulk_virtuallabels(i, j, side, gridsize), - ) - end - end - - ket, bra = layers - return [ket..., map(x -> Expr(:call, :conj, x), bra)...] -end - -function bulk_contraction_expr( - state::Type{<:Tuple{Vararg{InfinitePEPO}}}, rowrange, colrange, open - ) - return throw( - ArgumentError( - "Cannot contract a patch of $(length(state.parameters)) PEPO layers; only a \ - single layer or a two-layer sandwich are supported." - ) - ) -end - - -## Operator - -""" - operator_contraction_expr(::Type{Operator}, nsites) - -Build the contraction expressions for the operator factors inserted into a patch acting on -`nsites` sites, dispatching on the type of the operator. - -The factors refer to the enclosing generated function's `operator` argument by name, and -attach to the open physical legs of the bulk: `physicallabel(:O, 2, k)` is the bra index of -site `k` and `physicallabel(:O, 1, k)` its ket index. Every state type labels its open legs -with that same pair, so this generator does not depend on the state. Any label it introduces -itself must be closed among its own factors. - -!!! note - The expression generator assumes the operator variable name is `operator`. -""" -operator_contraction_expr(operator::Type, nsites) = throw( - ArgumentError("No patch operator contraction defined for operators of type $operator.") -) - -function operator_contraction_expr(::Type{<:AbstractTensorMap}, nsites) - return [ - tensorexpr( - :operator, - ntuple(i -> physicallabel(:O, 2, i), nsites), - ntuple(i -> physicallabel(:O, 1, i), nsites), - ), - ] -end - - -# Low-level patch contractions using generated contraction expressions -# -------------------------------------------------------------------- - -""" -$(SIGNATURES) - -Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with -`operator` inserted on those sites, leaving no physical leg open, and return the resulting -scalar. - -The sites are carried as `Val` parameters so that the patch geometry is available while the -contraction is generated. The expression is assembled from [`boundary_contraction_expr`](@ref), -[`bulk_contraction_expr`](@ref) and [`operator_contraction_expr`](@ref), which dispatch on the -types of `env`, `state` and `operator` respectively, so this single method covers every -combination those three have methods for. -""" -@generated function _contract_local_operator( - inds::NTuple{N, Val}, operator, state, env - ) where {N} - sites = _patch_inds(inds) - rowrange, colrange = _patch_ranges(sites) - - multiplication_ex = Expr( - :call, :*, - boundary_contraction_expr(env, rowrange, colrange)..., - bulk_contraction_expr(state, rowrange, colrange, sites)..., - operator_contraction_expr(operator, N)..., - ) - - returnex = _tensor_expr(multiplication_ex) - return macroexpand(@__MODULE__, returnex) -end - -""" -$(SIGNATURES) - -Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, with no -operator inserted, pairing the physical legs of the layers on every site, and return the -resulting scalar. - -Assembled as [`_contract_local_operator`](@ref), except that -[`bulk_contraction_expr`](@ref) is passed `nothing` in place of the open sites, so no physical -leg is left open and no operator factor is needed. Note this is the norm of the patch within -the given environment, not the physical norm of the state. -""" -@generated function _contract_local_norm( - inds::NTuple{N, Val}, state, env - ) where {N} - sites = _patch_inds(inds) - rowrange, colrange = _patch_ranges(sites) - - multiplication_ex = Expr( - :call, :*, - boundary_contraction_expr(env, rowrange, colrange)..., - bulk_contraction_expr(state, rowrange, colrange, nothing)..., # legs paired, not open - ) - - returnex = _tensor_expr(multiplication_ex) - return macroexpand(@__MODULE__, returnex) -end - -""" -$(SIGNATURES) - -Contract the rectangular patch of a network encoded by `state`, spanned by `inds`, leaving -the physical legs of those sites open, and return the resulting reduced density matrix, normalized by its supertrace. - -Assembled as [`_contract_local_operator`](@ref) but without an operator factor, so the open -legs become the indices of `ρ`: `physicallabel(:O, 1, k)` in its codomain and -`physicallabel(:O, 2, k)` in its domain, for each site `k`. Normalization uses `str` rather -than `tr`, since the supertrace carries the fermionic signs. -""" -@generated function _contract_densitymatrix( - inds::NTuple{N, Val}, state, env - ) where {N} - sites = _patch_inds(inds) - rowrange, colrange = _patch_ranges(sites) - - multiplication_ex = Expr( - :call, :*, - boundary_contraction_expr(env, rowrange, colrange)..., - bulk_contraction_expr(state, rowrange, colrange, sites)..., - ) - result = tensorexpr( - :ρ, - ntuple(i -> physicallabel(:O, 1, i), N), - ntuple(i -> physicallabel(:O, 2, i), N), - ) - - multex = _tensor_expr(multiplication_ex, result) - return quote - $(macroexpand(@__MODULE__, multex)) - return ρ / str(ρ) - end -end - -# Fast path specializations -# ------------------------- - -# Special case 1x1 density matrix: -# Keep contraction order but try to optimize intermediate permutations: -# EE_SWA is largest object so keep largest legs to the front there -function reduced_densitymatrix1x1( - inds::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) - row, col = Tuple(inds) - - # Unpack variables and absorb corners - A = ket[row, col] - Ā = bra[row, col] - - E_north = absorb_right( - edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) - ) - E_east = absorb_right( - edge(env, EAST, row, col + 1), corner(env, SOUTHEAST, row + 1, col + 1) - ) - E_south = absorb_right( - edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) - ) - E_west = absorb_right( - edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) - ) - - @tensor EE_SW[χSE χNW DSb DWb; DSt DWt] := - E_south[χSE DSt DSb; χSW] * E_west[χSW DWt DWb; χNW] - - @tensor EE_SWA[χSE χNW DNt DEt; dt DSb DWb] := - EE_SW[χSE χNW DSb DWb; DSt DWt] * A[dt; DNt DEt DSt DWt] - - @tensor EE_NE[DNb DEb; χSE χNW DNt DEt] := - E_north[χNW DNt DNb; χNE] * E_east[χNE DEt DEb; χSE] - - @tensor EEAEE[dt; DNb DEb DSb DWb] := - EE_NE[DNb DEb; χSE χNW DNt DEt] * EE_SWA[χSE χNW DNt DEt; dt DSb DWb] - - @tensor ρ[dt; db] := EEAEE[dt; DNb DEb DSb DWb] * conj(Ā[db; DNb DEb DSb DWb]) - - return ρ / str(ρ) -end - -# Special case 2x1 density matrix: -# Keep contraction order but try to optimize intermediate permutations: -function reduced_densitymatrix2x1( - ind::CartesianIndex, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) - row, col = Tuple(ind) - - # Unpack variables and absorb corners - A_north = ket[row, col] - Ā_north = bra[row, col] - A_south = ket[row + 1, col] - Ā_south = bra[row + 1, col] - - E_north = absorb_right( - edge(env, NORTH, row - 1, col), corner(env, NORTHEAST, row - 1, col + 1) - ) - E_northeast = edge(env, EAST, row, col + 1) - E_southeast = absorb_right( - edge(env, EAST, row + 1, col + 1), corner(env, SOUTHEAST, row + 2, col + 1) - ) - E_south = absorb_right( - edge(env, SOUTH, row + 2, col), corner(env, SOUTHWEST, row + 2, col - 1) - ) - E_southwest = edge(env, WEST, row + 1, col - 1) - E_northwest = absorb_right( - edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) - ) - - @tensor EE_NW[χW χNE DNWt DNt; DNWb DNb] := - E_northwest[χW DNWt DNWb; χNW] * E_north[χNW DNt DNb; χNE] - @tensor EEA_NW[χW DMb dNb χNE DNEb; DNWt DNt] := - EE_NW[χW χNE DNWt DNt; DNWb DNb] * conj(Ā_north[dNb; DNb DNEb DMb DNWb]) - @tensor EEAA_NW[χW DMb dNb dNt DMt; χNE DNEt DNEb] := - EEA_NW[χW DMb dNb χNE DNEb; DNWt DNt] * A_north[dNt; DNt DNEt DMt DNWt] - @tensor EEEAA_N[dNt dNb; χW DMt DMb χE] := - EEAA_NW[χW DMb dNb dNt DMt; χNE DNEt DNEb] * E_northeast[χNE DNEt DNEb; χE] - - @tensor EE_SE[χE χSW DSEt DSt; DSEb DSb] := - E_southeast[χE DSEt DSEb; χSE] * E_south[χSE DSt DSb; χSW] - @tensor EEA_SE[χE DMb dSb χSW DSWb; DSEt DSt] := - EE_SE[χE χSW DSEt DSt; DSEb DSb] * conj(Ā_south[dSb; DMb DSEb DSb DSWb]) - @tensor EEAA_SE[χE DMb dSb dSt DMt; χSW DSWt DSWb] := - EEA_SE[χE DMb dSb χSW DSWb; DSEt DSt] * A_south[dSt; DMt DSEt DSt DSWt] - @tensor EEEAA_S[χW DMt DMb χE; dSt dSb] := - EEAA_SE[χE DMb dSb dSt DMt; χSW DSWt DSWb] * E_southwest[χSW DSWt DSWb; χW] - - @tensor ρ[dNt dSt; dNb dSb] := - EEEAA_N[dNt dNb; χW DMt DMb χE] * EEEAA_S[χW DMt DMb χE; dSt dSb] - - return ρ / str(ρ) -end - -function reduced_densitymatrix1x2( - ind::CartesianIndex, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) - row, col = Tuple(ind) - - # Unpack variables and absorb corners - A_west = ket[row, col] - Ā_west = bra[row, col] - A_east = ket[row, col + 1] - Ā_east = bra[row, col + 1] - - E_northwest = edge(env, NORTH, row - 1, col) - E_northeast = absorb_right( - edge(env, NORTH, row - 1, col + 1), corner(env, NORTHEAST, row - 1, col + 2) - ) - E_east = absorb_right( - edge(env, EAST, row, col + 2), corner(env, SOUTHEAST, row + 1, col + 2) - ) - E_southeast = edge(env, SOUTH, row + 1, col + 1) - E_southwest = absorb_right( - edge(env, SOUTH, row + 1, col), corner(env, SOUTHWEST, row + 1, col - 1) - ) - E_west = absorb_right( - edge(env, WEST, row, col - 1), corner(env, NORTHWEST, row - 1, col - 1) - ) - - @tensor EE_SW[χS χNW DSWt DWt; DSWb DWb] := - E_southwest[χS DSWt DSWb; χSW] * E_west[χSW DWt DWb; χNW] - @tensor EEA_SW[χS DMb dWb χNW DNWb; DSWt DWt] := - EE_SW[χS χNW DSWt DWt; DSWb DWb] * conj(Ā_west[dWb; DNWb DMb DSWb DWb]) - @tensor EEAA_SW[χS DMb dWb dWt DMt; χNW DNWt DNWb] := - EEA_SW[χS DMb dWb χNW DNWb; DSWt DWt] * A_west[dWt; DNWt DMt DSWt DWt] - @tensor EEEAA_W[dWt dWb; χS DMt DMb χN] := - EEAA_SW[χS DMb dWb dWt DMt; χNW DNWt DNWb] * E_northwest[χNW DNWt DNWb; χN] - - @tensor EE_NE[χN χSE DNEt DEt; DNEb DEb] := - E_northeast[χN DNEt DNEb; χNE] * E_east[χNE DEt DEb; χSE] - @tensor EEA_NE[χN DMb dEb χSE DSEb; DNEt DEt] := - EE_NE[χN χSE DNEt DEt; DNEb DEb] * conj(Ā_east[dEb; DNEb DEb DSEb DMb]) - @tensor EEAA_NE[χN DMb dEb dEt DMt; χSE DSEt DSEb] := - EEA_NE[χN DMb dEb χSE DSEb; DNEt DEt] * A_east[dEt; DNEt DEt DSEt DMt] - @tensor EEEAA_E[χS DMt DMb χN; dEt dEb] := - EEAA_NE[χN DMb dEb dEt DMt; χSE DSEt DSEb] * E_southeast[χSE DSEt DSEb; χS] - - @tensor ρ[dWt dEt; dWb dEb] := - EEEAA_W[dWt dWb; χS DMt DMb χN] * EEEAA_E[χS DMt DMb χN; dEt dEb] - - return ρ / str(ρ) -end diff --git a/src/algorithms/expectation_value.jl b/src/algorithms/expectation_value.jl deleted file mode 100644 index a8f5048b8..000000000 --- a/src/algorithms/expectation_value.jl +++ /dev/null @@ -1,324 +0,0 @@ -# Expectation value of a LocalOperator -# ------------------------------------ - -""" - expectation_value(state, O::LocalOperator, env) - expectation_value(bra, O::LocalOperator, ket, env) - -Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state) of a [`LocalOperator`](@ref) `O`. -This can be done either for a PEPS, or alternatively for a density matrix PEPO. -In the latter case the first signature corresponds to a single layer PEPO contraction, while -the second signature yields a bilayer contraction instead. -""" -function MPSKit.expectation_value( - bra::S, O::LocalOperator, ket::S, env - ) where {S <: InfiniteState} - checklattice(bra, O, ket) - term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - return local_expectation_value(inds, bra, operator, ket, env) - end - return sum(term_vals) -end -MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) -function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) - checklattice(state, O) - term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - return local_expectation_value(inds, state, operator, env) - end - return sum(term_vals) -end - - -# Expectation value of an individual local term -# --------------------------------------------- - -""" - local_expectation_value(inds, bra, operator, ket, env) - local_expectation_value(inds, state, operator, env) - -Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation -value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state), where `operator` is the local -term acting on the sites `inds`. - -The implementation is overloaded based on the type of operator to be evaluated -""" -function local_expectation_value end - -# AbstractTensorMap evaluation goes through reduced density matrix -function local_expectation_value(inds, bra, operator::AbstractTensorMap, ket, env) - ρ = reduced_densitymatrix(inds, ket, bra, env) - return trmul(operator, ρ) -end -function local_expectation_value(inds, state, operator::AbstractTensorMap, env) - ρ = reduced_densitymatrix(inds, state, env) - return trmul(operator, ρ) -end - - -# Local patch contractions -# ------------------------ - -@doc """ - reduced_densitymatrix(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) - reduced_densitymatrix(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) - reduced_densitymatrix(inds, state::InfinitePEPO, env) - -Construct the reduced density matrix `ρ` of `|ket⟩⟨bra|`, where both `ket` and `bra` -correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. -Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the -density matrix PEPO `state`. The reduced density matrix is contracted over the virtual -indices surrounding the open physical indices at sites `inds` using the environment `env`, -and is normalized such that `str(ρ) = 1`. - -See also [`str`](@ref). -""" reduced_densitymatrix - -# PEPS case has fast-path specializations -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) - length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) - - if length(inds) == 2 - if inds[2] - inds[1] == CartesianIndex(1, 0) - return reduced_densitymatrix2x1(inds[1], ket, bra, env) - elseif inds[2] - inds[1] == CartesianIndex(0, 1) - return reduced_densitymatrix1x2(inds[1], ket, bra, env) - end - end - - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, (ket, bra), env) -end -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env - ) - size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, state, env) -end -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env - ) - size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) - size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_densitymatrix(static_inds, (ket, bra), env) -end -reduced_densitymatrix(inds, ket::InfinitePEPS, env) = - reduced_densitymatrix(inds, ket, ket, env) - -# handle deprecations of Tuple inds specifications -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, args... - ) where {N}, - reduced_densitymatrix(collect(inds), args...) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, args... - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), args...) -) - -""" - contract_local_operator(inds, O, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) - contract_local_operator(inds, O, ket::InfinitePEPO, bra::InfinitePEPO, env) - contract_local_operator(inds, O, state::InfinitePEPO, env) - -Contract a local operator `O` between `ket` and `bra` states, computing `⟨bra|O|ket⟩`, where -`ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary -legs. Alternatively, contract a local operator `O` with a density matrix PEPO `state`, -computing `tr(O * state)`. `O` is applied to the open physical indices at sites `inds`, and -the result is contracted over the surrounding virtual indices using the environment `env`. -""" -function contract_local_operator( - inds::Vector{CartesianIndex{2}}, O, - ket::InfinitePEPS, bra::InfinitePEPS, env, - ) - static_inds = Tuple(Val.(inds)) - return _contract_local_operator(static_inds, O, (ket, bra), env) -end -function contract_local_operator( - inds::Vector{CartesianIndex{2}}, O, state::InfinitePEPO, env - ) - size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_local_operator(static_inds, O, state, env) -end -function contract_local_operator( - inds::Vector{CartesianIndex{2}}, O, ket::InfinitePEPO, bra::InfinitePEPO, env - ) - size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) - size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_local_operator(static_inds, O, (ket, bra), env) -end -function contract_local_operator(inds::Vector{Tuple{Int, Int}}, O, args...) - return contract_local_operator(CartesianIndex.(inds), O, args...) -end - -Base.@deprecate( - contract_local_operator( - inds::NTuple, args... - ), - contract_local_operator(collect(inds), args...) -) - -""" - contract_local_norm(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) - contract_local_norm(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) - contract_local_norm(inds, state::InfinitePEPO, env) - -Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch -of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a -PEPS with ancillary legs. Alternatively, contract a local norm patch of a density matrix -PEPO `state`, computing a patch of `tr(state)`. The contracted rectangular norm patch is -determined by the open physical indices `inds` and is contracted over the surrounding virtual -indices using the environment `env`. In particular, the patch location is precisely the same -as that of the patch used in [`contract_local_operator`](@ref). -""" -function contract_local_norm( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env - ) - static_inds = Tuple(Val.(inds)) - return _contract_local_norm(static_inds, (ket, bra), env) -end -function contract_local_norm(inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env) - size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_local_norm(static_inds, state, env) -end -function contract_local_norm( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env - ) - size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) - size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) - static_inds = Tuple(Val.(inds)) - return _contract_local_norm(static_inds, (ket, bra), env) -end -function contract_local_norm(inds::Vector{Tuple{Int, Int}}, args...) - return contract_local_norm(CartesianIndex.(inds), args...) -end - -Base.@deprecate( - contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env), - contract_local_norm(collect(inds), ket, bra, env) -) - - -# Expectation value of a local partition function tensor -# ------------------------------------------------------ - -""" - expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) - -Compute the expectation value corresponding to inserting a local tensor(s) `O` at -position `inds` in the partition function `pf` and contracting the whole using a given CTMRG -environment `env`. - -Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O` -should be a rank-4 tensor conforming to the [`PartitionFunctionTensor`](@ref) indexing -convention. -""" -function MPSKit.expectation_value( - pf::InfinitePartitionFunction, - op::Pair{CartesianIndex{2}, <:AbstractTensorMap{T, S, 2, 2}}, - env, - ) where {T, S} - return contract_local_tensor(op[1], op[2], env) / - contract_local_tensor(op[1], pf[op[1]], env) -end -function MPSKit.expectation_value( - pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env - ) - return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) -end - -# Network values -# -------------- - -""" - network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) - -Return the value (per unit cell) of a given contractible network contracted using a given -CTMRG environment. -""" -function network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) - return prod(Iterators.product(axes(network)...)) do (r, c) - return _contract_site((r, c), network, env) * _contract_corners((r, c), env) / - _contract_vertical_edges((r, c), env) / _contract_horizontal_edges((r, c), env) - end -end -network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state), env) - -""" - _contract_site(ind::Tuple{Int,Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) - -Contract around a single site `ind` of a square network using a given CTMRG environment. -""" -function _contract_site(ind::Tuple{Int, Int}, network::InfiniteSquareNetwork, env::CTMRGEnv) - r, c = ind - return _contract_site( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r + 1, c + 1), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), - edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), - network[r, c], - ) -end - -""" - _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract all corners around the south-east at position `ind` of the CTMRG -environment `env`. -""" -function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_corners( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c), - corner(env, SOUTHEAST, r, c), - corner(env, SOUTHWEST, r, c - 1), - ) -end - -""" - _contract_vertical_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract the vertical edges and corners around the east edge at position `ind` of the -CTMRG environment `env`. -""" -function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_vertical_edges( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c), - corner(env, SOUTHEAST, r + 1, c), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, EAST, r, c), - edge(env, WEST, r, c - 1), - ) -end - -""" - _contract_horizontal_edges(ind::Tuple{Int,Int}, env::CTMRGEnv) - -Contract the horizontal edges and corners around the south edge at position `ind` of the -CTMRG environment `env`. -""" -function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) - r, c = ind - return _contract_horizontal_edges( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r, c + 1), - corner(env, SOUTHWEST, r, c - 1), - edge(env, NORTH, r - 1, c), - edge(env, SOUTH, r, c), - ) -end diff --git a/src/algorithms/correlator_adapters.jl b/src/algorithms/expectation_value/correlator_adapters.jl similarity index 100% rename from src/algorithms/correlator_adapters.jl rename to src/algorithms/expectation_value/correlator_adapters.jl diff --git a/src/algorithms/correlators.jl b/src/algorithms/expectation_value/correlators.jl similarity index 100% rename from src/algorithms/correlators.jl rename to src/algorithms/expectation_value/correlators.jl diff --git a/src/algorithms/expectation_value/expectation_value.jl b/src/algorithms/expectation_value/expectation_value.jl new file mode 100644 index 000000000..d2ea498f6 --- /dev/null +++ b/src/algorithms/expectation_value/expectation_value.jl @@ -0,0 +1,83 @@ +# Expectation value of a LocalOperator +# ------------------------------------ + +""" + expectation_value(state, O::LocalOperator, env) + expectation_value(bra, O::LocalOperator, ket, env) + +Compute the expectation value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state) of a [`LocalOperator`](@ref) `O`. +This can be done either for a PEPS, or alternatively for a density matrix PEPO. +In the latter case the first signature corresponds to a single layer PEPO contraction, while +the second signature yields a bilayer contraction instead. +""" +function MPSKit.expectation_value( + bra::S, O::LocalOperator, ket::S, env + ) where {S <: InfiniteState} + checklattice(bra, O, ket) + term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly + return local_expectation_value(inds, bra, operator, ket, env) + end + return sum(term_vals) +end +MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env) = expectation_value(peps, O, peps, env) +function MPSKit.expectation_value(state::InfinitePEPO, O::LocalOperator, env) + checklattice(state, O) + term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly + return local_expectation_value(inds, state, operator, env) + end + return sum(term_vals) +end + + +# Expectation value of an individual local term +# --------------------------------------------- + +""" + local_expectation_value(inds, bra, operator, ket, env) + local_expectation_value(inds, state, operator, env) + +Compute the contribution of a single term of a [`LocalOperator`](@ref) to the expectation +value ⟨bra|O|ket⟩ / ⟨bra|ket⟩ or tr(O * state) / tr(state), where `operator` is the local +term acting on the sites `inds`. + +The implementation is overloaded based on the type of operator to be evaluated +""" +function local_expectation_value end + +# AbstractTensorMap evaluation goes through reduced density matrix +function local_expectation_value(inds, bra, operator::AbstractTensorMap, ket, env) + ρ = reduced_densitymatrix(inds, ket, bra, env) + return trmul(operator, ρ) +end +function local_expectation_value(inds, state, operator::AbstractTensorMap, env) + ρ = reduced_densitymatrix(inds, state, env) + return trmul(operator, ρ) +end + +# Expectation value of a local partition function tensor +# ------------------------------------------------------ + +""" + expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv) + +Compute the expectation value corresponding to inserting a local tensor(s) `O` at +position `inds` in the partition function `pf` and contracting the whole using a given CTMRG +environment `env`. + +Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O` +should be a rank-4 tensor conforming to the [`PartitionFunctionTensor`](@ref) indexing +convention. +""" +function MPSKit.expectation_value( + pf::InfinitePartitionFunction, + op::Pair{CartesianIndex{2}, <:AbstractTensorMap{T, S, 2, 2}}, + env, + ) where {T, S} + return contract_local_tensor(op[1], op[2], env) / + contract_local_tensor(op[1], pf[op[1]], env) +end +function MPSKit.expectation_value( + pf::InfinitePartitionFunction, op::Pair{Tuple{Int, Int}}, env + ) + return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) +end diff --git a/src/algorithms/expectation_value/network_value.jl b/src/algorithms/expectation_value/network_value.jl new file mode 100644 index 000000000..0bac1e033 --- /dev/null +++ b/src/algorithms/expectation_value/network_value.jl @@ -0,0 +1,73 @@ +# Network values +# -------------- + +""" + network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) + +Return the value (per unit cell) of a given contractible network contracted using a given +CTMRG environment. +""" +function network_value(network::InfiniteSquareNetwork, env::CTMRGEnv) + return prod(Iterators.product(axes(network)...)) do (r, c) + return _contract_site((r, c), network, env) * _contract_corners((r, c), env) / + _contract_vertical_edges((r, c), env) / _contract_horizontal_edges((r, c), env) + end +end +network_value(state, env::CTMRGEnv) = network_value(InfiniteSquareNetwork(state), env) + +function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) + return network_value(InfiniteSquareNetwork(peps), env) +end + +# Local tensor insertions +# ----------------------- + +""" + contract_local_tensor(inds, O::PFTensor, env) + +Contract a local tensor `O` inserted into a partition function `pf` at position `inds`, +using the environment `env`. +""" +function contract_local_tensor( + inds::Tuple{Int, Int}, O::PFTensor, env::CTMRGEnv{C, <:CTMRG_PF_EdgeTensor} + ) where {C} + r, c = inds + return _contract_site( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), + O, + ) +end + +""" + contract_local_tensor(inds, O::PEPOTensor, network, env) + +Contract a local tensor `O` inserted into the PEPO of a given `network` at position `inds`, +using the environment `env`. +""" +function contract_local_tensor( + ind::Tuple{Int, Int, Int}, + O::PEPOTensor, + network::InfiniteSquareNetwork{<:PEPOSandwich}, + env::CTMRGEnv, + ) + r, c, h = ind + sandwich´ = Base.setindex(network[r, c], O, h + 2) + return _contract_site( + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), + sandwich´, + ) +end + +function contract_local_tensor(inds::CartesianIndex, O::AbstractTensorMap, env::CTMRGEnv) + return contract_local_tensor(Tuple(inds), O, env) +end diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl new file mode 100644 index 000000000..01d61e8c4 --- /dev/null +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -0,0 +1,87 @@ +# Direct local patch contractions +# ------------------------------- + +""" + contract_local_operator(inds, O, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_operator(inds, O, ket::InfinitePEPO, bra::InfinitePEPO, env) + contract_local_operator(inds, O, state::InfinitePEPO, env) + +Contract a local operator `O` between `ket` and `bra` states, computing `⟨bra|O|ket⟩`, where +`ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary +legs. Alternatively, contract a local operator `O` with a density matrix PEPO `state`, +computing `tr(O * state)`. `O` is applied to the open physical indices at sites `inds`, and +the result is contracted over the surrounding virtual indices using the environment `env`. +""" +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, + ket::InfinitePEPS, bra::InfinitePEPS, env, + ) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, state::InfinitePEPO, env + ) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, state, env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator(inds::Vector{Tuple{Int, Int}}, O, args...) + return contract_local_operator(CartesianIndex.(inds), O, args...) +end + +Base.@deprecate( + contract_local_operator( + inds::NTuple, args... + ), + contract_local_operator(collect(inds), args...) +) + +""" + contract_local_norm(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_norm(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) + contract_local_norm(inds, state::InfinitePEPO, env) + +Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch +of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a +PEPS with ancillary legs. Alternatively, contract a local norm patch of a density matrix +PEPO `state`, computing a patch of `tr(state)`. The contracted rectangular norm patch is +determined by the open physical indices `inds` and is contracted over the surrounding virtual +indices using the environment `env`. In particular, the patch location is precisely the same +as that of the patch used in [`contract_local_operator`](@ref). +""" +function contract_local_norm( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env + ) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, (ket, bra), env) +end +function contract_local_norm(inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, state, env) +end +function contract_local_norm( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_local_norm(static_inds, (ket, bra), env) +end +function contract_local_norm(inds::Vector{Tuple{Int, Int}}, args...) + return contract_local_norm(CartesianIndex.(inds), args...) +end + +Base.@deprecate( + contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env), + contract_local_norm(collect(inds), ket, bra, env) +) diff --git a/src/algorithms/expectation_value/reduced_densitymatrix.jl b/src/algorithms/expectation_value/reduced_densitymatrix.jl new file mode 100644 index 000000000..6bb5fe6cd --- /dev/null +++ b/src/algorithms/expectation_value/reduced_densitymatrix.jl @@ -0,0 +1,109 @@ +# Reduced density matrices +# ------------------------ + +@doc """ + reduced_densitymatrix(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + reduced_densitymatrix(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) + reduced_densitymatrix(inds, state::InfinitePEPO, env) + +Construct the reduced density matrix `ρ` of `|ket⟩⟨bra|`, where both `ket` and `bra` +correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. +Alternatively, construct the reduced density matrix `ρ` of a mixed state specified by the +density matrix PEPO `state`. The reduced density matrix is contracted over the virtual +indices surrounding the open physical indices at sites `inds` using the environment `env`, +and is normalized such that `str(ρ) = 1`. + +See also [`str`](@ref). +""" reduced_densitymatrix + +# PEPS case has fast-path specializations +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == 1 && return reduced_densitymatrix1x1(only(inds), ket, bra, env) + + if length(inds) == 2 + if inds[2] - inds[1] == CartesianIndex(1, 0) + return reduced_densitymatrix2x1(inds[1], ket, bra, env) + elseif inds[2] - inds[1] == CartesianIndex(0, 1) + return reduced_densitymatrix1x2(inds[1], ket, bra, env) + end + end + + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, (ket, bra), env) +end +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, state::InfinitePEPO, env + ) + size(state, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, state, env) +end +function reduced_densitymatrix( + inds::Vector{CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env + ) + size(ket) == size(bra) || throw(DimensionMismatch("incompatible bra and ket dimensions")) + size(ket, 3) == 1 || throw(DimensionMismatch("only single-layer densitymatrices are supported")) + static_inds = Tuple(Val.(inds)) + return _contract_densitymatrix(static_inds, (ket, bra), env) +end +reduced_densitymatrix(inds, ket::InfinitePEPS, env) = + reduced_densitymatrix(inds, ket, ket, env) + +# handle deprecations of Tuple inds specifications +Base.@deprecate( + reduced_densitymatrix( + inds::NTuple{N, CartesianIndex{2}}, args... + ) where {N}, + reduced_densitymatrix(collect(inds), args...) +) +Base.@deprecate( + reduced_densitymatrix( + inds::NTuple{N, Tuple{Int, Int}}, args... + ) where {N}, + reduced_densitymatrix(collect(CartesianIndex.(inds)), args...) +) + + +# Fixed-size fast paths +# --------------------- +# +# These carry the docstrings and reject unsupported environments; the implementations live in +# `algorithms/contractions/local_patch/densitymatrix/`. + +""" + reduced_densitymatrix1x1(ind, ket, bra, env) + +Construct the reduced density matrix of `|ket⟩⟨bra|` on the single site `ind`, using an +optimized contraction for the environment `env`. +""" +reduced_densitymatrix1x1(ind, ket, bra, env) = throw( + ArgumentError( + "No 1x1 reduced density matrix contraction defined for environments of type $(typeof(env))." + ) +) + +""" + reduced_densitymatrix2x1(ind, ket, bra, env) + +Construct the reduced density matrix of `|ket⟩⟨bra|` on the vertical pair of sites starting +at `ind`, using an optimized contraction for the environment `env`. +""" +reduced_densitymatrix2x1(ind, ket, bra, env) = throw( + ArgumentError( + "No 2x1 reduced density matrix contraction defined for environments of type $(typeof(env))." + ) +) + +""" + reduced_densitymatrix1x2(ind, ket, bra, env) + +Construct the reduced density matrix of `|ket⟩⟨bra|` on the horizontal pair of sites starting +at `ind`, using an optimized contraction for the environment `env`. +""" +reduced_densitymatrix1x2(ind, ket, bra, env) = throw( + ArgumentError( + "No 1x2 reduced density matrix contraction defined for environments of type $(typeof(env))." + ) +) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 2b93d85c0..533fac679 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -1,7 +1,3 @@ -function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) - return network_value(InfiniteSquareNetwork(peps), env) -end - """ edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} @@ -129,55 +125,3 @@ function product_peps(peps_args...; unitcell = (1, 1), noise_amp = 1.0e-2, state ψ = prod_peps + noise_amp * noise_peps return ψ / norm(ψ) end - -# Contract local tensors - -""" - contract_local_tensor(inds, O::PFTensor, env) - -Contract a local tensor `O` inserted into a partition function `pf` at position `inds`, -using the environment `env`. -""" -function contract_local_tensor( - inds::Tuple{Int, Int}, O::PFTensor, env::CTMRGEnv{C, <:CTMRG_PF_EdgeTensor} - ) where {C} - r, c = inds - return _contract_site( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r + 1, c + 1), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), - edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), - O, - ) -end - -""" - contract_local_tensor(inds, O::PEPOTensor, network, env) - -Contract a local tensor `O` inserted into the PEPO of a given `network` at position `inds`, -using the environment `env`. -""" -function contract_local_tensor( - ind::Tuple{Int, Int, Int}, - O::PEPOTensor, - network::InfiniteSquareNetwork{<:PEPOSandwich}, - env::CTMRGEnv, - ) - r, c, h = ind - sandwich´ = Base.setindex(network[r, c], O, h + 2) - return _contract_site( - corner(env, NORTHWEST, r - 1, c - 1), - corner(env, NORTHEAST, r - 1, c + 1), - corner(env, SOUTHEAST, r + 1, c + 1), - corner(env, SOUTHWEST, r + 1, c - 1), - edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), - edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), - sandwich´, - ) -end - -function contract_local_tensor(inds::CartesianIndex, O::AbstractTensorMap, env::CTMRGEnv) - return contract_local_tensor(Tuple(inds), O, env) -end diff --git a/src/utility/tensor_traces.jl b/src/utility/tensor_traces.jl new file mode 100644 index 000000000..7f52da297 --- /dev/null +++ b/src/utility/tensor_traces.jl @@ -0,0 +1,27 @@ +# Tensor traces +# ------------- + +""" + str(t) + +Fermionic supertrace by using `@tensor`. +""" +str(t::AbstractTensorMap) = _str(BraidingStyle(sectortype(t)), t) +_str(::Bosonic, t::AbstractTensorMap) = tr(t) +@generated function _str(::Fermionic, t::AbstractTensorMap{<:Any, <:Any, N, N}) where {N} + tex = tensorexpr(:t, ntuple(identity, N), ntuple(identity, N)) + return macroexpand(@__MODULE__, :(@tensor $tex)) +end + +""" + trmul(H, ρ) + +Compute `tr(H * ρ)` without forming `H * ρ`. +""" +@generated function trmul( + H::AbstractTensorMap{<:Any, S, N, N}, ρ::AbstractTensorMap{<:Any, S, N, N} + ) where {S, N} + Hex = tensorexpr(:H, ntuple(identity, N), ntuple(i -> i + N, N)) + ρex = tensorexpr(:ρ, ntuple(i -> i + N, N), ntuple(identity, N)) + return macroexpand(@__MODULE__, :(@tensor $Hex * $ρex)) +end diff --git a/src/utility/util.jl b/src/utility/util.jl index b79e1644b..603aec675 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -62,96 +62,8 @@ function absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorM return U * sqrt_S, sqrt_S * V end -""" - absorb_left( - A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - -Absorb a matrix `C` into the left of a tensor map `A` by contracting the first leg of `A` -with the second leg of `C`. -""" -function absorb_left( - A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - pC = (codomainind(C), domainind(C)) - pA = ((codomainind(A)[1],), (codomainind(A)[2:end]..., domainind(A)...)) - pCA = (codomainind(A), domainind(A)) - return tensorcontract(C, pC, false, A, pA, false, pCA) -end -function absorb_left( - P::AbstractTensorMap{<:Any, S, 1, N}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S, N} - return twistnondual(C, 2) * P -end - -""" - absorb_right( - A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - -Absorb a matrix `C` into the right of a tensor map `A` by contracting the last leg of `A` -with the first leg of `C`. -""" -function absorb_right( - A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - pA = ((codomainind(A)..., domainind(A)[2:end]...), (domainind(A)[1],)) - pC = (codomainind(C), domainind(C)) - pAC = (codomainind(A), (domainind(A)[end], domainind(A)[1:(end - 1)]...)) - return tensorcontract(A, pA, false, C, pC, false, pAC) -end -function absorb_right( - E::AbstractTensorMap{<:Any, S, N, 1}, C::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S, N} - return E * twistdual(C, 1) -end - -""" - absorb_left_right( - A::AbstractTensorMap{<:Any, S}, - CL::AbstractTensorMap{<:Any, S, 1, 1}, - CR::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - -Absorb matrices `CL` and `CR` into the left and right of a tensor map `A` by contracting the -first leg of `A` with the second leg of `CL` and the last leg of `A` with the first leg of -`CR`. -""" -function absorb_left_right( - A::AbstractTensorMap{<:Any, S}, - CL::AbstractTensorMap{<:Any, S, 1, 1}, - CR::AbstractTensorMap{<:Any, S, 1, 1} - ) where {S} - return absorb_right(absorb_left(A, CL), CR) -end - _fliptwist_s(s::DiagonalTensorMap) = twist!(DiagonalTensorMap(flip(s, 1:2)), 1) -""" - str(t) - -Fermionic supertrace by using `@tensor`. -""" -str(t::AbstractTensorMap) = _str(BraidingStyle(sectortype(t)), t) -_str(::Bosonic, t::AbstractTensorMap) = tr(t) -@generated function _str(::Fermionic, t::AbstractTensorMap{<:Any, <:Any, N, N}) where {N} - tex = tensorexpr(:t, ntuple(identity, N), ntuple(identity, N)) - return macroexpand(@__MODULE__, :(@tensor $tex)) -end - -""" - trmul(H, ρ) - -Compute `tr(H * ρ)` without forming `H * ρ`. -""" -@generated function trmul( - H::AbstractTensorMap{<:Any, S, N, N}, ρ::AbstractTensorMap{<:Any, S, N, N} - ) where {S, N} - Hex = tensorexpr(:H, ntuple(identity, N), ntuple(i -> i + N, N)) - ρex = tensorexpr(:ρ, ntuple(i -> i + N, N), ntuple(identity, N)) - return macroexpand(@__MODULE__, :(@tensor $Hex * $ρex)) -end - # Check whether diagonals contain degenerate values up to absolute or relative tolerance function is_degenerate_spectrum( S; atol::Real = 0, rtol::Real = atol > 0 ? 0 : sqrt(eps(scalartype(S))) From 42c39ab84afffa96ca60a31cb7e6c7a6868e3a5f Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 18 Sep 2026 11:08:13 +0200 Subject: [PATCH 25/32] Make `absorb` docstrings more accurate --- src/algorithms/contractions/absorb.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/algorithms/contractions/absorb.jl b/src/algorithms/contractions/absorb.jl index b0387ef4c..a7ac8b66e 100644 --- a/src/algorithms/contractions/absorb.jl +++ b/src/algorithms/contractions/absorb.jl @@ -6,8 +6,9 @@ A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the left of a tensor map `A` by contracting the first leg of `A` -with the second leg of `C`. +Absorb a matrix `C` into the left of a tensor map `A` by contracting the first index in the codomain of `A` +with the (only) index in the domain of `C`. This can be interpreted as contracting the first +leg of `A` with the last leg of `C`. """ function absorb_left( A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} @@ -28,7 +29,9 @@ end A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the right of a tensor map `A` by contracting the last leg of `A` +Absorb a matrix `C` into the right of a tensor map `A` by contracting the first index in +the domain of `A` with the (only) index in the codomain of `C`. In the case where `A` has +only one space in its domain, this can be interpreted as contracting the last leg of `A` with the first leg of `C`. """ function absorb_right( From 741bd39b1fa9ef68c9b12eb5299612620d03c1b3 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 18 Sep 2026 11:08:34 +0200 Subject: [PATCH 26/32] Be explicit about supertrace usage --- src/utility/tensor_traces.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/utility/tensor_traces.jl b/src/utility/tensor_traces.jl index 7f52da297..f99409180 100644 --- a/src/utility/tensor_traces.jl +++ b/src/utility/tensor_traces.jl @@ -16,7 +16,9 @@ end """ trmul(H, ρ) -Compute `tr(H * ρ)` without forming `H * ρ`. +Compute `str(H * ρ)` without forming `H * ρ`. + +See also [`str`](@ref). """ @generated function trmul( H::AbstractTensorMap{<:Any, S, N, N}, ρ::AbstractTensorMap{<:Any, S, N, N} From 6c2475803b08e6eeda62cc20f81f2fc45e2c0a36 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 18 Sep 2026 11:09:45 +0200 Subject: [PATCH 27/32] Move basic label generators up into `util` --- src/PEPSKit.jl | 1 + .../contractions/local_patch/expr_utils.jl | 10 ---------- src/utility/contraction_labels.jl | 15 +++++++++++++++ 3 files changed, 16 insertions(+), 10 deletions(-) create mode 100644 src/utility/contraction_labels.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 44c7c7f50..827a16689 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -52,6 +52,7 @@ using DocStringExtensions include("Defaults.jl") # Include first to allow for docstring interpolation with Defaults values include("utility/util.jl") +include("utility/contraction_labels.jl") include("utility/tensor_traces.jl") include("utility/indexing.jl") include("utility/diffable_threads.jl") diff --git a/src/algorithms/contractions/local_patch/expr_utils.jl b/src/algorithms/contractions/local_patch/expr_utils.jl index 6c0125614..37e96457a 100644 --- a/src/algorithms/contractions/local_patch/expr_utils.jl +++ b/src/algorithms/contractions/local_patch/expr_utils.jl @@ -1,16 +1,6 @@ # Contraction expression utilities for local patches # -------------------------------------------------- -# Contraction label helpers -# ------------------------- - -function tensorlabel(args...) - return Symbol(ntuple(i -> iseven(i) ? :_ : args[(i + 1) >> 1], 2 * length(args) - 1)...) -end -envlabel(args...) = tensorlabel(:χ, args...) -virtuallabel(args...) = tensorlabel(:D, args...) -physicallabel(args...) = tensorlabel(:d, args...) - """ $(SIGNATURES) diff --git a/src/utility/contraction_labels.jl b/src/utility/contraction_labels.jl new file mode 100644 index 000000000..2ae3020a6 --- /dev/null +++ b/src/utility/contraction_labels.jl @@ -0,0 +1,15 @@ +# Contraction index label helpers +# ------------------------------- +# +# These are shared by every contraction expression builder in the package (the CTMRG +# contractions as well as the local-patch ones). Since most of those builders are called from +# the generators of `@generated` functions, the helpers have to be defined *before* any of +# them: a generator can only call methods that already existed when the generated function was +# defined. + +function tensorlabel(args...) + return Symbol(ntuple(i -> iseven(i) ? :_ : args[(i + 1) >> 1], 2 * length(args) - 1)...) +end +envlabel(args...) = tensorlabel(:χ, args...) +virtuallabel(args...) = tensorlabel(:D, args...) +physicallabel(args...) = tensorlabel(:d, args...) From da753321cd4018a038bca5e8d4683e8ea8de6cae Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 18 Sep 2026 11:12:07 +0200 Subject: [PATCH 28/32] Remove non-existing default value from docstring --- src/algorithms/expectation_value/patch_contractions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl index 01d61e8c4..47fd7fc80 100644 --- a/src/algorithms/expectation_value/patch_contractions.jl +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -2,7 +2,7 @@ # ------------------------------- """ - contract_local_operator(inds, O, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_operator(inds, O, ket::InfinitePEPS, bra::InfinitePEPS, env) contract_local_operator(inds, O, ket::InfinitePEPO, bra::InfinitePEPO, env) contract_local_operator(inds, O, state::InfinitePEPO, env) From c5a0ead4a997246a7d2b6397675de4b850354e29 Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 18 Sep 2026 11:23:04 +0200 Subject: [PATCH 29/32] Always point to generic `_contract_densitymatrix` fallback --- .../reduced_densitymatrix.jl | 38 +++++++++---------- test/toolbox/densitymatrices.jl | 28 ++++++++++++++ 2 files changed, 46 insertions(+), 20 deletions(-) diff --git a/src/algorithms/expectation_value/reduced_densitymatrix.jl b/src/algorithms/expectation_value/reduced_densitymatrix.jl index 6bb5fe6cd..2767c2bc0 100644 --- a/src/algorithms/expectation_value/reduced_densitymatrix.jl +++ b/src/algorithms/expectation_value/reduced_densitymatrix.jl @@ -69,41 +69,39 @@ Base.@deprecate( # Fixed-size fast paths # --------------------- # -# These carry the docstrings and reject unsupported environments; the implementations live in -# `algorithms/contractions/local_patch/densitymatrix/`. +# These carry the docstrings and provide the generic fallback; the optimized implementations +# live in `algorithms/contractions/local_patch/densitymatrix/`. +# +# The fallbacks defer to the generic patch contraction instead of erroring, so that an +# environment only needs to define the specializations for the patch shapes where a +# hand-optimized contraction actually pays off. """ reduced_densitymatrix1x1(ind, ket, bra, env) Construct the reduced density matrix of `|ket⟩⟨bra|` on the single site `ind`, using an -optimized contraction for the environment `env`. +optimized contraction for the environment `env`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. """ -reduced_densitymatrix1x1(ind, ket, bra, env) = throw( - ArgumentError( - "No 1x1 reduced density matrix contraction defined for environments of type $(typeof(env))." - ) -) +reduced_densitymatrix1x1(ind, ket, bra, env) = + _contract_densitymatrix((Val(ind),), (ket, bra), env) """ reduced_densitymatrix2x1(ind, ket, bra, env) Construct the reduced density matrix of `|ket⟩⟨bra|` on the vertical pair of sites starting -at `ind`, using an optimized contraction for the environment `env`. +at `ind`, using an optimized contraction for the environment `env`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. """ -reduced_densitymatrix2x1(ind, ket, bra, env) = throw( - ArgumentError( - "No 2x1 reduced density matrix contraction defined for environments of type $(typeof(env))." - ) -) +reduced_densitymatrix2x1(ind, ket, bra, env) = + _contract_densitymatrix((Val(ind), Val(ind + CartesianIndex(1, 0))), (ket, bra), env) """ reduced_densitymatrix1x2(ind, ket, bra, env) Construct the reduced density matrix of `|ket⟩⟨bra|` on the horizontal pair of sites starting -at `ind`, using an optimized contraction for the environment `env`. +at `ind`, using an optimized contraction for the environment `env`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. """ -reduced_densitymatrix1x2(ind, ket, bra, env) = throw( - ArgumentError( - "No 1x2 reduced density matrix contraction defined for environments of type $(typeof(env))." - ) -) +reduced_densitymatrix1x2(ind, ket, bra, env) = + _contract_densitymatrix((Val(ind), Val(ind + CartesianIndex(0, 1))), (ket, bra), env) diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl index 8c3307a88..0c7e7a4f4 100644 --- a/test/toolbox/densitymatrices.jl +++ b/test/toolbox/densitymatrices.jl @@ -101,3 +101,31 @@ end inds = Tuple(Val.([CartesianIndex(1, 1)])) @test_throws ArgumentError PEPSKit._contract_densitymatrix(inds, (ρ, ρ, ρ), env) end + +@testset "Fixed-size fast paths agree with the generic fallback ($I)" for I in keys(ds) + d, D, χ = ds[I], Ds[I], χs[I] + peps = InfinitePEPS(d, D; unitcell = (2, 2)) + env = CTMRGEnv(peps, χ) + + # environments without a hand-optimized contraction of a given shape fall through to + # `_contract_densitymatrix`, which should agree with the specializations that do exist + for ind in CartesianIndices((2, 2)) + ρ_fast = reduced_densitymatrix([ind], peps, env) + ρ_generic = invoke( + PEPSKit.reduced_densitymatrix1x1, Tuple{Any, Any, Any, Any}, ind, peps, peps, env + ) + @test ρ_fast ≈ ρ_generic + + ρ_fast = reduced_densitymatrix([ind, ind + CartesianIndex(1, 0)], peps, env) + ρ_generic = invoke( + PEPSKit.reduced_densitymatrix2x1, Tuple{Any, Any, Any, Any}, ind, peps, peps, env + ) + @test ρ_fast ≈ ρ_generic + + ρ_fast = reduced_densitymatrix([ind, ind + CartesianIndex(0, 1)], peps, env) + ρ_generic = invoke( + PEPSKit.reduced_densitymatrix1x2, Tuple{Any, Any, Any, Any}, ind, peps, peps, env + ) + @test ρ_fast ≈ ρ_generic + end +end From 80a3efbb9f347b621126054053bae9a481ca91dc Mon Sep 17 00:00:00 2001 From: leburgel Date: Sun, 20 Sep 2026 11:36:13 +0200 Subject: [PATCH 30/32] Be more explicit about supertrace again --- .../expectation_value/patch_contractions.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl index 47fd7fc80..aa03968ae 100644 --- a/src/algorithms/expectation_value/patch_contractions.jl +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -52,11 +52,15 @@ Base.@deprecate( Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a -PEPS with ancillary legs. Alternatively, contract a local norm patch of a density matrix -PEPO `state`, computing a patch of `tr(state)`. The contracted rectangular norm patch is -determined by the open physical indices `inds` and is contracted over the surrounding virtual -indices using the environment `env`. In particular, the patch location is precisely the same -as that of the patch used in [`contract_local_operator`](@ref). +PEPS with ancillary legs. +Alternatively, contract a local norm patch of a density matrix PEPO `state`, computing a +patch of `str(state)`, where [`str`](@ref) is the fermionic supertrace. + +The contracted rectangular norm patch is determined by the open physical indices `inds` and +is contracted over the surrounding virtual indices using the environment `env`. In +particular, the patch location is precisely the same as that of the patch used in +[`contract_local_operator`](@ref). + """ function contract_local_norm( inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env From 49496f4793f41a9f533de40671fd3a427bf98386 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 21 Sep 2026 10:14:22 +0800 Subject: [PATCH 31/32] Docstring cleanups --- src/algorithms/contractions/absorb.jl | 14 ++++---------- .../expectation_value/patch_contractions.jl | 6 ++---- 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/src/algorithms/contractions/absorb.jl b/src/algorithms/contractions/absorb.jl index a7ac8b66e..5101bd4bc 100644 --- a/src/algorithms/contractions/absorb.jl +++ b/src/algorithms/contractions/absorb.jl @@ -6,9 +6,7 @@ A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the left of a tensor map `A` by contracting the first index in the codomain of `A` -with the (only) index in the domain of `C`. This can be interpreted as contracting the first -leg of `A` with the last leg of `C`. +Absorb a matrix `C` into the left of a tensor map `A` by contracting the first index of `A` with the last index of `C`. """ function absorb_left( A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} @@ -29,10 +27,8 @@ end A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb a matrix `C` into the right of a tensor map `A` by contracting the first index in -the domain of `A` with the (only) index in the codomain of `C`. In the case where `A` has -only one space in its domain, this can be interpreted as contracting the last leg of `A` -with the first leg of `C`. +Absorb a matrix `C` into the right of a tensor map `A` by contracting the first index in the domain of `A` with the first index of `C`. +In the case where `A` has only one space in its domain, this can be interpreted as contracting the last leg of `A` with the first leg of `C`. """ function absorb_right( A::AbstractTensorMap{<:Any, S}, C::AbstractTensorMap{<:Any, S, 1, 1} @@ -55,9 +51,7 @@ end CR::AbstractTensorMap{<:Any, S, 1, 1} ) where {S} -Absorb matrices `CL` and `CR` into the left and right of a tensor map `A` by contracting the -first leg of `A` with the second leg of `CL` and the last leg of `A` with the first leg of -`CR`. +Absorb matrices `CL` and `CR` into the left and right of a tensor map `A` by contracting the first leg of `A` with the second leg of `CL` and the first domain index of `A` with the first leg of `CR`. """ function absorb_left_right( A::AbstractTensorMap{<:Any, S}, diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl index aa03968ae..324d2ebb5 100644 --- a/src/algorithms/expectation_value/patch_contractions.jl +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -46,21 +46,19 @@ Base.@deprecate( ) """ - contract_local_norm(inds, ket::InfinitePEPS, bra::InfinitePEPS = ket, env) + contract_local_norm(inds, ket::InfinitePEPS, bra::InfinitePEPS, env) contract_local_norm(inds, ket::InfinitePEPO, bra::InfinitePEPO, env) contract_local_norm(inds, state::InfinitePEPO, env) Contract a local norm corresponding to the overlap `ket` and `bra` states, computing a patch of `⟨bra|ket⟩`, where `ket` and `bra` correspond to either a PEPS or a PEPO representing a PEPS with ancillary legs. -Alternatively, contract a local norm patch of a density matrix PEPO `state`, computing a -patch of `str(state)`, where [`str`](@ref) is the fermionic supertrace. +Alternatively, contract a local norm patch of a density matrix PEPO `state`, computing a patch of `tr(state)`. The contracted rectangular norm patch is determined by the open physical indices `inds` and is contracted over the surrounding virtual indices using the environment `env`. In particular, the patch location is precisely the same as that of the patch used in [`contract_local_operator`](@ref). - """ function contract_local_norm( inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env From 814667f70a5fe1ad8cb695378e7b414675c0513c Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 21 Sep 2026 10:14:51 +0800 Subject: [PATCH 32/32] Latest runic v1.10 formatting --- src/operators/transfermatrix.jl | 8 ++++---- test/ctmrg/unitcell.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl index a033817a4..76cfc54e3 100644 --- a/src/operators/transfermatrix.jl +++ b/src/operators/transfermatrix.jl @@ -150,10 +150,10 @@ function initialize_mps( return InfiniteMPS( [ f( - T, - virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)), - virtualspaces[mod1(i, end)], - ) for i in 1:length(O) + T, + virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)), + virtualspaces[mod1(i, end)], + ) for i in 1:length(O) ] ) end diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index f87548469..577d65205 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -30,9 +30,9 @@ function test_unitcell( Pspaces, [ (c,) => randn( - scalartype(peps), - Pspaces[c], Pspaces[c], - ) for c in CartesianIndices(unitcell) + scalartype(peps), + Pspaces[c], Pspaces[c], + ) for c in CartesianIndices(unitcell) ]..., ) @test expectation_value(peps, random_op, env) isa Number