Add Neighborhood Tensor Update (aka Cluster Update)

.github/workflows/Tests.yml

@@ -30,6 +30,7 @@ jobs:
           - boundarymps
           - examples
           - utility
+          - bondenv
         os:
           - ubuntu-latest
           - macOS-latest

examples/heisenberg_evol/heis_tools.jl

@@ -0,0 +1,61 @@
+@static if Sys.isapple()
+    using AppleAccelerate
+end
+using Test
+using Printf
+using Random
+import Statistics: mean
+using TensorKit
+using PEPSKit
+
+module MeasureHeis
+
+export measure_heis
+
+using TensorKit
+import MPSKitModels: S_x, S_y, S_z, S_exchange
+using PEPSKit
+
+"""
+Measure magnetization on each site
+"""
+function cal_mags(peps::InfinitePEPS, envs::CTMRGEnv)
+    N1, N2 = size(peps)
+    lattice = collect(space(t, 1) for t in peps.A)
+    # detect symmetry on physical axis
+    symm = sectortype(space(peps.A[1, 1]))
+    if symm == Trivial
+        Sas = real.([S_x(symm), im * S_y(symm), S_z(symm)])
+    elseif symm == U1Irrep
+        # only Sz preserves <Sz>
+        Sas = real.([S_z(symm)])
+    else
+        throw(ArgumentError("Unrecognized symmetry on physical axis"))
+    end
+    return [
+        collect(
+            expectation_value(
+                peps, LocalOperator(lattice, (CartesianIndex(r, c),) => Sa), envs
+            ) for (r, c) in Iterators.product(1:N1, 1:N2)
+        ) for Sa in Sas
+    ]
+end
+
+"""
+Measure physical quantities for Heisenberg model
+"""
+function measure_heis(peps::InfinitePEPS, H::LocalOperator, envs::CTMRGEnv)
+    results = Dict{String,Any}()
+    N1, N2 = size(peps)
+    results["e_site"] = costfun(peps, envs, H) / (N1 * N2)
+    results["mag"] = cal_mags(peps, envs)
+    results["mag_norm"] = collect(
+        norm([mags[r, c] for mags in results["mag"]]) for
+        (r, c) in Iterators.product(1:N1, 1:N2)
+    )
+    return results
+end
+
+end
+
+import .MeasureHeis: measure_heis

examples/heisenberg_evol/ntupdate.jl

@@ -0,0 +1,29 @@
+include("simpleupdate.jl")
+
+# continue SU with NTU
+dts = [1e-2, 5e-3]
+maxiter = 2000
+for dt in dts
+    alg = NTUpdate(;
+        dt,
+        maxiter,
+        tol=1e-8,
+        bondenv_alg=NTUEnvNN(),
+        opt_alg=FullEnvTruncation(; trscheme=trscheme_peps),
+        # opt_alg=ALSTruncation(; trscheme=trscheme_peps),
+        ctm_alg=SequentialCTMRG(;
+            tol=1e-7, verbosity=2, maxiter=30, trscheme=trscheme_envs
+        ),
+    )
+    result = ntupdate(peps, envs, ham, alg, ctm_alg; bipartite=true)
+    global peps = result[1]
+    global envs = result[2]
+end
+# measure physical quantities with CTMRG
+peps_ = InfinitePEPS(peps)
+meas = measure_heis(peps_, ham, envs)
+display(meas)
+@info @sprintf("Energy = %.8f\n", meas["e_site"])
+@info @sprintf("Staggered magnetization = %.8f\n", mean(meas["mag_norm"]))
+@test isapprox(meas["e_site"], -0.66884; atol=1e-3)
+@test isapprox(mean(meas["mag_norm"]), 0.34; atol=1e-2)

examples/heisenberg_evol/simpleupdate.jl

