Do not move phys leg to bond tensor for middle cluster sites

Author: Yue-Zhengyuan

src/algorithms/contractions/bondenv/benv_ctm.jl

@@ -24,8 +24,8 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in
 ```
 """
 function bondenv_ctm(
-        row::Int, col::Int, X::T, Y::T, env::CTMRGEnv
-    ) where {T}
+        row::Int, col::Int, X::TX, Y::TY, env::CTMRGEnv
+    ) where {TX, TY}
     Nr, Nc = size(env.corners)[[2, 3]]
     cm1 = _prev(col, Nc)
     cp1 = _next(col, Nc)

src/algorithms/contractions/bondenv/benv_ntu.jl

@@ -34,8 +34,8 @@ 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, S <: InfiniteState}
+        row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNEnv
+    ) where {TX, TY, S <: InfiniteState}
     neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)]
     m = collect_neighbors(state, row, col, neighbors)
     X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y)
@@ -78,8 +78,8 @@ Calculate the bond environment within "NTU-NN+" approximation.
 Dotted lines and ○ are splitted using SVD with `truncrank(1)`.
 """
 function bondenv_ntu(
-        row::Int, col::Int, X::T, Y::T, state::S, alg::NNpEnv
-    ) where {T, S <: InfiniteState}
+        row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNpEnv
+    ) where {TX, TY, S <: InfiniteState}
     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),
@@ -191,8 +191,8 @@ 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, S <: InfiniteState}
+        row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNNEnv
+    ) where {TX, TY, S <: InfiniteState}
     neighbors = [
         (-1, -1), (0, -1), (1, -1),
         (1, 0), (1, 1), (1, 2), (0, 2),

src/algorithms/contractions/bondenv/gaugefix.jl

@@ -90,40 +90,59 @@ function fixgauge_benv(
 end
 
 """
-When the (half) bond environment `Z` consists of 
-two `PEPSOrth` or `PEPOOrth` tensors `X`, `Y` as
+Apply the gauge transformation `Rinv` for `Z`
 ```
     ┌-----------------------┐
-    |                       |
-    └---Z---(X)--   --(Y)---┘
+    └---Z--(X)--Rinv--   ---┘
         ↓
 ```
-apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`:
+to `X`. For example, when `X` is a `PEPSTensor`,
 ```
-        -1                                     -1
-         |                                      |
-    -4 - X - 1 - Rinv - -2      -4 - Linv - 1 - Y - -2
-         |                                      |
-        -3                                     -3
-    
-        -2                                     -2
-         |                                      |
-    -5 - X - 1 - Rinv - -3      -5 - Linv - 1 - Y - -3
-         | ╲                                    | ╲
-        -4  -1                                 -4  -1
+        -2
+         |
+    -5 - X - 1 - Rinv - -3
+         | ╲ 
+        -4  -1
 ```
 """
-function _fixgauge_benvXY(
-        X::PEPSOrth, Y::PEPSOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor,
-    )
-    @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2]
-    @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4]
-    return X, Y
+function _fixgauge_benvX(X::PEPSOrth, Rinv::MPSBondTensor)
+    return @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2]
 end
-function _fixgauge_benvXY(
-        X::PEPOOrth, Y::PEPOOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor,
-    )
-    @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3]
-    @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5]
-    return X, Y
+function _fixgauge_benvX(X::PEPSTensor, Rinv::MPSBondTensor)
+    return @plansor X[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * Rinv[1; -3]
+end
+function _fixgauge_benvX(X::PEPOOrth, Rinv::MPSBondTensor)
+    return @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3]
+end
+function _fixgauge_benvX(X::PEPOTensor, Rinv::MPSBondTensor)
+    return @plansor X[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * Rinv[1; -4]
+end
+
+"""
+Apply the gauge transformation `Linv` for `Z`
+```
+    ┌-----------------------┐
+    └---Z---  ---Linv--(Y)--┘
+        ↓
+```
+to `Y`. For example, when `Y` is a `PEPSTensor`,
+```
+                   -2
+                    |
+    -5 - Linv - 1 - Y - -3
+                    | ╲
+                   -4  -1
+```
+"""
+function _fixgauge_benvY(Y::PEPSOrth, Linv::MPSBondTensor)
+    return @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4]
+end
+function _fixgauge_benvY(Y::PEPSTensor, Linv::MPSBondTensor)
+    return @plansor Y[-1; -2 -3 -4 -5] := Y[-1; -2 -3 -4 1] * Linv[1; -5]
+end
+function _fixgauge_benvY(Y::PEPOOrth, Linv::MPSBondTensor)
+    return @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5]
+end
+function _fixgauge_benvY(Y::PEPOTensor, Linv::MPSBondTensor)
+    return @plansor Y[-1 -2; -3 -4 -5 -6] := Y[-1 -2; -3 -4 -5 1] * Linv[1; -6]
 end

src/algorithms/time_evolution/ntupdate.jl

@@ -79,7 +79,7 @@ function _ntu_iter(
     Nr, Nc = size(state)
     trunc = only(_get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc)))
     alg′ = (@set alg.opt_alg.trunc = trunc)
-    return _bond_truncate(state, wts, Tuple(sites), alg′; gate)
+    return _bond_truncate(state, wts, Tuple(sites), (:first, :last), alg′; gate)
 end
 
 """

