From 257d3e2d9a49b0015c15a9a541dcaeb10496ced7 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 27 Feb 2025 10:12:16 +0800 Subject: [PATCH 1/2] Add NTU bond environments --- src/PEPSKit.jl | 2 + .../contractions/bondenv/benv_ntu.jl | 150 +++++++++++ .../contractions/bondenv/benv_tools.jl | 240 ++++++++++++++++++ src/operators/infinitepepo.jl | 2 + test/bondenv/benv_ntu.jl | 38 +++ test/runtests.jl | 3 + 6 files changed, 435 insertions(+) create mode 100644 src/algorithms/contractions/bondenv/benv_ntu.jl create mode 100644 test/bondenv/benv_ntu.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 0b7a817a2..423d3ed16 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -58,6 +58,7 @@ include("algorithms/contractions/vumps_contractions.jl") include("algorithms/contractions/bondenv/benv_tools.jl") include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") +include("algorithms/contractions/bondenv/benv_ntu.jl") include("algorithms/contractions/bondenv/benv_ctm.jl") include("algorithms/contractions/correlator/peps.jl") include("algorithms/contractions/correlator/pepo_1layer.jl") @@ -104,6 +105,7 @@ export fixedpoint export absorb_weight export ALSTruncation, FullEnvTruncation +export NNEnv, NNNEnv export su_iter, su3site_iter, simpleupdate, SimpleUpdate export InfiniteSquareNetwork diff --git a/src/algorithms/contractions/bondenv/benv_ntu.jl b/src/algorithms/contractions/bondenv/benv_ntu.jl new file mode 100644 index 000000000..df8f1a7c3 --- /dev/null +++ b/src/algorithms/contractions/bondenv/benv_ntu.jl @@ -0,0 +1,150 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + +""" +Algorithms to construct bond environment for Neighborhood Tensor Update (NTU). +""" +abstract type NeighbourEnv end + +""" +Construct the "NTU-NN" bond environment. +``` + (-1 +0)══(-1 +1) + ║ ║ + (+0 -1)═════X══ ═══Y═══(+0 +2) + ║ ║ + (+1 +0)══(+1 +1) +``` +""" +struct NNEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN" approximation. +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNEnv + ) where {T <: StateTensor, S <: InfinitePEPS} + neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)] + m = collect_neighbors(state, row, col, neighbors) + #= contraction indices + + (-1 +0) ══ Dn ══ (-1 +1) + ║ ║ + ........Dnw...... Dne + ║ : ║ + (+0 -1) ═══ X ══ Dw : De ══ Y ═══ (+0 +2) + ║ : ║ + Dsw :.......Dse........ + ║ ║ + (+1 +0) ══ Ds ══ (+1 +1) + =# + # bottom-left half + @autoopt @tensor benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] := + cor_se(m[1, 1])[Dse1 Dse0 Ds1 Ds0] * + cor_sw(m[1, 0])[Dsw1 Dsw0 Ds1 Ds0] * + edge_w(X, hair_w(m[0, -1]))[Dnw1 Dnw0 Dw1 Dw0 Dsw1 Dsw0] + normalize!(benv_sw, Inf) + # top-right half + @autoopt @tensor benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] := + cor_nw(m[-1, 0])[Dn1 Dn0 Dnw1 Dnw0] * + cor_ne(m[-1, 1])[Dne1 Dne0 Dn1 Dn0] * + edge_e(Y, hair_e(m[0, 2]))[Dne1 Dne0 Dse1 Dse0 De1 De0] + normalize!(benv_ne, Inf) + @tensor benv[Dw1 De1; Dw0 De0] := + benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] * benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] + normalize!(benv, Inf) + return benv +end + +""" +Construct the "NTU-NN+" bond environment. +``` + (-2 +0)┈┈(-2 +1) + ║ ║ + (-1 -1)┈(-1 +0)══(-1 +1)┈(-1 +2) + ┊ ║ ║ ┊ + (+0 -2)=(+0 -1)═════X══ ═══Y═══(+0 +2)=(+0 +3) + ┊ ║ ║ ┊ + (+1 -1)┈(+1 +0)══(+1 +1)┈(+1 +2) + ║ ║ + (+2 +0)┈┈(+2 +1) +``` +The tensors connected with dotted lines (e.g. (-1 +2)) +are splitted into two hair tensors. +""" +struct NNpEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN+" approximation. +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNpEnv + ) where {T <: StateTensor, S <: InfinitePEPS} + neighbors = [ + (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (1, 2), (0, 2), (-1, 2), + (-1, 1), (-1, 0), (0, -2), (2, 0), (2, 1), (0, 3), (-2, 1), (-2, 0), + ] + m = collect_neighbors(state, row, col, neighbors) + error("To be implemented.") +end + +""" +Construct the "NTU-NNN" bond environment. +``` + (-1 -1)=(-1 +0)══(-1 +1)=(-1 +2) + ║ ║ ║ ║ + (+0 -1)═════X══ ═══Y═══(+0 +2) + ║ ║ ║ ║ + (+1 -1)=(+1 +0)══(+1 +1)=(+1 +2) +``` +""" +struct NNNEnv <: NeighbourEnv end +""" +Calculates the bond environment within "NTU-NNN" approximation. +""" +function bondenv_ntu( + row::Int, col::Int, X::T, Y::T, state::S, alg::NNNEnv + ) where {T <: StateTensor, S <: InfinitePEPS} + neighbors = [ + (-1, -1), (0, -1), (1, -1), + (1, 0), (1, 1), (1, 2), (0, 2), + (-1, 2), (-1, 1), (-1, 0), + ] + m = collect_neighbors(state, row, col, neighbors) + #= left half + (-1 -1)══════(-1 +0)═ -1/-2 + ║ ║ + (+0 -1)════════ X ═══ -3/-4 + ║ ║ + ....D1..........D2......... + ║ ║ + (+1 -1)═ D3 ═(+1 +0)═ -5/-6 + =# + vecl = enlarge_corner_nw(cor_nw(m[-1, -1]), edge_n(m[-1, 0]), edge_w(m[0, -1]), X) + @tensor vecl[:] := + cor_sw(m[1, -1])[D11 D10 D31 D30] * + edge_s(m[1, 0])[D21 D20 -5 -6 D31 D30] * + vecl[D11 D10 D21 D20 -1 -2 -3 -4] + normalize!(vecl, Inf) + #= right half + -1/-2 ══ (-1 +1)═ D1 ═(-1 +2) + ║ ║ + ............D2..........D3... + ║ ║ + -3/-4 ═════ Y ═══════(+0 +2) + ║ ║ + -5/-6 ══ (+1 +1)═════(+1 +2) + =# + vecr = enlarge_corner_se(cor_se(m[1, 2]), edge_s(m[1, 1]), edge_e(m[0, 2]), Y) + @tensor vecr[:] := + edge_n(m[-1, 1])[D11 D10 D21 D20 -1 -2] * + cor_ne(m[-1, 2])[D31 D30 D11 D10] * + vecr[D21 D20 D31 D30 -3 -4 -5 -6] + normalize!(vecr, Inf) + # combine left and right part + @tensor benv[-1 -2; -3 -4] := vecl[1 2 -1 -3 3 4] * vecr[1 2 -2 -4 3 4] + normalize!(benv, Inf) + return benv +end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 91ad5958e..1c30740c2 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,3 +1,243 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace} +const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} +# Orthogonal tensors obtained PEPSTensor/PEPOTensor +# with one physical leg being factored out by `_qr_bond` const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace} const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace} +const StateTensor = Union{PEPSTensor, PEPOTensor, PEPSOrth, PEPOOrth} + +""" +Extract tensors in an InfinitePEPS or 1-layer InfinitePEPO +at positions `neighbors` relative to `(row, col)` +""" +function collect_neighbors( + state::InfiniteState, row::Int, col::Int, neighbors::Vector{Tuple{Int, Int}} + ) + Nr, Nc = size(state) + return Dict( + nb => state.A[mod1(row + nb[1], Nr), mod1(col + nb[2], Nc)] + for nb in neighbors + ) +end + +""" + benv_tensor(ket::T, bra::T, open_axs::Vector{Int}) where + {T <: Union{PEPSTensor, PEPOTensor, PEPSOrth, PEPOOrth}} + benv_tensor(ket::T, bra::T, open_axs::Vector{Int}, hairs::Vector{H}) where + {T <: Union{PEPSTensor, PEPOTensor, PEPSOrth, PEPOOrth}, H <: Union{Nothing, Hair}} + +Contract the physical axes and the virtual axes of `ket` with `bra` to obtain the tensor on the boundary of the bond environment. +Virtual axes specified by `open_axs` (in ascending order) are not contracted. + +# Examples + +- West "hair" tensor when `ket`, `bra` are `PEPSTensor` or `PEPOTensor` + (`open_axs = [2]`) + ``` + : + ╱| | ╱| + ┌-----ket----- 2 ┌-----ket----- 2 + | ╱ | | | ╱ | | + | | | | | | | | + | | | ╱ | | | ╱ + └---|-bra----- 1 └---|-bra----- 1 + |╱ |╱ | + : + ``` + For `PEPOTensor`, the remaining physical indices are traced out. + +- Northwest corner tensor when `ket`, `bra` are `PEPSOrth` or `PEPOOrth` + (`open_axs = [2, 3]`, `hairs = [h, nothing]`) + ``` + : + ╱| | ╱| + ┌-----ket----- 2 ┌-----ket----- 2 + | ╱ h | ╱ h + | 4 | | 4 | + | | | | + | ╱ | ╱ + └-----bra----- 1 └-----bra----- 1 + ╱ ╱ | + 3 3 : + ``` + +- West edge tensor when `ket`, `bra` are `PEPSTensor` or `PEPOTensor` + (`open_axs = [1, 2, 3]`) + ``` + 2 : 2 + ╱ | ╱ + ┌-----ket----- 4 ┌-----ket----- 4 + | ╱ | | ╱ | + | 6 | 1 | 6 | 1 + | | ╱ | | ╱ + └-----bra----- 3 └-----bra----- 3 + ╱ ╱ | + 5 5 : + ``` +""" +function benv_tensor( + ket::T, bra::T, open_axs::Vector{Int} + ) where {T <: StateTensor} + # no hair tensors to be attached to virtual legs + return benv_tensor(ket, bra, open_axs, fill(nothing, 4 - length(open_axs))) +end +function benv_tensor( + ket::T, bra::T, open_axs::Vector{Int}, hairs::Vector{H} + ) where {T <: StateTensor, H <: Union{Nothing, Hair}} + @assert length(hairs) == 4 - length(open_axs) + ket2 = if T <: PEPOTensor + first(fuse_physicalspaces(ket)) + elseif T <: PEPSOrth + deepcopy(ket) + else # PEPSTensor, PEPOOrth + twistdual(ket, 1) + end + nax = numind(ket2) + axs, open_axs2 = if T <: PEPSOrth + (1:4), open_axs + else # PEPSTensor, PEPOTensor (with physical legs fused), PEPOOrth + (2:5), open_axs .+ 1 + end + # contract with hair tensors + hair_axs = Tuple(ax for ax in axs if ax ∉ open_axs2) + for (h, ax) in zip(hairs, hair_axs) + if h === nothing + twistdual!(ket2, ax) + continue + end + # ensure the hair doesn't change the virtual spaces + @assert space(h, 1) == space(h, 2)' + ket_indices = collect(-1:-1:-nax) + ket_indices[ax] = 1 + ket2 = ncon([h, ket2], [[-ax, 1], ket_indices]) + end + # combine bra and ket + indexlist = [-collect(1:2:(2 * nax)), -collect(2:2:(2 * nax))] + for ax in 1:nax + if ax ∉ open_axs2 + indexlist[1][ax] = indexlist[2][ax] = ax + end + end + bra2 = (T <: PEPOTensor) ? first(fuse_physicalspaces(bra)) : bra + return ncon([bra2, ket2], indexlist, [true, false]) +end + +#= Free axes of different boundary tensors +(C/E/H mean corner/edge/hair) + + H_n + | + + C_nw - - E_n - - C_ne + | | | + + | | | + H_w - E_w - - ket - - E_e - H_e + | | | + + | | | + C_sw - - E_s - - C_se + + | + H_s +=# +const open_axs_hair = Dict(:n => [SOUTH], :e => [WEST], :s => [NORTH], :w => [EAST]) +const open_axs_cor = Dict( + :nw => [EAST, SOUTH], :ne => [SOUTH, WEST], :se => [NORTH, WEST], :sw => [NORTH, EAST] +) +const open_axs_edge = Dict( + :n => [EAST, SOUTH, WEST], + :e => [NORTH, SOUTH, WEST], + :s => [NORTH, EAST, WEST], + :w => [NORTH, EAST, SOUTH], +) + +# construction of hairs +for (dir, open_axs) in open_axs_hair + fname = Symbol("hair_", dir) + @eval begin + $(fname)(ket::StateTensor) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket::StateTensor, h1, h2, h3) = benv_tensor(ket, ket, $open_axs, [h1, h2, h3]) + end +end + +# construction of corners +for (dir, open_axs) in open_axs_cor + fname = Symbol("cor_", dir) + @eval begin + $(fname)(ket::StateTensor) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket::StateTensor, h1, h2) = benv_tensor(ket, ket, $open_axs, [h1, h2]) + end +end + +# construction of edges +for (dir, open_axs) in open_axs_edge + fname = Symbol("edge_", dir) + @eval begin + $(fname)(ket::StateTensor) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket::StateTensor, h) = benv_tensor(ket, ket, $open_axs, [h]) + end +end + +""" +Enlarge the northwest corner +``` + ctl══ D1 ══ et ══ -5/-6 + ║ ║ + D2 D3 + ║ ║ + el ══ D4 ══ X ═══ -7/-8 + ║ ║ + -1/-2 -3/-4 +``` +""" +function enlarge_corner_nw( + ctl::AbstractTensor{E, S, 4}, + et::AbstractTensor{E, S, 6}, el::AbstractTensor{E, S, 6}, + ket::T, bra::T = ket + ) where {E, S, T <: StateTensor} + if T <: PEPSOrth + return @tensoropt ctl2[:] := ctl[D11 D10 D21 D20] * et[-5 -6 D31 D30 D11 D10] * + el[D21 D20 D41 D40 -1 -2] * conj(bra[D31 -7 -3 D41]) * ket[D30 -8 -4 D40] + else + ket2 = (T <: PEPOTensor) ? first(fuse_physicalspaces(ket)) : twistdual(ket, 1) + bra2 = (T <: PEPOTensor) ? first(fuse_physicalspaces(bra)) : bra + return @tensoropt ctl2[:] := ctl[D11 D10 D21 D20] * et[-5 -6 D31 D30 D11 D10] * + el[D21 D20 D41 D40 -1 -2] * conj(bra2[d D31 -7 -3 D41]) * ket2[d D30 -8 -4 D40] + end +end + +""" +Enlarge the southeast corner +``` + -1/-2 -3/-4 + ║ ║ + -5/-6 ═════ Y ══ D1 ═══ er + ║ ║ + D2 D3 + ║ ║ + -7/-8 ═════ eb ═ D4 ══ cbr +``` +""" +function enlarge_corner_se( + cbr::AbstractTensor{E, S, 4}, + eb::AbstractTensor{E, S, 6}, er::AbstractTensor{E, S, 6}, + ket::T, bra::T = ket + ) where {E, S, T <: StateTensor} + if T <: PEPSOrth + return @tensoropt cbr2[:] := cbr[D31 D30 D41 D40] * eb[D21 D20 D41 D40 -7 -8] * + er[-3 -4 D31 D30 D11 D10] * conj(bra[-1 D11 D21 -5]) * ket[-2 D10 D20 -6] + else + ket2 = (T <: PEPOTensor) ? first(fuse_physicalspaces(ket)) : twistdual(ket, 1) + bra2 = (T <: PEPOTensor) ? first(fuse_physicalspaces(bra)) : bra + return @tensoropt cbr2[:] := cbr[D31 D30 D41 D40] * eb[D21 D20 D41 D40 -7 -8] * + er[-3 -4 D31 D30 D11 D10] * conj(bra2[d -1 D11 D21 -5]) * ket2[d -2 D10 D20 -6] + end +end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index e771cefb6..6f1b21e76 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -27,6 +27,8 @@ struct InfinitePEPO{T <: PEPOTensor} end end +const InfiniteState = Union{InfinitePEPS, InfinitePEPO} + ## Constructors """ diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl new file mode 100644 index 000000000..cc2f76808 --- /dev/null +++ b/test/bondenv/benv_ntu.jl @@ -0,0 +1,38 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random + +Nr, Nc = 2, 2 +Random.seed!(20) +Pspace = Vect[FermionParity](0 => 1, 1 => 1) +V2 = Vect[FermionParity](0 => 1, 1 => 1) +V3 = Vect[FermionParity](0 => 1, 1 => 2) +V4 = Vect[FermionParity](0 => 2, 1 => 2) +V5 = Vect[FermionParity](0 => 3, 1 => 2) +W1 = Vect[FermionParity](0 => 2, 1 => 3) +W2 = Vect[FermionParity](0 => 4, 1 => 1) +Pspaces = fill(Pspace, (Nr, Nc)) +Nspaces = [V2 V2; V4 V4] +Espaces = [V3 V5; V5 V3] + +peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) +# The NTU bond environments are constructed exactly +# and should be positive definite +for env_alg in (NNEnv(), NNNEnv()) + @info "Testing $(typeof(env_alg))" + for row in 1:Nr, col in 1:Nc + cp1 = PEPSKit._next(col, Nc) + A, B = peps.A[row, col], peps.A[row, cp1] + X, a, b, Y = PEPSKit._qr_bond(A, B) + benv = PEPSKit.bondenv_ntu(row, col, X, Y, peps, env_alg) + # benv should be Hermitian + @test benv' ≈ benv + benv = (benv + benv') / 2 + # benv should be positive definite + D, U = eigh(benv) + @test all(all(x -> x >= 0, diag(b)) for (k, b) in blocks(D)) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index f18007942..520d178c2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -59,6 +59,9 @@ end @time @safetestset "Full update bond environment" begin include("bondenv/benv_fu.jl") end + @time @safetestset "Bond environments for NTU" begin + include("bondenv/benv_ntu.jl") + end end if GROUP == "ALL" || GROUP == "TIMEEVOL" @time @safetestset "Cluster truncation with projectors" begin From ecbbee4dd7e3554049589751931ca529937d3578 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 13 May 2025 15:48:07 +0800 Subject: [PATCH 2/2] Add `add_bwt` option to NTU bondenv --- test/bondenv/benv_ntu.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl index cc2f76808..9249e1cfd 100644 --- a/test/bondenv/benv_ntu.jl +++ b/test/bondenv/benv_ntu.jl @@ -8,12 +8,10 @@ using Random Nr, Nc = 2, 2 Random.seed!(20) Pspace = Vect[FermionParity](0 => 1, 1 => 1) -V2 = Vect[FermionParity](0 => 1, 1 => 1) -V3 = Vect[FermionParity](0 => 1, 1 => 2) +V2 = Vect[FermionParity](0 => 4, 1 => 1) +V3 = Vect[FermionParity](0 => 3, 1 => 2) V4 = Vect[FermionParity](0 => 2, 1 => 2) -V5 = Vect[FermionParity](0 => 3, 1 => 2) -W1 = Vect[FermionParity](0 => 2, 1 => 3) -W2 = Vect[FermionParity](0 => 4, 1 => 1) +V5 = Vect[FermionParity](0 => 2, 1 => 3) Pspaces = fill(Pspace, (Nr, Nc)) Nspaces = [V2 V2; V4 V4] Espaces = [V3 V5; V5 V3]