Add full update algorithm

examples/hubbard_su.jl

@@ -4,16 +4,16 @@ using PEPSKit
 using TensorKit
 
 # random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
-Dcut, symm = 8, Trivial
-N1, N2 = 2, 2
+Dbond, symm = 8, Trivial
+Nr, Nc = 2, 2
 Random.seed!(10)
 if symm == Trivial
     Pspace = Vect[fℤ₂](0 => 2, 1 => 2)
-    Vspace = Vect[fℤ₂](0 => Dcut / 2, 1 => Dcut / 2)
+    Vspace = Vect[fℤ₂](0 => Dbond / 2, 1 => Dbond / 2)
 else
     error("Not implemented")
 end
-peps = InfiniteWeightPEPS(rand, Float64, Pspace, Vspace; unitcell=(N1, N2))
+peps = InfiniteWeightPEPS(rand, Float64, Pspace, Vspace; unitcell=(Nr, Nc))
 
 # normalize vertex tensors
 for ind in CartesianIndices(peps.vertices)
@@ -22,14 +22,14 @@ end
 
 # Hubbard model Hamiltonian at half-filling
 t, U = 1, 6
-ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(N1, N2); t, U, mu=U / 2)
+ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(Nr, Nc); t, U, mu=U / 2)
 
 # simple update
 dts = [1e-2, 1e-3, 4e-4, 1e-4]
 tols = [1e-6, 1e-8, 1e-8, 1e-8]
 maxiter = 10000
 for (n, (dt, tol)) in enumerate(zip(dts, tols))
-    trscheme = truncerr(1e-10) & truncdim(Dcut)
+    trscheme = truncerr(1e-10) & truncdim(Dbond)
     alg = SimpleUpdate(dt, tol, maxiter, trscheme)
     peps, = simpleupdate(peps, ham, alg; bipartite=false)
 end