src/algorithms/time_evolution/ntupdate3site.jl

@@ -23,9 +23,13 @@ function _ntu_iter(
 
     # truncate each bond sequentially along the path
     info = (; fid = 1.0)
-    for (bondsites, trunc) in zip(zip(sites, Iterators.drop(sites, 1)), truncs)
+    nbond = length(sites) - 1
+    for (i, bondsites) in enumerate(zip(sites, Iterators.drop(sites, 1)))
+        trunc = truncs[i]
         alg′ = (@set alg.opt_alg.trunc = trunc)
-        state, wts, info′ = _bond_truncate(state, wts, bondsites, alg′)
+        stype1 = (i == 1) ? :first : :middle
+        stype2 = (i == nbond) ? :last : :middle
+        state, wts, info′ = _bond_truncate(state, wts, bondsites, (stype1, stype2), alg′)
         # record the worst fidelity
         (info′.fid < info.fid) && (info = info′)
     end
@@ -35,10 +39,14 @@ end
 """
 Truncate a nearest neighbor bond between `site1` and `site2`
 after rotating the bond to standard x direction `A ← B`.
+
+`bondtype` takes values in (1, 2, 3), meaning that the current bond is
+(the first, a middle, the last) bond in the updated cluster.
 """
 function _bond_truncate(
         state::InfiniteState, wts::SUWeight,
         (site1, site2)::NTuple{2, CartesianIndex{2}},
+        (stype1, stype2)::NTuple{2, Symbol},
         alg::NeighbourUpdate; gate::Union{NNGate, Nothing} = nothing
     )
     # rotate bond to standard x direction `A ← B`
@@ -54,14 +62,26 @@ function _bond_truncate(
     A, B = state2[row, col], state2[row, cp1]
 
     # create bond environment
-    a, X = bond_tensor_first(A; trunc = trunctol(; rtol = 1.0e-12))
-    b, Y = bond_tensor_last(B; trunc = trunctol(; rtol = 1.0e-12))
+    qrtrunc = trunctol(; rtol = 1.0e-12)
+    a, X = if stype1 == :first
+        bond_tensor_first(A; trunc = qrtrunc)
+    else
+        @assert stype1 == :middle
+        bond_tensor_midnext(A; trunc = qrtrunc)
+    end
+    b, Y = if stype2 == :last
+        bond_tensor_last(B; trunc = qrtrunc)
+    else
+        @assert stype2 == :middle
+        bond_tensor_midprev(B; trunc = qrtrunc)
+    end
     benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg)
     @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))"
     if alg.fixgauge
         Z = positive_approx(benv)
         Z, a, b, (Linv, Rinv) = fixgauge_benv(Z, a, b)
-        X, Y = _fixgauge_benvXY(X, Y, Linv, Rinv)
+        X = _fixgauge_benvX(X, Rinv)
+        Y = _fixgauge_benvY(Y, Linv)
         benv = Z' * Z
         @debug "cond(L) = $(LinearAlgebra.cond(Linv)); cond(R): $(LinearAlgebra.cond(Rinv))"
         @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))"
@@ -70,14 +90,20 @@ function _bond_truncate(
     # (optional) apply the NN gate without truncation
     if !(gate === nothing)
         a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15))
+    end
+    a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg)
+
+    A = if stype1 == :first
+        undo_bond_tensor_first(a, X)
     else
-        a = permute(a, ((1, 2), (3,)))
-        b = permute(b, ((1,), (2, 3)))
+        undo_bond_tensor_midnext(a, X)
+    end
+    B = if stype2 == :last
+        undo_bond_tensor_last(b, Y)
+    else
+        undo_bond_tensor_midprev(b, Y)
     end
 
-    a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg)
-    A = undo_bond_tensor_first(a, X)
-    B = undo_bond_tensor_last(b, Y)
     normalize!(A, Inf)
     normalize!(B, Inf)
     normalize!(s, Inf)

src/algorithms/truncation/bond_tensor.jl

