diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index be6dec4d7..827a16689 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -52,6 +52,8 @@ 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") include("utility/twistdual.jl") @@ -95,15 +97,21 @@ 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") +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") @@ -145,9 +153,13 @@ include("algorithms/bp/beliefpropagation.jl") include("algorithms/bp/gaugefix.jl") include("algorithms/transfermatrix.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..5101bd4bc --- /dev/null +++ b/src/algorithms/contractions/absorb.jl @@ -0,0 +1,62 @@ +# 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 index of `A` with the last index 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 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} + ) 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 first domain index 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_contractions.jl b/src/algorithms/contractions/bp_contractions.jl deleted file mode 100644 index b751cbd75..000000000 --- a/src/algorithms/contractions/bp_contractions.jl +++ /dev/null @@ -1,257 +0,0 @@ -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 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( - 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) - - 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) - elseif ind_relative == CartesianIndex(0, 1) - return contract_local_operator1x2(inds[1], O, 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) - - 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, - ) - 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][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] - - 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_operator2x1( - coord::CartesianIndex{2}, - O::AbstractTensorMap{<:Any, <:Any, 2, 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][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]) * - M_north[DNt; DNb] * - M_northeast[DNEt; DNEb] * - M_southeast[DSEt; DSEb] * - M_south[DSb; DSt] * - M_southwest[DSWb; DSWt] * - M_northwest[DNWb; DNWt] * - O[dNb dSb; dNt dSt] -end - -function contract_local_operator1x2( - coord::CartesianIndex{2}, - O::AbstractTensorMap{<:Any, <:Any, 2, 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[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] * - M_east[DEt; DEb] * - M_southeast[DSEb; DSEt] * - M_southwest[DSWb; DSWt] * - O[dWb dEb; dWt dEt] - 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/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/characteristic_equations.jl b/src/algorithms/contractions/ctmrg/characteristic_equations.jl index fa852187f..c554f95e1 100644 --- a/src/algorithms/contractions/ctmrg/characteristic_equations.jl +++ b/src/algorithms/contractions/ctmrg/characteristic_equations.jl @@ -317,30 +317,6 @@ function eachcoordinate(tensor_unitcell::Array{<:AbstractTensorMap, 3}) return collect(Iterators.product(axes(tensor_unitcell)...)) 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..3961cc061 --- /dev/null +++ b/src/algorithms/contractions/ctmrg/network_value.jl @@ -0,0 +1,201 @@ +## 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 + +""" + _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, + 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 + +""" + _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, + 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_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, + 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 + +""" + _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/local_patch/densitymatrix/bp.jl b/src/algorithms/contractions/local_patch/densitymatrix/bp.jl new file mode 100644 index 000000000..16c087b52 --- /dev/null +++ b/src/algorithms/contractions/local_patch/densitymatrix/bp.jl @@ -0,0 +1,91 @@ +# Belief Propagation reduced density matrices +# ------------------------------------------- + +# 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) + return throw( + ArgumentError( + "Cannot contract a $(_patch_shape_string(sites)) patch using a `BPEnv`; + only 1x1, 2x1, 1x2 patches are supported." + ) + ) +end + +function reduced_densitymatrix1x1( + 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] + + @autoopt @tensor ρ[dt; db] := + ket[row, col][dt; DNt DEt DSt DWt] * + 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 ρ / str(ρ) +end + +function reduced_densitymatrix2x1( + 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] + + @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]) * + M_north[DNt; DNb] * + M_northeast[DNEt; DNEb] * + M_southeast[DSEt; DSEb] * + M_south[DSb; DSt] * + M_southwest[DSWb; DSWt] * + M_northwest[DNWb; DNWt] + + return ρ / str(ρ) +end + +function reduced_densitymatrix1x2( + 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] + + @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] * + M_east[DEt; DEb] * + M_southeast[DSEb; DSEt] * + M_southwest[DSWb; DSWt] + + return ρ / str(ρ) +end 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..37e96457a --- /dev/null +++ b/src/algorithms/contractions/local_patch/expr_utils.jl @@ -0,0 +1,80 @@ +# Contraction expression utilities for local patches +# -------------------------------------------------- + +""" +$(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 56364cdde..000000000 --- a/src/algorithms/contractions/localoperator.jl +++ /dev/null @@ -1,621 +0,0 @@ -# Contraction of local operators on arbitrary lattice locations -# ------------------------------------------------------------- -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) - -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`. -# The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance). - -function _contract_corner_expr(rowrange, colrange) - 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])) - - 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( - 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 edges_N, edges_E, edges_S, edges_W -end - -function _contract_state_expr(rowrange, colrange, height, cartesian_inds = nothing) - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - - return map(1:height) 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 - 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,), - ( - 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, - ), - ) - end - end -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." - rmin, rmax = extrema(rowrange) - cmin, cmax = extrema(colrange) - gridsize = (rmax - rmin + 1, cmax - cmin + 1) - return map(1:height) 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 - 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)) - 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 - 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, - ), - ) - end - 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 -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 - ) - 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 - ) - 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 - ) 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) - - 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)..., - ) - - returnex = quote - @autoopt @tensor $multiplication_ex - end - return macroexpand(@__MODULE__, returnex) -end - -Base.@deprecate( - contract_local_norm(inds::NTuple, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv), - 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`. - -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 - -function reduced_densitymatrix( - inds::Vector{CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) - 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::CTMRGEnv - ) - 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 - ) - 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, ket, env) - -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv - ) where {N}, - reduced_densitymatrix(collect(inds), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, state::InfinitePEPO, env::CTMRGEnv - ) where {N}, - reduced_densitymatrix(collect(inds), state, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, CartesianIndex{2}}, ket::InfinitePEPO, bra::InfinitePEPO, env::CTMRGEnv - ) 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 - ) 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 - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), ket, bra, env) -) -Base.@deprecate( - reduced_densitymatrix( - inds::NTuple{N, Tuple{Int, Int}}, state::InfinitePEPO, env::CTMRGEnv - ) where {N}, - reduced_densitymatrix(collect(CartesianIndex.(inds)), state, env) -) - -# 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 = - 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) - - @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 - -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( - 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 = - edge(env, NORTH, row - 1, col) * - twistdual(corner(env, NORTHEAST, row - 1, col + 1), 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_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) - - @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 = - 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_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) - - @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 - -@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/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..324d2ebb5 --- /dev/null +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -0,0 +1,89 @@ +# Direct local patch contractions +# ------------------------------- + +""" + 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) + +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, 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..2767c2bc0 --- /dev/null +++ b/src/algorithms/expectation_value/reduced_densitymatrix.jl @@ -0,0 +1,107 @@ +# 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 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`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. +""" +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`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. +""" +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`. Falls back to +[`_contract_densitymatrix`](@ref) for environments without such a specialization. +""" +reduced_densitymatrix1x2(ind, ket, bra, env) = + _contract_densitymatrix((Val(ind), Val(ind + CartesianIndex(0, 1))), (ket, bra), env) 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 b53e0f443..533fac679 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -1,243 +1,3 @@ -""" - expectation_value(state, O::LocalOperator, env::CTMRGEnv) - expectation_value(bra, O::LocalOperator, ket, env::CTMRGEnv) - -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::CTMRGEnv - ) where {S <: InfiniteState} - checklattice(bra, O, ket) - term_vals = dtmap(collect(O.terms)) do (inds, operator) # OhMyThreads can't iterate over O.terms directly - ρ = reduced_densitymatrix(inds, ket, bra, env) - return trmul(operator, ρ) - 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 - ) - 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, ρ) - end - return sum(term_vals) -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::CTMRGEnv, - ) 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 - ) - return expectation_value(pf, CartesianIndex(op[1]) => op[2], env) -end - -""" - cost_function(peps::InfinitePEPS, env::CTMRGEnv, 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) - 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 -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) - -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 - 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] -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_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) - -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 -@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} @@ -365,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/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...) diff --git a/src/utility/tensor_traces.jl b/src/utility/tensor_traces.jl new file mode 100644 index 000000000..f99409180 --- /dev/null +++ b/src/utility/tensor_traces.jl @@ -0,0 +1,29 @@ +# 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 `str(H * ρ)` without forming `H * ρ`. + +See also [`str`](@ref). +""" +@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 dc7931ce7..603aec675 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -64,31 +64,6 @@ 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))) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index f31a994ce..018a68a82 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -235,6 +235,7 @@ using .TimeEvolTFIsingFiniteT module ToolboxDensityMatrices include("toolbox/densitymatrices.jl") export toolbox_single_layer_densitymatrix, toolbox_double_layer_densitymatrix + export toolbox_densitymatrix_too_many_layers, toolbox_densitymatrix_generic_fallback end using .ToolboxDensityMatrices diff --git a/test/testsuite/toolbox/densitymatrices.jl b/test/testsuite/toolbox/densitymatrices.jl index 63e7cca53..225240dcd 100644 --- a/test/testsuite/toolbox/densitymatrices.jl +++ b/test/testsuite/toolbox/densitymatrices.jl @@ -1,6 +1,6 @@ using TensorKit using PEPSKit -using PEPSKit: contract_local_operator, contract_local_norm +using PEPSKit: contract_local_operator, contract_local_norm, _contract_site using Test using TestExtras: @testinferred using Adapt @@ -23,10 +23,16 @@ function toolbox_single_layer_densitymatrix(AT) @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( @@ -36,6 +42,9 @@ function toolbox_single_layer_densitymatrix(AT) 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 end @@ -63,7 +72,12 @@ function toolbox_double_layer_densitymatrix(AT) @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 + # 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( @@ -78,6 +92,50 @@ function toolbox_double_layer_densitymatrix(AT) 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 +end + +function toolbox_densitymatrix_too_many_layers(AT) + return @testset "Three or more PEPO layers are rejected ($AT)" begin + d, D, χ = ds[Trivial], Ds[Trivial], χs[Trivial] + ρ = adapt(AT, 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 +end + +function toolbox_densitymatrix_generic_fallback(AT) + return @testset "Fixed-size fast paths agree with the generic fallback ($I) ($AT)" for I in keys(ds) + d, D, χ = ds[I], Ds[I], χs[I] + peps = adapt(AT, 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 end diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl index 8d0f3d08e..3030ec187 100644 --- a/test/toolbox/densitymatrices.jl +++ b/test/toolbox/densitymatrices.jl @@ -9,4 +9,6 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true" if !is_buildkite TestSuite.toolbox_single_layer_densitymatrix(Vector) TestSuite.toolbox_double_layer_densitymatrix(Vector) + TestSuite.toolbox_densitymatrix_too_many_layers(Vector) + TestSuite.toolbox_densitymatrix_generic_fallback(Vector) end