@@ -53,7 +53,7 @@ Es_exact = Dict(0 => -1.62, 2 => -0.176, 4 => 0.8603, 6 => -0.6567, 8 => -0.5243
 E_exact = Es_exact[U] - U / 2
 
 # measure energy
-E = costfun(peps, envs, ham) / (N1 * N2)
+E = costfun(peps, envs, ham) / (Nr * Nc)
 @info "Energy           = $E"
 @info "Benchmark energy = $E_exact"
 @test isapprox(E, E_exact; atol=5e-2)

src/PEPSKit.jl

@@ -3,14 +3,14 @@ module PEPSKit
 using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Base: @kwdef
 using Compat
-using Accessors: @set
 using VectorInterface
 using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations
 using ChainRulesCore, Zygote
 using LoggingExtras
 using MPSKit: loginit!, logiter!, logfinish!, logcancel!
 using MPSKitModels
 using FiniteDifferences
+using Accessors: @set
 using OhMyThreads: tmap
 
 include("utility/util.jl")
@@ -49,7 +49,9 @@ include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 
 include("algorithms/time_evolution/gatetools.jl")
+include("algorithms/time_evolution/bondenv.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
+include("algorithms/time_evolution/fullupdate.jl")
 
 include("algorithms/toolbox.jl")
 
@@ -193,6 +195,7 @@ export fixedpoint
 
 export absorb_weight
 export su_iter, simpleupdate, SimpleUpdate
+export fu_iter, fullupdate, FullUpdate, ALSOptimize
 
 export InfinitePEPS, InfiniteTransferPEPS
 export SUWeight, InfiniteWeightPEPS

src/algorithms/contractions/ctmrg_contractions.jl

@@ -230,7 +230,7 @@ Right projector:
 ```
 """
 function contract_projectors(U, S, V, Q, Q_next)
-    isqS = sdiag_pow(S, -0.5)
+    isqS = sdiag_inv_sqrt(S)
     P_left = Q_next * V' * isqS  # use * to respect fermionic case
     P_right = isqS * U' * Q
     return P_left, P_right

src/algorithms/ctmrg/sequential.jl

@@ -30,18 +30,36 @@ function SequentialCTMRG(;
     )
 end
 
+"""
+    ctmrg_leftmove(col::Int, state, envs::CTMRGEnv, alg::SequentialCTMRG)
+
+Perform sequential CTMRG left move on the `col`-th column.
+"""
+function ctmrg_leftmove(col::Int, state, envs::CTMRGEnv, alg::SequentialCTMRG)
+    #=
+        ----> left move
+        C1 ← T1 ←   r-1
+        ↓    ‖
+        T4 = M ==   r
+        ↓    ‖
+        C4 → T3 →   r+1
+        c-1  c 
+    =#
+    projectors, info = sequential_projectors(col, state, envs, alg.projector_alg)
+    envs = renormalize_sequentially(col, projectors, state, envs)
+    return envs, info
+end
+
 function ctmrg_iteration(state, envs::CTMRGEnv, alg::SequentialCTMRG)
     ϵ = zero(real(scalartype(state)))
     for _ in 1:4 # rotate
         for col in 1:size(state, 2) # left move column-wise
-            projectors, info = sequential_projectors(col, state, envs, alg.projector_alg)
-            envs = renormalize_sequentially(col, projectors, state, envs)
+            envs, info = ctmrg_leftmove(col, state, envs, alg)
             ϵ = max(ϵ, info.err)
         end
         state = rotate_north(state, EAST)
         envs = rotate_north(envs, EAST)
     end
-
     return envs, (; err=ϵ)
 end
 

src/algorithms/ctmrg/sparse_environments.jl

@@ -104,7 +104,7 @@
 
 # Compute left and right projectors sparsely without constructing enlarged corners explicitly 
 function left_and_right_projector(U, S, V, Q::EnlargedCorner, Q_next::EnlargedCorner)
-    isqS = sdiag_pow(S, -0.5)
+    isqS = sdiag_inv_sqrt(S)
     P_left = left_projector(Q.E_1, Q.C, Q.E_2, V, isqS, Q.ket, Q.bra)
     P_right = right_projector(
         Q_next.E_1, Q_next.C, Q_next.E_2, U, isqS, Q_next.ket, Q_next.bra

src/algorithms/time_evolution/bondenv.jl

@@ -0,0 +1,15 @@
+"""
+    BondEnvAlgorithm
+
+Abstract super type for all algorithms to construct the environment of a bond in the InfinitePEPS.
+"""
+abstract type BondEnvAlgorithm end
+
+const Hair{S} = AbstractTensorMap{S,1,1} where {S<:ElementarySpace}
+const BondEnv{S} = AbstractTensorMap{S,2,2} where {S<:ElementarySpace}
+const PEPSOrth{S} = AbstractTensor{S,4} where {S<:ElementarySpace}
+const BondPhys{S} = AbstractTensor{S,3} where {S<:ElementarySpace}
+const BondPhys2{S} = AbstractTensor{S,4} where {S<:ElementarySpace}
+
+include("bondenv/env_ctm.jl")
+include("bondenv/optimize.jl")

src/algorithms/time_evolution/bondenv/env_ctm.jl

@@ -0,0 +1,153 @@
+"""
+Construct the environment (norm) tensor
+```
+    left half                       right half
+    C1 -χ4 - T1 ------- χ6 ------- T1 - χ8 - C2     r-1
+    |        ‖                      ‖        |
+    χ2      DNX                    DNY      χ10
+    |        ‖                      ‖        |
+    T4 =DWX= XX = DX =       = DY = YY =DEY= T2     r
+    |        ‖                      ‖        |
+    χ1      DSX                    DSY       χ9
+    |        ‖                      ‖        |
+    C4 -χ3 - T3 ------- χ5 ------- T3 - χ7 - C3     r+1
+    c-1      c                      c+1     c+2
+```
+which can be more simply denoted as
+```
+    |------------|
+    |→ DX1  DY1 ←|   axis order
+    |← DX0  DX1 →|   (DX1, DY1, DX0, DY0)
+    |------------|
+```
+The axes 1, 2 (or 3, 4) come from X†, Y† (or X, Y)
+"""
+function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, envs::CTMRGEnv)
+    Nr, Nc = size(envs.corners)[[2, 3]]
+    cm1 = _prev(col, Nc)
+    cp1 = _next(col, Nc)
+    cp2 = _next(cp1, Nc)
+    rm1 = _prev(row, Nr)
+    rp1 = _next(row, Nr)
+    c1 = envs.corners[1, rm1, cm1]
+    c2 = envs.corners[2, rm1, cp2]
+    c3 = envs.corners[3, rp1, cp2]
+    c4 = envs.corners[4, rp1, cm1]
+    t1X, t1Y = envs.edges[1, rm1, col], envs.edges[1, rm1, cp1]
+    t2 = envs.edges[2, row, cp2]
+    t3X, t3Y = envs.edges[3, rp1, col], envs.edges[3, rp1, cp1]
+    t4 = envs.edges[4, row, cm1]
+    # left half
+    @autoopt @tensor lhalf[DX1, DX0, χ5, χ6] := (
+        c4[χ3, χ1] *
+        t4[χ1, DWX0, DWX1, χ2] *
+        c1[χ2, χ4] *
+        t3X[χ5, DSX0, DSX1, χ3] *
+        X[DNX0, DX0, DSX0, DWX0] *
+        conj(X[DNX1, DX1, DSX1, DWX1]) *
+        t1X[χ4, DNX0, DNX1, χ6]
+    )
+    # right half
+    @autoopt @tensor rhalf[DY1, DY0, χ5, χ6] := (
+        c3[χ9, χ7] *
+        t2[χ10, DEY0, DEY1, χ9] *
+        c2[χ8, χ10] *
+        t3Y[χ7, DSY0, DSY1, χ5] *
+        Y[DNY0, DEY0, DSY0, DY0] *
+        conj(Y[DNY1, DEY1, DSY1, DY1]) *
+        t1Y[χ6, DNY0, DNY1, χ8]
+    )
+    # combine
+    @autoopt @tensor env[DX1, DY1; DX0, DY0] := (
+        lhalf[DX1, DX0, χ5, χ6] * rhalf[DY1, DY0, χ5, χ6]
+    )
+    return env / norm(env, Inf)
+end
+
+"""
+Replace `env` by its positive/negative approximant `± Z Z†`
+(returns the sign and Z†)
+```
+                        |-→ 1   2 ←-|
+                        |           |
+    |----env----|       |←--- Z ---→|
+    |→ 1     2 ←|   =         ↑
+    |← 3     4 →|       |---→ Z† ←--|
+    |-----------|       |           |
+                        |←- 3   4 -→|
+```
+"""
+function positive_approx(env::BondEnv)
+    @assert [isdual(space(env, ax)) for ax in 1:4] == [0, 0, 1, 1]
+    # hermitize env, and perform eigen-decomposition
+    # env = U D U'
+    D, U = eigh((env + env') / 2)
+    # determine env is (mostly) positive or negative
+    sgn = sign(mean(vcat((diag(b) for (k, b) in blocks(D))...)))
+    if sgn == -1
+        D *= -1
+    end
+    # set negative eigenvalues to 0
+    for (k, b) in blocks(D)
+        for i in diagind(b)
+            if b[i] < 0
+                b[i] = 0.0
+            end
+        end
+    end
+    Zdg = sdiag_pow(D, 0.5) * U'
+    return sgn, Zdg
+end
+
+"""
+Fix local gauge of the env tensor around a bond
+"""
+function fu_fixgauge(
+    Zdg::AbstractTensorMap{S,1,2}, X::PEPSOrth, Y::PEPSOrth, aR::BondPhys, bL::BondPhys
+) where {S<:ElementarySpace}
+    #= 
+            1               1
+            ↑               ↑
+        2 → Z† ← 3  =   2 → QR ← 3  1 ← R ← 2
+
+                                        1
+                                        ↑
+                    =   2 → L → 1   3 → QL ← 2
+    =#
+    QR, R = leftorth(Zdg, ((1, 2), (3,)); alg=QRpos())
+    QL, L = leftorth(Zdg, ((1, 3), (2,)); alg=QRpos())
+    @assert !isdual(codomain(R)[1]) && !isdual(domain(R)[1])
+    @assert !isdual(codomain(L)[1]) && !isdual(domain(L)[1])
+    Rinv, Linv = inv(R), inv(L)
+    #= fix gauge of aR, bL, Z†
+
+                    ↑
+        |→-(Linv -→ Z† ← Rinv)←-|
+        |                       |
+        ↑                       ↑
+        |        ↑     ↑        |
+        |← (L ← aR) ← (bL → R) →|
+        |-----------------------|
+
+                     -2              -2
+                      ↑               ↑        
+        -1 ← L ← 1 ← aR2 ← -3   -1 ← bL2 → 1 → R → -3
+
+                        -1
+                        ↑
+        -2 → Linv → 1 → Z† ← 2 ← Rinv ← -3
+    =#
+    aR = ncon([L, aR], [[-1, 1], [1, -2, -3]])
+    bL = ncon([bL, R], [[-1, -2, 1], [-3, 1]])
+    Zdg = permute(ncon([Zdg, Linv, Rinv], [[-1, 1, 2], [1, -2], [2, -3]]), (1,), (2, 3))
+    #= fix gauge of X, Y
+            -1                                      -1
+             |                                      |
+        -4 - X ← 1 ← Linv ← -2      -4 → Rinv → 1 → Y - -2
+             |                                      |
+            -3                                      -3
+    =#
+    X = ncon([X, Linv], [[-1, 1, -3, -4], [1, -2]])
+    Y = ncon([Y, Rinv], [[-1, -2, -3, 1], [1, -4]])
+    return Zdg, X, Y, aR, bL
+end