Add add_bwt option to NTU bondenv

Author: Yue-Zhengyuan

src/algorithms/contractions/bondenv/benv_ntu.jl

@@ -20,12 +20,14 @@ Construct the "NTU-NN" bond environment.
             (+1 +0)══(+1 +1)
 ```
 """
-struct NTUEnvNN <: NTUEnvAlgorithm end
+@kwdef struct NTUEnvNN <: NTUEnvAlgorithm
+    add_bwt::Bool = true
+end
 """
 Calculate the bond environment within "NTU-NN" approximation.
 """
 function bondenv_ntu(
-    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, ::NTUEnvNN
+    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, alg::NTUEnvNN
 ) where {T<:Union{PEPSTensor,PEPSOrth}}
     neighbors = [
         (-1, 0, [NORTH, WEST]),
@@ -35,7 +37,7 @@ function bondenv_ntu(
         (0, 2, [NORTH, EAST, SOUTH]),
         (-1, 1, [NORTH, EAST]),
     ]
-    m = collect_neighbors(peps, row, col, neighbors)
+    m = collect_neighbors(peps, row, col, neighbors, alg.add_bwt)
     #= contraction indices
 
                 (-1 +0) ══ Dn ══ (-1 +1)
@@ -75,12 +77,14 @@ Construct the "NTU-NNN" bond environment.
     (+1 -1)=(+1 +0)══(+1 +1)=(+1 +2)
 ```
 """
-struct NTUEnvNNN <: NTUEnvAlgorithm end
+@kwdef struct NTUEnvNNN <: NTUEnvAlgorithm
+    add_bwt::Bool = true
+end
 """
 Calculates the bond environment within "NTU-NNN" approximation.
 """
 function bondenv_ntu(
-    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, ::NTUEnvNNN
+    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, alg::NTUEnvNNN
 ) where {T<:Union{PEPSTensor,PEPSOrth}}
     neighbors = [
         (-1, -1, [NORTH, WEST]),
@@ -94,7 +98,7 @@ function bondenv_ntu(
         (-1, 1, [NORTH]),
         (-1, 0, [NORTH]),
     ]
-    m = collect_neighbors(peps, row, col, neighbors)
+    m = collect_neighbors(peps, row, col, neighbors, alg.add_bwt)
     #= left half
         (-1 -1)══════(-1 +0)═ -1/-2
             ║           ║
@@ -144,12 +148,14 @@ Construct the "NTU-NNN+" bond environment.
             (+2 -1) (+2 +0)  (+2 +1) (+2 +2)
 ```
 """
-struct NTUEnvNNNp <: NTUEnvAlgorithm end
+@kwdef struct NTUEnvNNNp <: NTUEnvAlgorithm
+    add_bwt::Bool = true
+end
 """
 Calculates the bond environment within "NTU-NNN+" approximation.
 """
 function bondenv_ntu(
-    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, ::NTUEnvNNNp
+    row::Int, col::Int, X::T, Y::T, peps::InfiniteWeightPEPS, alg::NTUEnvNNNp
 ) where {T<:Union{PEPSTensor,PEPSOrth}}
     EMPTY = Vector{Int}()
     neighbors = [
@@ -178,7 +184,7 @@ function bondenv_ntu(
         (2, 1, [EAST, SOUTH, WEST]),
         (2, 2, [EAST, SOUTH, WEST]),
     ]
-    m = collect_neighbors(peps, row, col, neighbors)
+    m = collect_neighbors(peps, row, col, neighbors, alg.add_bwt)
     #= left half
                 (-2 -1)      (-2 +0)
                     ║           ║   

src/algorithms/contractions/bondenv/benv_tools.jl

@@ -31,13 +31,18 @@ function collect_neighbors(
     row::Int,
     col::Int,
     neighbors::Vector{Tuple{Int,Int,Vector{Int}}},
+    add_bwt::Bool=true,
 )
     Nr, Nc = size(peps)
     axs = Tuple(1:4)
     return Dict(
         (x, y) => begin
             r, c = mod1(row + x, Nr), mod1(col + y, Nc)
-            sqrts = Tuple(!(ax in open_axs) for ax in 1:4)
+            sqrts = if add_bwt
+                Tuple(!(ax in open_axs) for ax in 1:4)
+            else
+                ntuple(Returns(true), 4)
+            end
             _absorb_weights(peps.vertices[r, c], peps.weights, r, c, axs, sqrts, false)
         end for (x, y, open_axs) in neighbors
     )

test/bondenv/benv_ntu.jl

@@ -8,33 +8,37 @@ 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]
 
 peps = InfiniteWeightPEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
-for I in CartesianIndices(peps.vertices)
-    peps.vertices[I] /= norm(peps.vertices[I], Inf)
+normalize!.(peps.vertices, Inf)
+for wt in peps.weights
+    wt.data[:] = normalize(rand(length(wt.data)), Inf)
 end
-# The NTU bond environments are constructed exactly
-# and should be positive definite
-for env_alg in (NTUEnvNN(), NTUEnvNNN(), NTUEnvNNNp())
-    @info "Testing $(typeof(env_alg))"
-    for row in 1:Nr, col in 1:Nc
-        cp1 = PEPSKit._next(col, Nc)
-        A, B = peps.vertices[row, col], peps.vertices[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 should be positive definite
-        D, U = eigh(benv)
-        @test all(all(x -> x >= 0, diag(b)) for (k, b) in blocks(D))
+for add_bwt in (true, false)
+    for env_alg in (NTUEnvNN(; add_bwt), NTUEnvNNN(; add_bwt), NTUEnvNNNp(; add_bwt))
+        @info "Testing $(typeof(env_alg))"
+        for row in 1:Nr, col in 1:Nc
+            cp1 = PEPSKit._next(col, Nc)
+            A, B = peps.vertices[row, col], peps.vertices[row, cp1]
+            X, a, b, Y = PEPSKit._qr_bond(A, B)
+            @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY]
+            benv = PEPSKit.bondenv_ntu(row, col, X, Y, peps, env_alg)
+            # NTU bond environments are constructed exactly
+            # and should be positive definite
+            @test benv' ≈ benv
+            @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1]
+            nrm1 = PEPSKit.inner_prod(benv, ab, ab)
+            @test nrm1 ≈ real(nrm1)
+            D, U = eigh(benv)
+            @test all(all(x -> x >= 0, diag(b)) for (k, b) in blocks(D))
+            @assert benv ≈ U * D * U'
+        end
     end
 end