Add simple update algorithm (#97)

* remove CTMRG redundant logging * add simple update algorithm * add full update core algorithm * improve formatting * add TODO for svd initialization of ALS optimization * remove a redundant SVD in full update * prepare for addition of full-infinite env CTMRG * Add `length` for SUWeight * Define `Base.show` for `SUWeight` * add test for simple update * add test for full update * update formatting * Refactor calc_convergence Co-authored-by: Lukas Devos <[email protected]> * Update sdiag_pow for latest TensorKit Co-authored-by: Lukas Devos <[email protected]> * Rename folder "timeevol" to "time_evolution" * Replace `maxabs` by infinity norm * Focus on simple update * implement `InfiniteWeightPEPS` * Update to the Heisenberg tensors Using MPSKitModels, the functions gen_siteop and gen_bondop are not necessary anymore. The relevant tensors and their tensor product can be defined directly from MPSKitModels * Minor fix after using operators from MPSKitModels * Use Julia's logging system; shorten SU test * Solve deprecation warning for `permute` * Remove buggy in-place rotations and reflections * Minor refactoring * Print more message during simple update * Replace custom rho with existing exp. value calculation * Integrate SU with LocalOperator * remove accidentally added test * Add rotation and reflection of LocalOperator * Improve `get_gateterm` * fix format * fix format again * Introduce `SimpleUpdate` algorithm struct * Export `simpleupdate`; remove abbreviated `SU` * add FixedSpaceTruncation for simple update * Create heisenberg_sufu_onesite This tests whether one-body terms can be accurately handled by SU by rewriting them as one-body terms. This is just an example and should probably not be in the final test set. * format fix * Add spin U(1) symmetry to Heisenberg model SU test * Add SU-AD test for heisenberg * fix formatting * add check that all operators are twosite * update check two-body terms the code now also checks whether all interactions are defined on nearest neighbours * Change PEPSOptimize parameters in `suad` test * Hubbard model added the Hubbard model as in the pull request from lkdvos added test for simple update on the Hubbard model for t = 1, U = 6, half filling * Clean up test on Hubbard model * Add t-J model Hamiltonian * Squashed commit of the following: commit e2545a3 Merge: 36973e7 9532507 Author: Paul Brehmer <[email protected]> Date: Thu Dec 5 11:06:04 2024 +0100 Merge pull request #90 from QuantumKitHub/pb-improve-sequential Make `:sequential` act column-wise commit 9532507 Merge: 77fc207 36973e7 Author: Lukas Devos <[email protected]> Date: Wed Dec 4 15:32:45 2024 -0500 Merge branch 'master' into pb-improve-sequential commit 77fc207 Author: Lukas Devos <[email protected]> Date: Tue Dec 3 15:33:15 2024 -0500 reenable expansion for simultaneous ctmrg commit 51f2cc4 Author: Lukas Devos <[email protected]> Date: Tue Dec 3 14:18:36 2024 -0500 excise expansion step commit 3754b67 Merge: 8a09a82 08cf1dc Author: Paul Brehmer <[email protected]> Date: Fri Nov 8 14:57:56 2024 +0100 Merge branch 'master' into pb-improve-sequential * Update `ctmrg_leftmove` with latest upstream * Apply suggestions from code review Co-authored-by: Lukas Devos <[email protected]> * Rename for `absorb_weight` * Improve SUWeight construction * Move geometric operations on LocalOperator * Remove `sdiag_inv_sqrt` (replaced with `sdiag_pow`) * Fix rrule for sdiag_pow and formatting * Fix Heisenberg SU test * Improve some docstring * Add back `length` for `SUWeight` * update actions This should allow the tests to run, because the secrets are now explicitly passed on * Refactoring `SUWeight` * Update simple update test * disable multithreading * remove superfluous broadcats * Add function signatures * Merge Heisenberg tests * Fix docstring and error messages * Promote `SUWeight` to a custom `struct` * Fix formatting * Fix pow rrule and make eltype more consistent * Fix `similar` for `InfinitePEPO` after modifying `eltype` * Remove `eltype` of `InfiniteWeightPEPS` * Scrap eltype for CTMRGEnv, add args to similar(::PEPO) * Fix conj in sdiag_pow rrule * Cicrumvent Heisenberg tests errors by using GMRES to differentiate SVD * small improvements * update MPSKitModels compat --------- Co-authored-by: Lukas Devos <[email protected]> Co-authored-by: sanderdemeyer <[email protected]> Co-authored-by: Paul Brehmer <[email protected]>
QuantumKitHub · Dec 11, 2024 · 45d00e8 · 45d00e8
1 parent e2545a3
commit 45d00e8
Show file tree

Hide file tree

Showing 23 changed files with 955 additions and 83 deletions.
diff --git a/.github/workflows/CI.yml → .github/workflows/Tests.yml b/.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 }}
+
diff --git a/Project.toml b/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"

diff --git a/examples/hubbard_su.jl b/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)
diff --git a/src/PEPSKit.jl b/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
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
@@ -133,6 +133,26 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG)
     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 @@ function ctmrg_iter(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 = 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 @@ function build_projectors(
     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 @@ function build_projectors(
     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

diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl
@@ -173,6 +173,17 @@ function calc_convergence(envs, CSold, TSold)
     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...)
 
 """

diff --git a/src/algorithms/peps_opt.jl b/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
 

diff --git a/src/algorithms/time_evolution/gatetools.jl b/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