Add simple update algorithm

.github/workflows/CI.yml → .github/workflows/Tests.yml

@@ -1,4 +1,4 @@
-name: CI
+name: Tests
 on:
   push:
     branches:
@@ -23,8 +23,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'lts'
-          - '1'
+          - 'lts' # minimal supported version
+          - '1' # latest released Julia version
         group:
           - ctmrg
           - boundarymps
@@ -34,9 +34,11 @@ jobs:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
-    uses: "QuantumKitHub/.github/.github/workflows/tests.yml@main"
+    uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/Tests.yml@main"
     with:
       group: "${{ matrix.group }}"
       julia-version: "${{ matrix.version }}"
       os: "${{ matrix.os }}"
-    secrets: inherit
+    secrets:
+      CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+

Project.toml

@@ -32,7 +32,7 @@ KrylovKit = "0.8"
 LinearAlgebra = "1"
 LoggingExtras = "1"
 MPSKit = "0.11"
-MPSKitModels = "0.3"
+MPSKitModels = "0.3.5"
 OhMyThreads = "0.7"
 OptimKit = "0.3"
 Printf = "1"

examples/hubbard_su.jl

@@ -0,0 +1,62 @@
+using Test
+using Printf
+using Random
+using PEPSKit
+using TensorKit
+
+# random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
+Dcut, symm = 8, Trivial
+N1, N2 = 2, 2
+Random.seed!(10)
+if symm == Trivial
+    Pspace = Vect[fℤ₂](0 => 2, 1 => 2)
+    Vspace = Vect[fℤ₂](0 => Dcut / 2, 1 => Dcut / 2)
+else
+    error("Not implemented")
+end
+
+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
+# Hubbard model Hamiltonian at half-filling
+t, U = 1.0, 6.0
+ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(N1, N2); t=t, U=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)
+    alg = SimpleUpdate(dt, tol, maxiter, trscheme)
+    result = simpleupdate(peps, ham, alg; bipartite=false)
+    global peps = result[1]
+end
+
+# absort weight into site tensors
+peps = InfinitePEPS(peps)
+# CTMRG
+χenv0, χenv = 6, 20
+Espace = Vect[fℤ₂](0 => χenv0 / 2, 1 => χenv0 / 2)
+envs = CTMRGEnv(randn, Float64, peps, Espace)
+for χ in [χenv0, χenv]
+    trscheme = truncerr(1e-9) & truncdim(χ)
+    ctm_alg = CTMRG(;
+        maxiter=40, tol=1e-10, verbosity=3, trscheme=trscheme, ctmrgscheme=:sequential
+    )
+    global envs = leading_boundary(envs, peps, ctm_alg)
+end
+
+"""
+Benchmark values of the ground state energy from
+Qin, M., Shi, H., & Zhang, S. (2016). Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method. Physical Review B, 94(8), 085103.
+"""
+# measure physical quantities
+E = costfun(peps, envs, ham) / (N1 * N2)
+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
+@info @sprintf("Energy           = %.8f\n", E)
+@info @sprintf("Benchmark energy = %.8f\n", E_exact)
+@test isapprox(E, E_exact; atol=1e-2)

src/PEPSKit.jl

@@ -17,12 +17,14 @@ include("utility/util.jl")
 include("utility/diffable_threads.jl")
 include("utility/svd.jl")
 include("utility/rotations.jl")
+include("utility/mirror.jl")
 include("utility/diffset.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
 
 include("states/abstractpeps.jl")
 include("states/infinitepeps.jl")
+include("states/infiniteweightpeps.jl")
 
 include("operators/transferpeps.jl")
 include("operators/infinitepepo.jl")
@@ -43,6 +45,9 @@ include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 
+include("algorithms/time_evolution/gatetools.jl")
+include("algorithms/time_evolution/simpleupdate.jl")
+
 include("algorithms/toolbox.jl")
 
 include("algorithms/peps_opt.jl")
@@ -59,7 +64,7 @@ include("utility/symmetrization.jl")
         const ctmrgscheme = :simultaneous
         const reuse_env = true
         const trscheme = FixedSpaceTruncation()
-        const fwd_alg = TensorKit.SVD()
+        const fwd_alg = TensorKit.SDD()
         const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
         const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
         const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
@@ -167,13 +172,18 @@ export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver
 export fixedpoint
 
+export absorb_weight
+export su_iter, simpleupdate, SimpleUpdate
+
 export InfinitePEPS, InfiniteTransferPEPS
+export SUWeight, InfiniteWeightPEPS
 export InfinitePEPO, InfiniteTransferPEPO
 export initializeMPS, initializePEPS
 export ReflectDepth, ReflectWidth, Rotate, RotateReflect
 export symmetrize!, symmetrize_retract_and_finalize!
 export showtypeofgrad
 export InfiniteSquare, vertices, nearest_neighbours, next_nearest_neighbours
-export transverse_field_ising, heisenberg_XYZ, j1_j2, pwave_superconductor
+export transverse_field_ising, heisenberg_XYZ, j1_j2
+export pwave_superconductor, hubbard_model, tj_model
 
 end # module

src/algorithms/ctmrg/ctmrg.jl

@@ -133,6 +133,26 @@
     end
 end
 
+"""
+    ctmrg_leftmove(col::Int, state, envs::CTMRGEnv, alg::SequentialCTMRG)
+
+Perform 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 = ctmrg_projectors(col, state, envs, alg)
+    envs = ctmrg_renormalize(col, projectors, state, envs, alg)
+    return envs, info
+end
+
 """
     ctmrg_iter(state, envs::CTMRGEnv, alg::CTMRG) -> envs′, info
 
@@ -143,14 +163,12 @@
     ϵ = zero(real(scalartype(state)))
     for _ in 1:4 # rotate
         for col in 1:size(state, 2) # left move column-wise
-            projectors, info = ctmrg_projectors(col, state, envs, alg)
-            envs = ctmrg_renormalize(col, projectors, state, envs, alg)
+            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
 function ctmrg_iter(state, envs::CTMRGEnv, alg::SimultaneousCTMRG)
@@ -300,7 +318,7 @@
     Q::AbstractTensorMap{E,3,3},
     Q_next::AbstractTensorMap{E,3,3},
 ) where {E<:ElementarySpace}
