diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index ed4fade63..89b803502 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -1574,6 +1574,30 @@ function _pepo_pepotensor_expr( ) end +# PEPO layers slice h +function _pepo_pepotracetensor_expr( + tensorname, + h::Int, + H::Int, + args...; + contract_north=nothing, + contract_east=nothing, + contract_south=nothing, + contract_west=nothing, +) + layer = Symbol(h) + return tensorexpr( + tensorname, + (physicallabel(mod1(h + 1, H), args...), physicallabel(h, args...)), + ( + virtuallabel(_north_labels(layer, args...; contract=contract_north)...), + virtuallabel(_east_labels(layer, args...; contract=contract_east)...), + virtuallabel(_south_labels(layer, args...; contract=contract_south)...), + virtuallabel(_west_labels(layer, args...; contract=contract_west)...), + ), + ) +end + # PEPOSandwich function _pepo_sandwich_expr(sandwichname, H::Int, args...; kwargs...) ket_e = _pepo_pepstensor_expr(:(ket($sandwichname)), :top, 1, args...; kwargs...) @@ -1585,6 +1609,15 @@ function _pepo_sandwich_expr(sandwichname, H::Int, args...; kwargs...) return ket_e, bra_e, pepo_es end +# PEPOTraceSandwich +function _pepotrace_sandwich_expr(sandwichname, H::Int, args...; kwargs...) + pepo_es = map(1:H) do h + return _pepo_pepotracetensor_expr(:(pepo($sandwichname, $h)), h, H, args...; kwargs...) + end + + return pepo_es +end + ## Corner expressions function _corner_expr(cornername, codom_label, dom_label, args...) @@ -1608,6 +1641,14 @@ function _pepo_edge_expr(edgename, codom_label, dom_label, dir, H::Int, args...) ) end +function _pepotrace_edge_expr(edgename, codom_label, dom_label, dir, H::Int, args...) + return tensorexpr( + edgename, + (envlabel(codom_label, args...), ntuple(i -> virtuallabel(dir, i, args...), H)...), + (envlabel(dom_label, args...),), + ) +end + ## Enlarged corner (quadrant) expressions function _pepo_enlarged_corner_expr( @@ -1630,6 +1671,22 @@ function _pepo_enlarged_corner_expr( ) end +function _pepotrace_enlarged_corner_expr( + cornername, codom_label, dom_label, codom_dir, dom_dir, H::Int, args... +) + return tensorexpr( + cornername, + ( + envlabel(codom_label, args...), + ntuple(i -> virtuallabel(codom_dir, i, args...), H)..., + ), + ( + envlabel(dom_label, args...), + ntuple(i -> virtuallabel(dom_dir, i, args...), H)..., + ), + ) +end + ## Environment expressions function _pepo_env_expr( @@ -1689,6 +1746,19 @@ function _pepo_codomain_projector_expr( ) end +function _pepotrace_codomain_projector_expr( + projname, codom_label, dom_label, dom_dir, H::Int, args... +) + return tensorexpr( + projname, + (envlabel(codom_label, args...),), + ( + envlabel(dom_label, args...), + ntuple(i -> virtuallabel(dom_dir, i, args...), H)..., + ), + ) +end + function _pepo_domain_projector_expr( projname, codom_label, codom_dir, dom_label, H::Int, args... ) @@ -1704,6 +1774,19 @@ function _pepo_domain_projector_expr( ) end +function _pepotrace_domain_projector_expr( + projname, codom_label, codom_dir, dom_label, H::Int, args... +) + return tensorexpr( + projname, + ( + envlabel(codom_label, args...), + ntuple(i -> virtuallabel(codom_dir, i, args...), H)..., + ), + (envlabel(dom_label, args...),), + ) +end + # # PEPO Contractions # @@ -1753,6 +1836,48 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) end +@generated function _contract_site( + C_northwest, + C_northeast, + C_southeast, + C_southwest, + E_north::CTMRGEdgeTensor{T,S,N}, + E_east::CTMRGEdgeTensor{T,S,N}, + E_south::CTMRGEdgeTensor{T,S,N}, + E_west::CTMRGEdgeTensor{T,S,N}, + O::PEPOTraceSandwich{H}, +) where {T,S,N,H} + @assert N == H + 1 + + 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 = _pepotrace_edge_expr(:E_north, :NNW, :NNE, :N, H) + E_east_e = _pepotrace_edge_expr(:E_east, :ENE, :ESE, :E, H) + E_south_e = _pepotrace_edge_expr(:E_south, :SSE, :SSW, :S, H) + E_west_e = _pepotrace_edge_expr(:E_west, :WSW, :WNW, :W, H) + + pepo_es = _pepotrace_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, + pepo_es..., + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $rhs)) +end + ## Enlarged corner contractions @generated function enlarge_northwest_corner( @@ -1784,6 +1909,26 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function enlarge_northwest_corner( + E_west::CTMRGEdgeTensor{T,S,N}, + C_northwest::CTMRGCornerTensor, + E_north::CTMRGEdgeTensor{T,S,N}, + O::PEPOTraceSandwich{H}, +) where {T,S,N,H} + @assert N == H + 1 + + E_west_e = _pepotrace_edge_expr(:E_west, :SW, :WNW, :W, H) + C_northwest_e = _corner_expr(:C_northwest, :WNW, :NNW) + E_north_e = _pepotrace_edge_expr(:E_north, :NNW, :NE, :N, H) + pepo_es = _pepotrace_sandwich_expr(:O, H) + + C_out_e = _pepotrace_enlarged_corner_expr(:C_northwest´, :SW, :NE, :S, :E, H) + + rhs = Expr(:call, :*, E_west_e, C_northwest_e, E_north_e, pepo_es...) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function enlarge_northeast_corner( E_north::CTMRGEdgeTensor{T,S,N}, C_northeast::CTMRGCornerTensor, @@ -1813,6 +1958,26 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function enlarge_northeast_corner( + E_north::CTMRGEdgeTensor{T,S,N}, + C_northeast::CTMRGCornerTensor, + E_east::CTMRGEdgeTensor{T,S,N}, + O::PEPOTraceSandwich{H}, +) where {T,S,N,H} + @assert N == H + 1 + + E_north_e = _pepotrace_edge_expr(:E_north, :NW, :NNE, :N, H) + C_northeast = _corner_expr(:C_northeast, :NNE, :ENE) + E_east_e = _pepotrace_edge_expr(:E_east, :ENE, :SE, :E, H) + pepo_es = _pepotrace_sandwich_expr(:O, H) + + C_out_e = _pepotrace_enlarged_corner_expr(:C_northeast´, :NW, :SE, :W, :S, H) + + rhs = Expr(:call, :*, E_north_e, C_northeast, E_east_e, pepo_es...) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function enlarge_southeast_corner( E_east::CTMRGEdgeTensor{T,S,N}, C_southeast::CTMRGCornerTensor, @@ -1842,6 +2007,26 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function enlarge_southeast_corner( + E_east::CTMRGEdgeTensor{T,S,N}, + C_southeast::CTMRGCornerTensor, + E_south::CTMRGEdgeTensor{T,S,N}, + O::PEPOTraceSandwich{H}, +) where {T,S,N,H} + @assert N == H + 1 + + E_east_e = _pepotrace_edge_expr(:E_east, :NE, :ESE, :E, H) + C_southeast_e = _corner_expr(:C_southeast, :ESE, :SSE) + E_south_e = _pepotrace_edge_expr(:E_south, :SSE, :SW, :S, H) + pepo_es = _pepotrace_sandwich_expr(:O, H) + + C_out_e = _pepotrace_enlarged_corner_expr(:C_southeast´, :NE, :SW, :N, :W, H) + + rhs = Expr(:call, :*, E_east_e, C_southeast_e, E_south_e, pepo_es...) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function enlarge_southwest_corner( E_south::CTMRGEdgeTensor{T,S,N}, C_southwest::CTMRGCornerTensor, @@ -1871,6 +2056,26 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function enlarge_southwest_corner( + E_south::CTMRGEdgeTensor{T,S,N}, + C_southwest::CTMRGCornerTensor, + E_west::CTMRGEdgeTensor{T,S,N}, + O::PEPOTraceSandwich{H}, +) where {T,S,N,H} + @assert N == H + 1 + + E_south_e = _pepotrace_edge_expr(:E_south, :SE, :SSW, :S, H) + C_southwest_e = _corner_expr(:C_southwest, :SSW, :WSW) + E_west_e = _pepotrace_edge_expr(:E_west, :WSW, :NW, :W, H) + pepo_es = _pepotrace_sandwich_expr(:O, H) + + C_out_e = _pepotrace_enlarged_corner_expr(:C_southwest´, :SE, :NW, :E, :N, H) + + rhs = Expr(:call, :*, E_south_e, C_southwest_e, E_west_e, pepo_es...) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + ## Projector contractions: skip, since these are somehow never used ## HalfInfiniteEnvironment contractions @@ -2159,6 +2364,25 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function renormalize_northwest_corner( + E_west, C_northwest, E_north, P_left, P_right, A::PEPOTraceSandwich{H} +) where {H} + C_out_e = _corner_expr(:corner, :out, :in) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :S, :S, H) + E_west_e = _pepotrace_edge_expr(:E_west, :S, :WNW, :W, H) + C_northwest_e = _corner_expr(:C_northwest, :WNW, :NNW) + E_north_e = _pepotrace_edge_expr(:E_north, :NNW, :E, :N, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :E, :E, :in, H) + + rhs = Expr( + :call, :*, P_right_e, E_west_e, C_northwest_e, E_north_e, pepo_es..., P_left_e + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function renormalize_northeast_corner( E_north, C_northeast, E_east, P_left, P_right, A::PEPOSandwich{H} ) where {H} @@ -2187,6 +2411,25 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function renormalize_northeast_corner( + E_north, C_northeast, E_east, P_left, P_right, A::PEPOTraceSandwich{H} +) where {H} + C_out_e = _corner_expr(:corner, :out, :in) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :W, :W, H) + E_north_e = _pepotrace_edge_expr(:E_north, :W, :NNE, :N, H) + C_northeast_e = _corner_expr(:C_northeast, :NNE, :ENE) + E_east_e = _pepotrace_edge_expr(:E_east, :ENE, :S, :E, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :S, :S, :in, H) + + rhs = Expr( + :call, :*, P_right_e, E_north_e, C_northeast_e, E_east_e, pepo_es..., P_left_e + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function renormalize_southeast_corner( E_east, C_southeast, E_south, P_left, P_right, A::PEPOSandwich{H} ) where {H} @@ -2215,6 +2458,25 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function renormalize_southeast_corner( + E_east, C_southeast, E_south, P_left, P_right, A::PEPOTraceSandwich{H} +) where {H} + C_out_e = _corner_expr(:corner, :out, :in) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :N, :N, H) + E_east_e = _pepotrace_edge_expr(:E_east, :N, :EWE, :E, H) + C_southeast_e = _corner_expr(:C_southeast, :ESE, :SSE) + E_south_e = _pepotrace_edge_expr(:E_south, :SSE, :W, :S, H) + ket_e, bra_e, pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :W, :W, :in, H) + + rhs = Expr( + :call, :*, P_right_e, E_east_e, C_southeast_e, E_south_e, pepo_es..., P_left_e + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + @generated function renormalize_southwest_corner( E_south, C_southwest, E_west, P_left, P_right, A::PEPOSandwich{H} ) where {H} @@ -2243,6 +2505,25 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) end +@generated function renormalize_southwest_corner( + E_south, C_southwest, E_west, P_left, P_right, A::PEPOTraceSandwich{H} +) where {H} + C_out_e = _corner_expr(:corner, :out, :in) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :E, :E, H) + E_south_e = _pepotrace_edge_expr(:E_south, :E, :SSW, :S, H) + C_southwest_e = _corner_expr(:C_southwest, :SSW, :WSW) + E_west_e = _pepotrace_edge_expr(:E_west, :WSW, :N, :W, H) + ket_e, bra_e, pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :N, :N, :in, H) + + rhs = Expr( + :call, :*, P_right_e, E_south_e, C_southwest_e, E_west_e, pepo_es..., P_left_e + ) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $C_out_e := $rhs)) +end + ## Edge renormalization contractions @generated function renormalize_north_edge( @@ -2271,6 +2552,23 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) end +@generated function renormalize_north_edge( + E_north::CTMRGEdgeTensor{T,S,N}, P_left, P_right, A::PEPOTraceSandwich{H} +) where {T,S,N,H} + @assert N == H + 1 + + E_out_e = _pepotrace_edge_expr(:edge, :out, :in, :S, H) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :W, :W, H) + E_north_e = _pepotrace_edge_expr(:E_north, :W, :E, :N, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :E, :E, :in, H) + + rhs = Expr(:call, :*, P_right_e, E_north_e, pepo_es..., P_left_e) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) +end + @generated function renormalize_east_edge( E_east::CTMRGEdgeTensor{T,S,N}, P_bottom, P_top, A::PEPOSandwich{H} ) where {T,S,N,H} @@ -2297,6 +2595,23 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) end +@generated function renormalize_east_edge( + E_east::CTMRGEdgeTensor{T,S,N}, P_bottom, P_top, A::PEPOTraceSandwich{H} +) where {T,S,N,H} + @assert N == H + 1 + + E_out_e = _pepotrace_edge_expr(:edge, :out, :in, :W, H) + + P_top_e = _pepotrace_codomain_projector_expr(:P_top, :out, :N, :N, H) + E_east_e = _pepotrace_edge_expr(:E_east, :N, :S, :E, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_bottom_e = _pepotrace_domain_projector_expr(:P_bottom, :S, :S, :in, H) + + rhs = Expr(:call, :*, P_top_e, E_east_e, pepo_es..., P_bottom_e) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) +end + @generated function renormalize_south_edge( E_south::CTMRGEdgeTensor{T,S,N}, P_left, P_right, A::PEPOSandwich{H} ) where {T,S,N,H} @@ -2323,6 +2638,23 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) end +@generated function renormalize_south_edge( + E_south::CTMRGEdgeTensor{T,S,N}, P_left, P_right, A::PEPOTraceSandwich{H} +) where {T,S,N,H} + @assert N == H + 1 + + E_out_e = _pepotrace_edge_expr(:edge, :out, :in, :N, H) + + P_right_e = _pepotrace_codomain_projector_expr(:P_right, :out, :E, :E, H) + E_south_e = _pepotrace_edge_expr(:E_south, :E, :W, :S, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_left_e = _pepotrace_domain_projector_expr(:P_left, :W, :W, :in, H) + + rhs = Expr(:call, :*, P_right_e, E_south_e, pepo_es..., P_left_e) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) +end + @generated function renormalize_west_edge( E_west::CTMRGEdgeTensor{T,S,N}, P_bottom, P_top, A::PEPOSandwich{H} ) where {T,S,N,H} @@ -2348,3 +2680,20 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) end + +@generated function renormalize_west_edge( + E_west::CTMRGEdgeTensor{T,S,N}, P_bottom, P_top, A::PEPOTraceSandwich{H} +) where {T,S,N,H} + @assert N == H + 1 + + E_out_e = _pepotrace_edge_expr(:edge, :out, :in, :E, H) + + P_top_e = _pepotrace_codomain_projector_expr(:P_top, :out, :S, :S, H) + E_west_e = _pepotrace_edge_expr(:E_west, :S, :N, :W, H) + pepo_es = _pepotrace_sandwich_expr(:A, H) + P_bottom_e = _pepotrace_domain_projector_expr(:P_bottom, :N, :N, :in, H) + + rhs = Expr(:call, :*, P_top_e, E_west_e, pepo_es..., P_bottom_e) + + return macroexpand(@__MODULE__, :(return @autoopt @tensor $E_out_e := $rhs)) +end diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 1e9324319..ae8a29ed0 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -384,6 +384,27 @@ function contract_local_tensor( ) end +function contract_local_tensor( + ind::Tuple{Int,Int,Int}, + O::PEPOTensor, + network::InfiniteSquareNetwork{<:PEPOTraceSandwich}, + env::CTMRGEnv, +) + r, c, h = ind + sandwich´ = Base.setindex(network[r, c], O, h) + return _contract_site( + env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], + env.corners[NORTHEAST, _prev(r, end), _next(c, end)], + env.corners[SOUTHEAST, _next(r, end), _next(c, end)], + env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], + env.edges[NORTH, _prev(r, end), c], + env.edges[EAST, r, _next(c, end)], + env.edges[SOUTH, _next(r, end), c], + env.edges[WEST, r, _prev(c, end)], + 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/networks/local_sandwich.jl b/src/networks/local_sandwich.jl index 431f67f48..d5711c22f 100644 --- a/src/networks/local_sandwich.jl +++ b/src/networks/local_sandwich.jl @@ -55,39 +55,6 @@ function virtualspace(O::PEPSSandwich, dir) return virtualspace(ket(O), dir) ⊗ virtualspace(bra(O), dir)' end -# not overloading MPOTensor because that defines AbstractTensorMap{<:Any,S,2,2}(::PEPSTensor, ::PEPSTensor) -# ie type piracy -mpotensor(top::PEPSTensor) = mpotensor((top, top)) -function mpotensor((top, bot)::PEPSSandwich) - @assert virtualspace(top, NORTH) == dual(virtualspace(top, SOUTH)) && - virtualspace(bot, NORTH) == dual(virtualspace(bot, SOUTH)) && - virtualspace(top, EAST) == dual(virtualspace(top, WEST)) && - virtualspace(bot, EAST) == dual(virtualspace(bot, WEST)) && - isdual(virtualspace(top, NORTH)) && - isdual(virtualspace(bot, NORTH)) && - isdual(virtualspace(top, EAST)) && - isdual(virtualspace(bot, EAST)) "Method not yet implemented for given virtual spaces" - - F_west = isomorphism( - storagetype(top), - fuse(virtualspace(top, WEST), virtualspace(bot, WEST)'), - virtualspace(top, WEST) ⊗ virtualspace(bot, WEST)', - ) - F_south = isomorphism( - storagetype(top), - fuse(virtualspace(top, SOUTH), virtualspace(bot, SOUTH)'), - virtualspace(top, SOUTH) ⊗ virtualspace(bot, SOUTH)', - ) - @tensor O[west south; north east] := - top[phys; top_north top_east top_south top_west] * - conj(bot[phys; bot_north bot_east bot_south bot_west]) * - twist(F_west, 3)[west; top_west bot_west] * - twist(F_south, 3)[south; top_south bot_south] * - conj(F_west[east; top_east bot_east]) * - conj(F_south[north; top_north bot_north]) - return O -end - ## PEPO const PEPOSandwich{N,T<:PEPSTensor,P<:PEPOTensor} = Tuple{T,T,Vararg{P,N}} @@ -104,3 +71,154 @@ function virtualspace(O::PEPOSandwich, dir) virtualspace(bra(O), dir)', ]) end + +## PEPOTraceSandwich +# In a CTMRG contraction, the top physical leg of the top PEPOTensor is contracted with the bottom physical leg of the bottom PEPOTensor + +const PEPOTraceSandwich{N,P<:PEPOTensor} = Tuple{Vararg{P,N}} + +pepo(O::PEPOTraceSandwich) = O[1:end] +pepo(O::PEPOTraceSandwich, i::Int) = O[i] + +function virtualspace(O::PEPOTraceSandwich, dir) + return prod([virtualspace.(pepo(O), Ref(dir))...]) +end + +# +# Expressions for mpotensor +# + +function _pepo_fuser_tensor_expr(tensorname, H::Int, dir, args...;) + return tensorexpr( + tensorname, + (virtuallabel(dir, :fuser, args...),), + ( + virtuallabel(dir, :top), + ntuple(h -> virtuallabel(_virtual_labels(dir, :mid, h, args...)...), H)..., + virtuallabel(dir, :bot), + ), + ) +end + +function _pepotrace_fuser_tensor_expr(tensorname, H::Int, dir, args...;) + return tensorexpr( + tensorname, + (virtuallabel(dir, :fuser, args...),), + (ntuple(h -> virtuallabel(_virtual_labels(dir, h, args...)...), H)...,), + ) +end + +@generated function _mpotensor_contraction( + F_north::AbstractTensorMap{T,S}, + F_east::AbstractTensorMap{T,S}, + F_south::AbstractTensorMap{T,S}, + F_west::AbstractTensorMap{T,S}, + O::PEPOSandwich{H}, +) where {T,S,H} + fuser_north = _pepo_fuser_tensor_expr(:F_north, H, :N) + fuser_east = _pepo_fuser_tensor_expr(:F_east, H, :E) + fuser_south = _pepo_fuser_tensor_expr(:F_south, H, :S) + fuser_west = _pepo_fuser_tensor_expr(:F_west, H, :W) + ket_e, bra_e, pepo_es = _pepo_sandwich_expr(:O, H) + + result = tensorexpr(:res, (:D_W_fuser, :D_S_fuser), (:D_N_fuser, :D_E_fuser)) + + rhs = Expr( + :call, + :*, + Expr(:call, :conj, fuser_north), + Expr(:call, :conj, fuser_east), + fuser_south, + fuser_west, + ket_e, + Expr(:call, :conj, bra_e), + pepo_es..., + ) + return macroexpand(@__MODULE__, :(return @autoopt @tensor $result := $rhs)) +end + +@generated function _mpotensor_contraction( + F_north::AbstractTensorMap{T,S}, + F_east::AbstractTensorMap{T,S}, + F_south::AbstractTensorMap{T,S}, + F_west::AbstractTensorMap{T,S}, + O::PEPOTraceSandwich{H}, +) where {T,S,H} + fuser_north = _pepotrace_fuser_tensor_expr(:F_north, H, :N) + fuser_east = _pepotrace_fuser_tensor_expr(:F_east, H, :E) + fuser_south = _pepotrace_fuser_tensor_expr(:F_south, H, :S) + fuser_west = _pepotrace_fuser_tensor_expr(:F_west, H, :W) + pepo_es = _pepotrace_sandwich_expr(:O, H) + + result = tensorexpr(:res, (:D_W_fuser, :D_S_fuser), (:D_N_fuser, :D_E_fuser)) + + rhs = Expr( + :call, + :*, + Expr(:call, :conj, fuser_north), + Expr(:call, :conj, fuser_east), + fuser_south, + fuser_west, + pepo_es..., + ) + return macroexpand(@__MODULE__, :(return @autoopt @tensor $result := $rhs)) +end + +# not overloading MPOTensor because that defines AbstractTensorMap{<:Any,S,2,2}(::PEPSTensor, ::PEPSTensor) +# ie type piracy +mpotensor(top::PEPSTensor) = mpotensor((top, top)) +function mpotensor(network::PEPOSandwich{H}) where {H} + @assert virtualspace(ket(network), NORTH) == dual(virtualspace(ket(network), SOUTH)) && + virtualspace(bra(network), NORTH) == dual(virtualspace(bra(network), SOUTH)) && + virtualspace(ket(network), EAST) == dual(virtualspace(ket(network), WEST)) && + virtualspace(bra(network), EAST) == dual(virtualspace(bra(network), WEST)) "Method not yet implemented for given virtual spaces" + F_west = isomorphism( + storagetype(network[1]), + fuse(virtualspace(network, WEST)), + virtualspace(network, WEST), + ) + F_south = isomorphism( + storagetype(network[1]), + fuse(virtualspace(network, SOUTH)), + virtualspace(network, SOUTH), + ) + F_east = copy(F_west) + F_north = copy(F_south) + + isdual(virtualspace(ket(network), NORTH)) || twist!(F_north, 2) + isdual(virtualspace(bra(network), NORTH)) || twist!(F_north, H + 3) + isdual(virtualspace(ket(network), EAST)) || twist!(F_east, 2) + isdual(virtualspace(bra(network), EAST)) || twist!(F_east, H + 3) + for h in 1:H + @assert virtualspace(network[h], NORTH) == dual(virtualspace(network[h], SOUTH)) && + virtualspace(network[h], EAST) == dual(virtualspace(network[h], WEST)) "Method not yet implemented for given virtual spaces" + isdual(virtualspace(network[h], NORTH)) || twist!(F_north, h + 2) + isdual(virtualspace(network[h], EAST)) || twist!(F_east, h + 2) + end + twist!(F_west, H + 3) + twist!(F_south, H + 3) + return _mpotensor_contraction(F_north, F_east, F_south, F_west, network) +end + +function mpotensor(network::PEPOTraceSandwich{H}) where {H} + F_west = isomorphism( + storagetype(network[1]), + fuse(virtualspace(network, WEST)), + virtualspace(network, WEST), + ) + F_south = isomorphism( + storagetype(network[1]), + fuse(virtualspace(network, SOUTH)), + virtualspace(network, SOUTH), + ) + F_east = copy(F_west) + F_north = copy(F_south) + + for h in 1:H + @assert virtualspace(network[h], NORTH) == dual(virtualspace(network[h], SOUTH)) && + virtualspace(network[h], EAST) == dual(virtualspace(network[h], WEST)) "Method not yet implemented for given virtual spaces" + isdual(virtualspace(network[h], NORTH)) || twist!(F_north, h) + isdual(virtualspace(network[h], EAST)) || twist!(F_east, h) + end + return _mpotensor_contraction(F_north, F_east, F_south, F_west, network) +end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 054d4cfdf..048b21c31 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -175,6 +175,33 @@ function InfiniteSquareNetwork(top::InfinitePEPS, mid::InfinitePEPO, bot::Infini ) end +function InfiniteSquareNetwork(mid::InfinitePEPO) + return InfiniteSquareNetwork(map(tuple, eachslice(unitcell(mid); dims=3)...)) +end + +""" + _dag(O::PEPOTensor) + +Calculate the conjugate of an operator O, while permuting the physical indices. +Twists are included to ensure the correct result for the adjoint function. +""" +function _dag(O::PEPOTensor) + @tensor O_conj[-1 -2; -3 -4 -5 -6] := conj(O[-2 -1; -3 -4 -5 -6]) + isdual(codomain(O_conj)[1]) && twist!(O_conj, 1) + isdual(codomain(O_conj)[2]) || twist!(O_conj, 2) + return O_conj +end + +""" + adjoint(O::InfinitePEPO) + +Create the adjoint of an InfinitePEPO. +This is defined such that `dot(psi, O, phi) == dot(psi, O * phi) == dot(O' * psi, phi)` for any two states `psi` and `phi`. +""" +function Base.adjoint(O::InfinitePEPO) + return InfinitePEPO(_dag.(unitcell(O))) +end + ## Vector interface function VectorInterface.scalartype(::Type{NT}) where {NT<:InfinitePEPO} diff --git a/test/ctmrg/pepotrace.jl b/test/ctmrg/pepotrace.jl new file mode 100644 index 000000000..abddf90b7 --- /dev/null +++ b/test/ctmrg/pepotrace.jl @@ -0,0 +1,190 @@ +using Test +using Random +using LinearAlgebra +using TensorKit +using PEPSKit + +function fuse_PEPO_PEPS(O, A, fuser) + @tensor t[-1; -3 -4 -5 -6] := + O[-1 1; 2 4 6 8] * + A[1; 3 5 7 9] * + fuser[2 3; -3] * + fuser[4 5; -4] * + conj(fuser[6 7; -5]) * + conj(fuser[8 9; -6]) + return t +end + +function fuse_PEPO_PEPO(O₁, O₂, fuser) + @tensor t[-1 -2; -3 -4 -5 -6] := + O₁[-1 1; 2 4 6 8] * + O₂[1 -2; 3 5 7 9] * + fuser[2 3; -3] * + fuser[4 5; -4] * + conj(fuser[6 7; -5]) * + conj(fuser[8 9; -6]) + return t +end + +χenv = 24 +T = ComplexF64 + +I = fℤ₂ +pspace = Vect[I](0 => 1, 1 => 1) +vspace = Vect[I](0 => 1, 1 => 1) +envspace = Vect[I](0 => χenv, 1 => χenv) + +Random.seed!(97646475) + +# Construct random PEPO tensors +O = randn(T, pspace ⊗ pspace', vspace ⊗ vspace ⊗ vspace' ⊗ vspace') +O = O + twist(flip(PEPSKit._dag(O), 3:6), [4 6]) +M = randn(T, pspace ⊗ pspace', vspace ⊗ vspace ⊗ vspace' ⊗ vspace') +M = M + twist(flip(PEPSKit._dag(M), 3:6), [4 6]) + +# Fuse a layer consisting of O-O and MO together +fuser = isomorphism(T, vspace ⊗ vspace, fuse(vspace, vspace)) +fuser_adj = isomorphism(T, vspace' ⊗ vspace, fuse(vspace', vspace)) + +O2 = fuse_PEPO_PEPO(O, O, fuser) +MO = fuse_PEPO_PEPO(M, O, fuser) + +# Create `InfiniteSquareNetwork`s of of both options +O_stack = fill(O, 1, 1, 2) +network = InfiniteSquareNetwork(InfinitePEPO(O_stack)); +network_fused = InfiniteSquareNetwork(InfinitePEPO(O2)); + +# Construct random PEPS tensors +ψ = randn(T, pspace, vspace ⊗ vspace ⊗ vspace' ⊗ vspace') +ψ = ψ / norm(ψ) +ϕ = randn(T, pspace, vspace ⊗ vspace ⊗ vspace' ⊗ vspace') +ϕ = ϕ / norm(ϕ) + +# cover all different flavors +ctm_styles = [:sequential, :simultaneous] +projector_algs = [:halfinfinite, :fullinfinite] + +@testset "mpotensor for PEPOTraceSandwich using $alg with $projector_alg" for ( + alg, projector_alg +) in Iterators.product( + ctm_styles, projector_algs +) + # Test whether the mpotensor of a double layer of PEPOs is the same as in the case where the layers are fused + mpo = InfiniteSquareNetwork(map(PEPSKit.mpotensor, PEPSKit.unitcell(network))) + mpo_fused = InfiniteSquareNetwork( + map(PEPSKit.mpotensor, PEPSKit.unitcell(network_fused)) + ) + env, = leading_boundary(CTMRGEnv(mpo, envspace), mpo; alg, maxiter=250, projector_alg) + env_fused, = leading_boundary( + CTMRGEnv(mpo_fused, envspace), mpo_fused; alg, maxiter=250, projector_alg + ) + + @test network_value(mpo, env) ≈ network_value(mpo_fused, env_fused) +end + +@testset "PEPO layers CTMRG contraction using $alg with $projector_alg" for ( + alg, projector_alg +) in Iterators.product( + ctm_styles, projector_algs +) + # Test whether calculating the environment of a double layer of PEPOs results in the same expectation value as in the case where the layers are fused + env, = leading_boundary( + CTMRGEnv(network, envspace), network; alg, maxiter=250, projector_alg + ) + env_fused, = leading_boundary( + CTMRGEnv(network_fused, envspace), network_fused; alg, maxiter=250, projector_alg + ) + + m = PEPSKit.contract_local_tensor((1, 1, 1), M, network, env) + m_fused = PEPSKit.contract_local_tensor((1, 1, 1), MO, network_fused, env_fused) + + nrm = PEPSKit._contract_site((1, 1), network, env) + nrm_fused = PEPSKit._contract_site((1, 1), network_fused, env_fused) + @test (m / nrm) ≈ (m_fused / nrm_fused) atol = 1e-9 +end + +projector_alg = projector_algs[1] # only use :halfinfinite for this test due to convergence issues +@testset "Test adjoint of an InfinitePEPO using $alg with $projector_alg" for alg in + ctm_styles + # Test the definition of the adjoint of an operator, i.e. dot(psi, O', phi) == dot(psi, O' * phi) == dot(O * psi, phi) for any two states psi and phi. + Oψ = fuse_PEPO_PEPS(O, ψ, fuser) + Odagϕ = fuse_PEPO_PEPS(twist(PEPSKit._dag(O), 1:4), ϕ, fuser_adj) + + network_ψOϕ = InfiniteSquareNetwork( + InfinitePEPS(ψ), InfinitePEPO(PEPSKit._dag(O)), InfinitePEPS(ϕ) + ) + network_ψ_Odagϕ = InfiniteSquareNetwork(InfinitePEPS(ψ), InfinitePEPS(Odagϕ)) + network_Oψ_ϕ = InfiniteSquareNetwork(InfinitePEPS(Oψ), InfinitePEPS(ϕ)) + + env_ψOϕ, = leading_boundary( + CTMRGEnv(network_ψOϕ, envspace), + network_ψOϕ; + alg, + maxiter=500, + projector_alg, + verbosity=2, + ) + env_ψ_Odagϕ, = leading_boundary( + CTMRGEnv(network_ψ_Odagϕ, envspace), + network_ψ_Odagϕ; + alg, + maxiter=500, + projector_alg, + verbosity=2, + ) + env_Oψ_ϕ, = leading_boundary( + CTMRGEnv(network_Oψ_ϕ, envspace), + network_Oψ_ϕ; + alg, + maxiter=500, + projector_alg, + verbosity=2, + ) + + overlap1 = network_value(network_ψOϕ, env_ψOϕ) + overlap2 = network_value(network_ψ_Odagϕ, env_ψ_Odagϕ) + overlap3 = network_value(network_Oψ_ϕ, env_Oψ_ϕ) + + @test overlap1 ≈ overlap2 atol = 5e-3 + @test overlap1 ≈ overlap3 atol = 5e-3 +end + +@testset "PEPO layers CTMRG contraction with its adjoint using $alg with $projector_alg" for ( + alg, projector_alg +) in Iterators.product( + ctm_styles, projector_algs +) + # Test whether calculating the environment of `PEPOSandwich` results in the same expectation value as when the PEPO is fused with the PEPS + + (Nr, Nc) = (1, 1) + Oinf = InfinitePEPO(O; unitcell=(Nr, Nc, 1)) + O_stack = fill(O, Nr, Nc, 2) + O_stack[:, :, 2] .= PEPSKit.unitcell(adjoint(Oinf)) + OOdag = InfinitePEPO(O_stack) + + Oψ = fuse_PEPO_PEPS(O, ψ, fuser) + Mψ = fuse_PEPO_PEPS(M, ψ, fuser) + + network_O = InfiniteSquareNetwork(InfinitePEPS(ψ), OOdag, InfinitePEPS(ψ)) + network_fused_MO = InfiniteSquareNetwork(InfinitePEPS(Mψ), InfinitePEPS(Oψ)) + network_fused_OO = InfinitePEPS(Oψ) + + env, = leading_boundary( + CTMRGEnv(network_O, envspace), network_O; alg, maxiter=400, projector_alg + ) + env_fused, = leading_boundary( + CTMRGEnv(network_fused_OO, envspace), + network_fused_OO; + alg, + maxiter=400, + projector_alg, + ) + + m = PEPSKit.contract_local_tensor((1, 1, 1), M, network_O, env) + m_fused = network_value(network_fused_MO, env_fused) + + nrm = PEPSKit._contract_site((1, 1), network_O, env) + nrm_fused = network_value(network_fused_OO, env_fused) + + @test (m / nrm) ≈ (m_fused / nrm_fused) atol = 5e-5 +end diff --git a/test/runtests.jl b/test/runtests.jl index 2ea46b4fc..47ffed2be 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,6 +32,9 @@ end @time @safetestset "PEPO" begin include("ctmrg/pepo.jl") end + @time @safetestset "PEPOTrace" begin + include("ctmrg/pepotrace.jl") + end @time @safetestset "correlation length" begin include("ctmrg/correlation_length.jl") end