@@ -113,3 +113,101 @@ function undo_bond_tensor_last(b::MPSTensor, Y::PEPOOrth; gate_ax::Integer = 1)
         return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1]
     end
 end
+
+"""
+Given a middle tensor `A` in the cluster acted on by a gate,
+obtain reduced tensor on its next bond.
+
+For PEPSTensor,
+```
+        2
+        |
+    5 - X ← 3   1 ← a - 3
+        | ↘          ↘
+        4   1         (2)
+```
+For PEPOTensor,
+```
+    2   3
+      ↘ |
+    6 - X ← 4   1 ← a - 3
+        | ↘          ↘
+        5   1         (2)
+```
+
+Here the physical leg on `a` is an auxiliary trivial leg.
+"""
+function bond_tensor_midnext(A::PEPSTensor; kwargs...)
+    X, a = left_orth!(permute(A, ((1, 2, 4, 5), (3,)); copy = true); kwargs...)
+    X = permute(X, ((1,), (2, 5, 3, 4)))
+    a = insertrightunit(a, 1)
+    return a, X
+end
+function bond_tensor_midnext(A::PEPOTensor; kwargs...)
+    X, a = left_orth!(permute(A, ((1, 2, 3, 5, 6), (4,)); copy = true); kwargs...)
+    X = permute(X, ((1, 2), (3, 6, 4, 5)))
+    a = insertrightunit(a, 1)
+    return a, X
+end
+
+"""
+Undo the decomposition in `bond_tensor_midprev`.
+"""
+function undo_bond_tensor_midnext(a::MPSTensor, X::PEPSTensor)
+    a = removeunit(a, 2)
+    return @tensor A[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * a[1; -3]
+end
+function undo_bond_tensor_midnext(a::MPSTensor, X::PEPOTensor)
+    a = removeunit(a, 2)
+    return @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * a[1; -4]
+end
+
+"""
+Given a middle tensor `A` in the cluster acted on by a gate,
+obtain reduced tensor on its previous bond.
+
+For PEPSTensor,
+```
+                    2
+                    |
+    1 - b → 3   5 → Y - 3
+          ↘         | ↘  
+           (2)      4   1
+```
+For PEPOTensor,
+```
+                2   3
+                  ↘ |
+    1 - b → 3   6 → Y - 4
+          ↘         | ↘  
+           (2)      5   1
+```
+
+Here the physical leg on `a` is an auxiliary trivial leg.
+"""
+function bond_tensor_midprev(A::PEPSTensor; kwargs...)
+    Y, b = left_orth!(permute(A, ((1, 2, 3, 4), (5,)); copy = true); kwargs...)
+    Y = permute(Y, ((1,), (2, 3, 4, 5)))
+    b = permute(b, ((2,), (1,)))
+    b = insertrightunit(b, 1)
+    return b, Y
+end
+function bond_tensor_midprev(A::PEPOTensor; kwargs...)
+    Y, b = left_orth!(permute(A, ((1, 2, 3, 4, 5), (6,)); copy = true); kwargs...)
+    Y = permute(Y, ((1, 2), (3, 4, 5, 6)))
+    b = permute(b, ((2,), (1,)))
+    b = insertrightunit(b, 1)
+    return b, Y
+end
+
+"""
+Undo the decomposition in `bond_tensor_midprev`.
+"""
+function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPSTensor)
+    b = removeunit(b, 2)
+    return @tensor A[-1; -2 -3 -4 -5] := b[-5; 1] * Y[-1; -2 -3 -4 1]
+end
+function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPOTensor)
+    b = removeunit(b, 2)
+    return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6; 1] * Y[-1 -2; -3 -4 -5 1]
+end

test/bondenv/benv_gaugefix.jl

@@ -37,7 +37,8 @@ for V1 in Vs, V2 in Vs, V3 in Vs
     @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3]
     @test half ≈ half2
     # test gauge transformation of X, Y
-    X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv)
+    X2 = _fixgauge_benvX(X, Rinv)
+    Y2 = _fixgauge_benvY(Y, Linv)
     @tensor Z2_[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X2[Xn Xe Xs Xw] * Y2[Yn Ye Ys Yw]
     @test Z2 ≈ Z2_
 end

test/bondenv/benv_ntu.jl

@@ -41,7 +41,8 @@ function test_ntu_env(
         # verify gauge transformation of X, Y
         @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY]
         nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2)
-        X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv)
+        X2 = _fixgauge_benvX(X, Rinv)
+        Y2 = _fixgauge_benvY(Y, Linv)
         benv3 = PEPSKit.bondenv_ntu(row, col, X2, Y2, state, env_alg)
         benv3 *= norm(benv2, Inf)
         nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2)