@@ -0,0 +1,52 @@
+include("heis_tools.jl")
+
+# benchmark data is from Phys. Rev. B 94, 035133 (2016)
+Dbond, χenv, symm = 4, 16, Trivial
+trscheme_peps = truncerr(1e-10) & truncdim(Dbond)
+trscheme_envs = truncerr(1e-10) & truncdim(χenv)
+N1, N2 = 2, 2
+# Heisenberg model Hamiltonian
+# (already only includes nearest neighbor terms)
+ham = heisenberg_XYZ(ComplexF64, symm, InfiniteSquare(N1, N2); Jx=1.0, Jy=1.0, Jz=1.0)
+# convert to real tensors
+ham = LocalOperator(ham.lattice, Tuple(ind => real(op) for (ind, op) in ham.terms)...)
+
+# random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
+if symm == Trivial
+    Pspace = ℂ^2
+    Vspace = ℂ^Dbond
+    Espace = ℂ^χenv
+elseif symm == U1Irrep
+    Pspace = ℂ[U1Irrep](1//2 => 1, -1//2 => 1)
+    Vspace = ℂ[U1Irrep](0 => Dbond ÷ 2, 1//2 => Dbond ÷ 4, -1//2 => Dbond ÷ 4)
+    Espace = ℂ[U1Irrep](0 => χenv ÷ 2, 1//2 => χenv ÷ 4, -1//2 => χenv ÷ 4)
+else
+    error("Not implemented")
+end
+Random.seed!(0)
+peps = InfiniteWeightPEPS(rand, Float64, Pspace, Vspace; unitcell=(N1, N2))
+# normalize vertex tensors
+for ind in CartesianIndices(peps.vertices)
+    peps.vertices[ind] /= norm(peps.vertices[ind], Inf)
+end
+
+# simple update
+dts = [1e-2, 1e-3, 4e-4]
+tols = [1e-6, 1e-8, 1e-8]
+maxiter = 10000
+for (dt, tol) in zip(dts, tols)
+    alg = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
+    result = simpleupdate(peps, ham, alg; bipartite=true)
+    global peps = result[1]
+end
+# measure physical quantities with CTMRG
+peps_ = InfinitePEPS(peps)
+envs = CTMRGEnv(rand, Float64, peps_, Espace)
+ctm_alg = SequentialCTMRG(; tol=1e-10, verbosity=2, trscheme=trscheme_envs)
+envs = leading_boundary(envs, peps_, ctm_alg)
+meas = measure_heis(peps_, ham, envs)
+display(meas)
+@info @sprintf("Energy = %.8f\n", meas["e_site"])
+@info @sprintf("Staggered magnetization = %.8f\n", mean(meas["mag_norm"]))
+@test isapprox(meas["e_site"], -0.6675; atol=1e-3)
+@test isapprox(mean(meas["mag_norm"]), 0.3767; atol=1e-3)

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")
@@ -21,6 +21,7 @@ include("utility/mirror.jl")
 include("utility/diffset.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
+include("utility/fullenv_truncation.jl")
 
 include("networks/tensors.jl")
 include("networks/infinitesquarenetwork.jl")
@@ -50,8 +51,9 @@ include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 
-include("algorithms/time_evolution/gatetools.jl")
+include("algorithms/time_evolution/evoltools.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
+include("algorithms/time_evolution/ntupdate.jl")
 
 include("algorithms/toolbox.jl")
 
@@ -194,7 +196,10 @@ export PEPSOptimize, GeomSum, ManualIter, LinSolver
 export fixedpoint
 
 export absorb_weight
+export ALSTruncation, FullEnvTruncation, fullenv_truncate
 export su_iter, simpleupdate, SimpleUpdate
+export NTUEnvNN, NTUEnvNNN, NTUEnvNNNp
+export ntu_iter, ntupdate, NTUpdate
 
 export InfinitePartitionFunction
 export InfinitePEPS, InfiniteTransferPEPS

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/time_evolution/bondenv/env_approx.jl

@@ -0,0 +1,11 @@
+"""
+Use CTMRG to approximately contract 
+a cluster around a bond to construct bond environment
+"""
+@kwdef struct CTMBondEnv <: BondEnvAlgorithm
+    pad::Int = 1
+end
+
+function bondenv_ctm(row::Int, col::Int, peps::InfinitePEPS, alg::CTMBondEnv)
+    return error("Not implemented")
+end