Merge branch 'master' of github.com:quantumghent/PEPSKit.jl into pb-f…

…ull-inf-proj

.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,59 @@
+using Test
+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, 6
+ham = hubbard_model(Float64, Trivial, Trivial, InfiniteSquare(N1, N2); 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)
+    alg = SimpleUpdate(dt, tol, maxiter, trscheme)
+    peps, = simpleupdate(peps, ham, alg; bipartite=false)
+end
+
+# absorb 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]
+    ctm_alg = SequentialCTMRG(; maxiter=300, tol=1e-7)
+    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.
+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)
+@info "Energy           = $E"
+@info "Benchmark energy = $E_exact"
+@test isapprox(E, E_exact; atol=5e-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")
@@ -46,6 +48,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/simpleupdate.jl")
+
 include("algorithms/toolbox.jl")
 
 include("algorithms/peps_opt.jl")
@@ -186,13 +191,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/contractions/ctmrg_contractions.jl

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

src/algorithms/ctmrg/ctmrg.jl

@@ -8,7 +8,7 @@ for contracting infinite PEPS.
 abstract type CTMRGAlgorithm end
 
 """
-    ctmrg_iteration(state, env, alg::CTMRG) -> env′, info
+    ctmrg_iteration(state, env, alg::CTMRGAlgorithm) -> env′, info
 
 Perform a single CTMRG iteration in which all directions are being grown and renormalized.
 """

src/algorithms/ctmrg/sparse_environments.jl

@@ -104,7 +104,7 @@ end
 
 # 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_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/peps_opt.jl

@@ -82,7 +82,7 @@ end
 
 Algorithm struct that represent PEPS ground-state optimization using AD.
 Set the algorithm to contract the infinite PEPS in `boundary_alg`;
-currently only `CTMRG` is supported. The `optimizer` computes the gradient directions
+currently only `CTMRGAlgorithm`s are supported. The `optimizer` computes the gradient directions
 based on the CTMRG gradient and updates the PEPS parameters. In this optimization,
 the CTMRG runs can be started on the converged environments of the previous optimizer
 step by setting `reuse_env` to true. Otherwise a random environment is used at each
@@ -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