-    isqS = sdiag_inv_sqrt(S)
+    isqS = sdiag_pow(S, -0.5)
     P_left = Q_next * V' * isqS
     P_right = isqS * U' * Q
     return P_left, P_right
@@ -312,7 +330,7 @@
     Q::EnlargedCorner,
     Q_next::EnlargedCorner,
 ) where {E<:ElementarySpace}
-    isqS = sdiag_inv_sqrt(S)
+    isqS = sdiag_pow(S, -0.5)
     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/ctmrg/gaugefix.jl

@@ -173,6 +173,17 @@
     return max(ΔCS, ΔTS), CSnew, TSnew
 end
 
+"""
+    calc_convergence(envsNew::CTMRGEnv, envsOld::CTMRGEnv)
+
+Calculate convergence of CTMRG by comparing the singular values of CTM tensors.
+"""
+function calc_convergence(envsNew::CTMRGEnv, envsOld::CTMRGEnv)
+    CSOld = map(x -> tsvd(x)[2], envsOld.corners)
+    TSOld = map(x -> tsvd(x)[2], envsOld.edges)
+    return calc_convergence(envsNew, CSOld, TSOld)
+end
+
 @non_differentiable calc_convergence(args...)
 
 """

src/algorithms/peps_opt.jl

@@ -159,8 +159,8 @@ function fixedpoint(
 
     if scalartype(env₀) <: Real
         env₀ = complex(env₀)
-        @warn "the provided real environment was converted to a complex environment since\
-        :fixed mode generally produces complex gauges; use :diffgauge mode instead to work\
+        @warn "the provided real environment was converted to a complex environment since \
+        :fixed mode generally produces complex gauges; use :diffgauge mode instead to work \
         with purely real environments"
     end
 

src/algorithms/time_evolution/gatetools.jl

@@ -0,0 +1,52 @@
+"""
+    get_gate(dt::Float64, H::LocalOperator)
+
+Compute `exp(-dt * H)` from the nearest neighbor Hamiltonian `H`.
+"""
+function get_gate(dt::Float64, H::LocalOperator)
+    @assert all([
+        length(op.dom) == 2 && norm(Tuple(terms[2] - terms[1])) == 1.0 for
+        (terms, op) in H.terms
+    ]) "Only nearest-neighbour terms allowed"
+    return LocalOperator(
+        H.lattice, Tuple(sites => exp(-dt * op) for (sites, op) in H.terms)...
+    )
+end
+
+"""
+    is_equivalent(bond1::NTuple{2,CartesianIndex{2}}, bond2::NTuple{2,CartesianIndex{2}}, (Nrow, Ncol)::NTuple{2,Int})
+
+Check if two 2-site bonds are related by a (periodic) lattice translation.
+"""
+function is_equivalent(
+    bond1::NTuple{2,CartesianIndex{2}},
+    bond2::NTuple{2,CartesianIndex{2}},
+    (Nrow, Ncol)::NTuple{2,Int},
+)
+    r1 = bond1[1] - bond1[2]
+    r2 = bond2[1] - bond2[2]
+    shift_row = bond1[1][1] - bond2[1][1]
+    shift_col = bond1[1][2] - bond2[1][2]
+    return r1 == r2 && mod(shift_row, Nrow) == 0 && mod(shift_col, Ncol) == 0
+end
+
+"""
+    get_gateterm(gate::LocalOperator, bond::NTuple{2,CartesianIndex{2}})
+
+Get the term of a 2-site gate acting on a certain bond.
+Input `gate` should only include one term for each nearest neighbor bond.
+"""
+function get_gateterm(gate::LocalOperator, bond::NTuple{2,CartesianIndex{2}})
+    label = findall(p -> is_equivalent(p.first, bond, size(gate.lattice)), gate.terms)
+    if length(label) == 0
+        # try reversed site order
+        label = findall(
+            p -> is_equivalent(p.first, reverse(bond), size(gate.lattice)), gate.terms
+        )
+        @assert length(label) == 1
+        return permute(gate.terms[label[1]].second, ((2, 1), (4, 3)))
+    else
+        @assert length(label) == 1
+        return gate.terms[label[1]].second
+    end
+end