Allow bond-dependent truncation scheme in Simple Update

src/PEPSKit.jl

@@ -82,7 +82,8 @@ using .Defaults: set_scheduler!
 export set_scheduler!
 export SVDAdjoint, FullSVDReverseRule, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
-export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector
+export FixedSpaceTruncation, SiteDependentTruncation
+export HalfInfiniteProjector, FullInfiniteProjector
 export LocalOperator, physicalspace
 export expectation_value, cost_function, product_peps, correlation_length, network_value
 export correlator

src/algorithms/time_evolution/simpleupdate.jl

@@ -34,7 +34,7 @@
     col::Int,
     gate::AbstractTensorMap{T,S,2,2},
     peps::InfiniteWeightPEPS,
-    alg::SimpleUpdate,
+    trscheme::TensorKit.TruncationScheme,
 ) where {T<:Number,S<:ElementarySpace}
     Nr, Nc = size(peps)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
@@ -47,7 +47,7 @@
     B = _absorb_weights(B, peps.weights, row, cp1, Tuple(1:4), sqrtsB, false)
     # apply gate
     X, a, b, Y = _qr_bond(A, B)
-    a, s, b, ϵ = _apply_gate(a, b, gate, alg.trscheme)
+    a, s, b, ϵ = _apply_gate(a, b, gate, trscheme)
     A, B = _qr_bond_undo(X, a, b, Y)
     # remove environment weights
     _allfalse = ntuple(Returns(false), 3)
@@ -94,7 +94,13 @@
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, 1), CartesianIndex(r, 2)),
                 )
-                ϵ = _su_bondx!(r, 1, term, peps2, alg)
+                ϵ = _su_bondx!(
+                    r,
+                    1,
+                    term,
+                    peps2,
+                    truncation_scheme(alg.trscheme; direction, r, c, mirror_antidiag=true),
+                )
                 peps2.vertices[rp1, 2] = deepcopy(peps2.vertices[r, 1])
                 peps2.vertices[rp1, 1] = deepcopy(peps2.vertices[r, 2])
                 peps2.weights[1, rp1, 2] = deepcopy(peps2.weights[1, r, 1])
@@ -106,7 +112,13 @@
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, c), CartesianIndex(r, c + 1)),
                 )
-                ϵ = _su_bondx!(r, c, term, peps2, alg)
+                ϵ = _su_bondx!(
+                    r,
+                    c,
+                    term,
+                    peps2,
+                    truncation_scheme(alg.trscheme; direction, r, c, mirror_antidiag=true),
+                )
             end
         end
         if direction == 2

src/algorithms/time_evolution/simpleupdate3site.jl

@@ -373,7 +373,11 @@
 end
 
 function _su3site_se!(
-    row::Int, col::Int, gs::Vector{T}, peps::InfiniteWeightPEPS, alg::SimpleUpdate
+    row::Int,
+    col::Int,
+    gs::Vector{T},
+    peps::InfiniteWeightPEPS,
+    trscheme::TensorKit.TruncationScheme,
 ) where {T<:AbstractTensorMap}
     Nr, Nc = size(peps)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
@@ -384,7 +388,7 @@
     coords = ((row, col), (row, cp1), (rm1, cp1))
     # weights in the cluster
     wt_idxs = ((1, row, col), (2, row, cp1))
-    wts, ϵ = apply_gatempo!(Ms, gs; trunc=alg.trscheme)
+    wts, ϵ = apply_gatempo!(Ms, gs; trunc=trscheme)
     for (wt, wt_idx) in zip(wts, wt_idxs)
         peps.weights[CartesianIndex(wt_idx)] = wt / norm(wt, Inf)
     end
@@ -418,7 +422,13 @@
         for site in CartesianIndices(peps2.vertices)
             r, c = site[1], site[2]
             gs = gatempos[i][r, c]
-            _su3site_se!(r, c, gs, peps2, alg)
+            _su3site_se!(
+                r,
+                c,
+                gs,
+                peps2,
+                truncation_scheme(alg.trscheme, direction, r, c; mirror_antidiag=true),
+            )
         end
         peps2 = (i == 1) ? rotl90(peps2) : rotr90(peps2)
     end

src/algorithms/truncation/truncationschemes.jl

@@ -7,6 +7,22 @@
 """
 struct FixedSpaceTruncation <: TensorKit.TruncationScheme end
 
+struct SiteDependentTruncation <: TensorKit.TruncationScheme
+    trschemes::Array{T,3} where {T<:TensorKit.TruncationScheme}
+end
+
+function SiteDependentTruncation(trscheme::TensorKit.TruncationScheme, Nr::Int, Nc::Int)
+    return SiteDependentTruncation(reshape(fill(trscheme, Nr, Nc), 2, Nr, Nc))
+end
+
+function SiteDependentTruncation(
+    trschemes::Tuple{T,S}, Nr::Int, Nc::Int
+) where {T<:TensorKit.TruncationScheme,S<:TensorKit.TruncationScheme}
+    return SiteDependentTruncation(
+        reshape([trschemes[mod1(i, 2)] for i in 1:(2 * Nr * Nc)], 2, Nr, Nc)
+    )
+end
+
 const TRUNCATION_SCHEME_SYMBOLS = IdDict{Symbol,Type{<:TruncationScheme}}(
     :fixedspace => FixedSpaceTruncation,
     :notrunc => TensorKit.NoTruncation,
@@ -25,3 +41,22 @@
 
     return isnothing(η) ? alg_type() : alg_type(η)
 end
+
+function truncation_scheme(trscheme::T; kwargs...) where {T<:TensorKit.TruncationScheme}
+    return trscheme
+end
+
+function truncation_scheme(
+    trscheme::SiteDependentTruncation;
+    direction::Int,
+    r::Int,
+    c::Int,
+    mirror_antidiag::Bool=false,
+)
+    if mirror_antidiag && direction == 2
+        depth, Nr, Nc = size(trscheme.trschemes)
+        @assert depth == 2
+        return trscheme.trschemes[direction, Nc - c + 1, Nr - r + 1]
+    end
+    return trscheme.trschemes[direction, r, c]
+end

test/examples/heisenberg.jl

@@ -71,7 +71,32 @@ end
     maxiter = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
-        trscheme = truncerr(1e-10) & truncdim(Dbond2)
+        if n == 2
+            trscheme = SiteDependentTruncation(
+                (truncerr(tol) & truncdim(Dbond2), truncerr(tol) & truncdim(Dbond2 + 1)),
+                N1,
+                N2,
+            )
+        elseif n == 3
+            trscheme = SiteDependentTruncation(
+                reshape(
+                    [truncerr(tol) & truncdim(Dbond2) truncerr(tol) & truncdim(Dbond2 + 1) truncerr(
+                        tol
+                    ) & truncdim(
+                        Dbond2
+                    ) truncerr(tol) & truncdim(Dbond2 + 1) truncerr(tol) & truncdim(Dbond2) truncerr(
+                        tol
+                    ) & truncdim(
+                        Dbond2
+                    ) truncerr(tol) & truncdim(Dbond2) truncerr(tol) & truncdim(Dbond2)],
+                    2,
+                    2,
+                    2,
+                ),
+            )
+        else
+            trscheme = truncerr(tol) & truncdim(Dbond2)
+        end
         alg = SimpleUpdate(dt, tol, maxiter, trscheme)
         result = simpleupdate(wpeps, ham, alg; bipartite=false)
         wpeps = result[1]