From 45d00e84f39f8dfaec4ed1876545eec9cc862086 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 12 Dec 2024 01:53:57 +0800 Subject: [PATCH] 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 * Update sdiag_pow for latest TensorKit Co-authored-by: Lukas Devos * 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 e2545a379f0efb0e70087bb2a526ce08a6986bcc Merge: 36973e7 9532507 Author: Paul Brehmer Date: Thu Dec 5 11:06:04 2024 +0100 Merge pull request #90 from QuantumKitHub/pb-improve-sequential Make `:sequential` act column-wise commit 9532507408401bb09916ff38834c19a085ff1639 Merge: 77fc207 36973e7 Author: Lukas Devos Date: Wed Dec 4 15:32:45 2024 -0500 Merge branch 'master' into pb-improve-sequential commit 77fc207184198bdc8a1aa9d7d9b0acf78129d576 Author: Lukas Devos Date: Tue Dec 3 15:33:15 2024 -0500 reenable expansion for simultaneous ctmrg commit 51f2cc4baefe57b8351e5ab3f6835f37726a5db3 Author: Lukas Devos Date: Tue Dec 3 14:18:36 2024 -0500 excise expansion step commit 3754b67cfe249ffb0aeff33cebd0d22ac322a9c9 Merge: 8a09a82 08cf1dc Author: Paul Brehmer 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 * 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 Co-authored-by: sanderdemeyer <80397440+Sander-De-Meyer@users.noreply.github.com> Co-authored-by: Paul Brehmer --- .github/workflows/{CI.yml => Tests.yml} | 12 +- Project.toml | 2 +- examples/hubbard_su.jl | 62 +++++ src/PEPSKit.jl | 14 +- src/algorithms/ctmrg/ctmrg.jl | 28 ++- src/algorithms/ctmrg/gaugefix.jl | 11 + src/algorithms/peps_opt.jl | 4 +- src/algorithms/time_evolution/gatetools.jl | 52 ++++ src/algorithms/time_evolution/simpleupdate.jl | 218 +++++++++++++++++ src/algorithms/toolbox.jl | 4 +- src/environments/ctmrg_environments.jl | 44 +++- src/environments/transferpeps_environments.jl | 3 +- src/operators/infinitepepo.jl | 7 +- src/operators/lattices/squarelattice.jl | 2 + src/operators/localoperator.jl | 100 ++++++++ src/operators/models.jl | 59 ++++- src/states/infinitepeps.jl | 7 +- src/states/infiniteweightpeps.jl | 223 ++++++++++++++++++ src/utility/mirror.jl | 12 + src/utility/svd.jl | 2 +- src/utility/util.jl | 56 +++-- test/heisenberg.jl | 114 +++++++-- test/runtests.jl | 2 +- 23 files changed, 955 insertions(+), 83 deletions(-) rename .github/workflows/{CI.yml => Tests.yml} (75%) create mode 100644 examples/hubbard_su.jl create mode 100644 src/algorithms/time_evolution/gatetools.jl create mode 100644 src/algorithms/time_evolution/simpleupdate.jl create mode 100644 src/states/infiniteweightpeps.jl create mode 100644 src/utility/mirror.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/Tests.yml similarity index 75% rename from .github/workflows/CI.yml rename to .github/workflows/Tests.yml index 85c75900..1f19b865 100644 --- a/.github/workflows/CI.yml +++ b/.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 \ No newline at end of file + secrets: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + diff --git a/Project.toml b/Project.toml index cffc9c1a..ba49bc98 100644 --- 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 new file mode 100644 index 00000000..4bd84846 --- /dev/null +++ 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 index c8a6b032..a69e5331 100644 --- 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 index 91d86159..317a7231 100644 --- 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 index 5e63d5ef..338395b2 100644 --- 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 index 39e138a1..8b927f93 100644 --- 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 new file mode 100644 index 00000000..91ebaa0f --- /dev/null +++ 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 diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl new file mode 100644 index 00000000..685d069b --- /dev/null +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -0,0 +1,218 @@ +""" + struct SimpleUpdate + +Algorithm struct for simple update (SU) of infinite PEPS with bond weights. +Each SU run is converged when the singular value difference becomes smaller than `tol`. +""" +struct SimpleUpdate + dt::Float64 + tol::Float64 + maxiter::Int + trscheme::TensorKit.TruncationScheme +end + +function truncation_scheme(alg::SimpleUpdate, v::ElementarySpace) + if alg.trscheme isa FixedSpaceTruncation + return truncspace(v) + else + return alg.trscheme + end +end + +""" +_su_bondx!(row::Int, col::Int, gate::AbstractTensorMap{S,2,2}, + peps::InfiniteWeightPEPS, alg::SimpleUpdate) where {S<:ElementarySpace} + +Simple update of the x-bond `peps.weights[1,r,c]`. + +``` + [2,r,c] [2,r,c+1] + ↓ ↓ + [1,r,c-1] ← T[r,c] ← [1,r,c] ←- T[r,c+1] ← [1,r,c+1] + ↓ ↓ + [2,r+1,c] [2,r+1,c+1] +``` +""" +function _su_bondx!( + row::Int, + col::Int, + gate::AbstractTensorMap{S,2,2}, + peps::InfiniteWeightPEPS, + alg::SimpleUpdate, +) where {S<:ElementarySpace} + Nr, Nc = size(peps) + @assert 1 <= row <= Nr && 1 <= col <= Nc + row2, col2 = row, _next(col, Nc) + T1, T2 = peps.vertices[row, col], peps.vertices[row2, col2] + # absorb environment weights + for ax in (2, 4, 5) + T1 = absorb_weight(T1, row, col, ax, peps.weights) + end + for ax in (2, 3, 4) + T2 = absorb_weight(T2, row2, col2, ax, peps.weights) + end + # absorb bond weight + T1 = absorb_weight(T1, row, col, 3, peps.weights; sqrtwt=true) + T2 = absorb_weight(T2, row2, col2, 5, peps.weights; sqrtwt=true) + #= QR and LQ decomposition + + 2 1 1 2 + ↓ ↗ ↓ ↗ + 5 ← T ← 3 ====> 3 ← X ← 4 ← 1 ← aR ← 3 + ↓ ↓ + 4 2 + + 2 1 2 2 + ↓ ↗ ↗ ↓ + 5 ← T ← 3 ====> 1 ← bL ← 3 ← 1 ← Y ← 3 + ↓ ↓ + 4 4 + =# + X, aR = leftorth(T1, ((2, 4, 5), (1, 3)); alg=QRpos()) + bL, Y = rightorth(T2, ((5, 1), (2, 3, 4)); alg=LQpos()) + #= apply gate + + -2 -3 + ↑ ↑ + |----gate---| + ↑ ↑ + 1 2 + ↑ ↑ + -1← aR -← 3 -← bL ← -4 + =# + @tensor tmp[-1 -2; -3 -4] := gate[-2, -3, 1, 2] * aR[-1, 1, 3] * bL[3, 2, -4] + # SVD + aR, s, bL, ϵ = tsvd!(tmp; trunc=truncation_scheme(alg, space(T1, 3))) + #= + -2 -1 -1 -2 + | ↗ ↗ | + -5- X ← 1 ← aR - -3 -5 - bL ← 1 ← Y - -3 + | | + -4 -4 + =# + @tensor T1[-1; -2 -3 -4 -5] := X[-2, -4, -5, 1] * aR[1, -1, -3] + @tensor T2[-1; -2 -3 -4 -5] := bL[-5, -1, 1] * Y[1, -2, -3, -4] + # remove environment weights + for ax in (2, 4, 5) + T1 = absorb_weight(T1, row, col, ax, peps.weights; invwt=true) + end + for ax in (2, 3, 4) + T2 = absorb_weight(T2, row2, col2, ax, peps.weights; invwt=true) + end + # update tensor dict and weight on current bond + # (max element of weight is normalized to 1) + peps.vertices[row, col], peps.vertices[row2, col2] = T1, T2 + peps.weights[1, row, col] = s / norm(s, Inf) + return ϵ +end + +""" + su_iter(gate::LocalOperator, peps::InfiniteWeightPEPS, alg::SimpleUpdate; bipartite::Bool=false) + +One round of simple update on `peps` applying the nearest neighbor `gate`. +""" +function su_iter( + gate::LocalOperator, peps::InfiniteWeightPEPS, alg::SimpleUpdate; bipartite::Bool=false +) + @assert size(gate.lattice) == size(peps) + Nr, Nc = size(peps) + if bipartite + @assert Nr == Nc == 2 + end + # TODO: make algorithm independent on the choice of dual in the network + for (r, c) in Iterators.product(1:Nr, 1:Nc) + @assert [isdual(space(peps.vertices[r, c], ax)) for ax in 1:5] == [0, 1, 1, 0, 0] + @assert [isdual(space(peps.weights[1, r, c], ax)) for ax in 1:2] == [0, 1] + @assert [isdual(space(peps.weights[2, r, c], ax)) for ax in 1:2] == [0, 1] + end + peps2 = deepcopy(peps) + gate_mirrored = mirror_antidiag(gate) + for direction in 1:2 + # mirror the y-weights to x-direction + # to update them using code for x-weights + if direction == 2 + peps2 = mirror_antidiag(peps2) + end + if bipartite + for r in 1:2 + rp1 = _next(r, 2) + term = get_gateterm( + direction == 1 ? gate : gate_mirrored, + (CartesianIndex(r, 1), CartesianIndex(r, 2)), + ) + ϵ = _su_bondx!(r, 1, term, peps2, alg) + peps2.vertices[rp1, 2] = deepcopy(peps2.vertices[r, 1]) + peps2.vertices[rp1, 1] = deepcopy(peps2.vertices[r, 2]) + peps2.weights[1, rp1, 2] = deepcopy(peps2.weights[1, r, 1]) + end + else + for site in CartesianIndices(peps2.vertices) + r, c = Tuple(site) + term = get_gateterm( + direction == 1 ? gate : gate_mirrored, + (CartesianIndex(r, c), CartesianIndex(r, c + 1)), + ) + ϵ = _su_bondx!(r, c, term, peps2, alg) + end + end + if direction == 2 + peps2 = mirror_antidiag(peps2) + end + end + return peps2 +end + +""" + simpleupdate(peps::InfiniteWeightPEPS, ham::LocalOperator, alg::SimpleUpdate; + bipartite::Bool=false, check_int::Int=500) + +Perform simple update with nearest neighbor Hamiltonian `ham`, where the evolution +information is printed every `check_int` steps. + +If `bipartite == true` (for square lattice), a unit cell size of `(2, 2)` is assumed, +as well as tensors and x/y weights which are the same across the diagonals, i.e. at +`(row, col)` and `(row+1, col+1)`. +""" +function simpleupdate( + peps::InfiniteWeightPEPS, + ham::LocalOperator, + alg::SimpleUpdate; + bipartite::Bool=false, + check_int::Int=500, +) + time_start = time() + N1, N2 = size(peps) + if bipartite + @assert N1 == N2 == 2 + end + # exponentiating the 2-site Hamiltonian gate + gate = get_gate(alg.dt, ham) + wtdiff = 1.0 + wts0 = deepcopy(peps.weights) + for count in 1:(alg.maxiter) + time0 = time() + peps = su_iter(gate, peps, alg; bipartite=bipartite) + wtdiff = compare_weights(peps.weights, wts0) + converge = (wtdiff < alg.tol) + cancel = (count == alg.maxiter) + wts0 = deepcopy(peps.weights) + time1 = time() + if ((count == 1) || (count % check_int == 0) || converge || cancel) + @info "Space of x-weight at [1, 1] = $(space(peps.weights[1, 1, 1], 1))" + label = (converge ? "conv" : (cancel ? "cancel" : "iter")) + message = @sprintf( + "SU %s %-7d: dt = %.0e, weight diff = %.3e, time = %.3f sec\n", + label, + count, + alg.dt, + wtdiff, + time1 - ((converge || cancel) ? time_start : time0) + ) + cancel ? (@warn message) : (@info message) + end + if converge || cancel + break + end + end + return peps, wtdiff +end diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 13e725e3..b89a072a 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -58,11 +58,11 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) end """ - correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2) + correlation_length(peps::InfinitePEPS, env::CTMRGEnv; num_vals=2) Compute the PEPS correlation length based on the horizontal and vertical transfer matrices. Additionally the (normalized) eigenvalue spectrum is -returned. Specify the number of computed eigenvalues with `howmany`. +returned. Specify the number of computed eigenvalues with `num_vals`. """ function MPSKit.correlation_length(peps::InfinitePEPS, env::CTMRGEnv; num_vals=2) T = scalartype(peps) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index c026762b..7bca4b90 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -160,7 +160,7 @@ end """ CTMRGEnv( - [f=randn, ComplexF64], D_north::S, D_south::S, chi_north::S, [chi_east::S], [chi_south::S], [chi_west::S]; unitcell::Tuple{Int,Int}=(1, 1), + [f=randn, ComplexF64], D_north::S, D_east::S, chi_north::S, [chi_east::S], [chi_south::S], [chi_west::S]; unitcell::Tuple{Int,Int}=(1, 1), ) where {S<:Union{Int,ElementarySpace}} Construct a CTMRG environment by specifying the north and east virtual spaces of the @@ -173,7 +173,7 @@ corresponding edge tensor for each direction. """ function CTMRGEnv( D_north::S, - D_south::S, + D_east::S, chi_north::S, chi_east::S=chi_north, chi_south::S=chi_north, @@ -184,7 +184,7 @@ function CTMRGEnv( randn, ComplexF64, fill(D_north, unitcell), - fill(D_south, unitcell), + fill(D_east, unitcell), fill(chi_north, unitcell), fill(chi_east, unitcell), fill(chi_south, unitcell), @@ -195,7 +195,7 @@ function CTMRGEnv( f, T, D_north::S, - D_south::S, + D_east::S, chi_north::S, chi_east::S=chi_north, chi_south::S=chi_north, @@ -206,7 +206,7 @@ function CTMRGEnv( f, T, fill(D_north, unitcell), - fill(D_south, unitcell), + fill(D_east, unitcell), fill(chi_north, unitcell), fill(chi_east, unitcell), fill(chi_south, unitcell), @@ -362,16 +362,44 @@ function Base.rotl90(env::CTMRGEnv{C,T}) where {C,T} Array{C,3}(undef, 4, size(env.corners, 3), size(env.corners, 2)) ) edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env.edges, 3), size(env.edges, 2))) + for dir in 1:4 + dir2 = _prev(dir, 4) + corners′[dir2, :, :] = rotl90(env.corners[dir, :, :]) + edges′[dir2, :, :] = rotl90(env.edges[dir, :, :]) + end + return CTMRGEnv(copy(corners′), copy(edges′)) +end +# Rotate corners & edges clockwise +function Base.rotr90(env::CTMRGEnv{C,T}) where {C,T} + # Initialize rotated corners & edges with rotated sizes + corners′ = Zygote.Buffer( + Array{C,3}(undef, 4, size(env.corners, 3), size(env.corners, 2)) + ) + edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env.edges, 3), size(env.edges, 2))) for dir in 1:4 - corners′[_prev(dir, 4), :, :] = rotl90(env.corners[dir, :, :]) - edges′[_prev(dir, 4), :, :] = rotl90(env.edges[dir, :, :]) + dir2 = _next(dir, 4) + corners′[dir2, :, :] = rotr90(env.corners[dir, :, :]) + edges′[dir2, :, :] = rotr90(env.edges[dir, :, :]) end + return CTMRGEnv(copy(corners′), copy(edges′)) +end +# Rotate corners & edges by 180 degrees +function Base.rot180(env::CTMRGEnv{C,T}) where {C,T} + # Initialize rotated corners & edges with rotated sizes + corners′ = Zygote.Buffer( + Array{C,3}(undef, 4, size(env.corners, 2), size(env.corners, 3)) + ) + edges′ = Zygote.Buffer(Array{T,3}(undef, 4, size(env.edges, 2), size(env.edges, 3))) + for dir in 1:4 + dir2 = _next(_next(dir, 4), 4) + corners′[dir2, :, :] = rot180(env.corners[dir, :, :]) + edges′[dir2, :, :] = rot180(env.edges[dir, :, :]) + end return CTMRGEnv(copy(corners′), copy(edges′)) end -Base.eltype(env::CTMRGEnv) = eltype(env.corners[1]) Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) function eachcoordinate(x::CTMRGEnv) return collect(Iterators.product(axes(x, 2), axes(x, 3))) diff --git a/src/environments/transferpeps_environments.jl b/src/environments/transferpeps_environments.jl index 4332c232..316eaf77 100644 --- a/src/environments/transferpeps_environments.jl +++ b/src/environments/transferpeps_environments.jl @@ -120,8 +120,7 @@ function MPSKit.transfer_spectrum( @assert size(below) == size(O) numrows = size(above, 1) - envtype = eltype(init[1]) - eigenvals = Vector{Vector{scalartype(envtype)}}(undef, numrows) + eigenvals = Vector{Vector{scalartype(init[1])}}(undef, numrows) @threads for cr in 1:numrows L0, = init[cr] diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index e136aa70..6fad9fd3 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -109,12 +109,13 @@ end Base.size(T::InfinitePEPO) = size(T.A) Base.size(T::InfinitePEPO, i) = size(T.A, i) Base.length(T::InfinitePEPO) = length(T.A) -Base.eltype(T::InfinitePEPO) = eltype(T.A) -VectorInterface.scalartype(T::InfinitePEPO) = scalartype(T.A) +Base.eltype(T::InfinitePEPO) = eltype(typeof(T)) +Base.eltype(::Type{<:InfinitePEPO{T}}) where {T} = T +VectorInterface.scalartype(::Type{T}) where {T<:InfinitePEPO} = scalartype(eltype(T)) ## Copy Base.copy(T::InfinitePEPO) = InfinitePEPO(copy(T.A)) -Base.similar(T::InfinitePEPO) = InfinitePEPO(similar(T.A)) +Base.similar(T::InfinitePEPO, args...) = InfinitePEPO(similar(T.A, args...)) Base.repeat(T::InfinitePEPO, counts...) = InfinitePEPO(repeat(T.A, counts...)) Base.getindex(T::InfinitePEPO, args...) = Base.getindex(T.A, args...) diff --git a/src/operators/lattices/squarelattice.jl b/src/operators/lattices/squarelattice.jl index 097a4f7b..a2f3d9b5 100644 --- a/src/operators/lattices/squarelattice.jl +++ b/src/operators/lattices/squarelattice.jl @@ -12,6 +12,8 @@ struct InfiniteSquare <: AbstractLattice{2} end end +Base.size(lattice::InfiniteSquare) = (lattice.Nrows, lattice.Ncols) + function vertices(lattice::InfiniteSquare) return CartesianIndices((1:(lattice.Nrows), 1:(lattice.Ncols))) end diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 2698ce8c..90845de0 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -108,3 +108,103 @@ end Base.:-(O::LocalOperator) = -1 * O Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2) + +# Rotation and mirroring +# ---------------------- + +""" + _mirror_antidiag_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} + ) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + +Get the position of `site` after reflection about the anti-diagonal line. +""" +function _mirror_antidiag_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} +) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + r, c = site[1], site[2] + return CartesianIndex(1 - c + Ncol, 1 - r + Nrow) +end + +""" + _rotr90_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} + ) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + +Get the position of `site` after clockwise (right) rotation by 90 degrees. +""" +function _rotr90_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} +) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + r, c = site[1], site[2] + return CartesianIndex(c, 1 + Nrow - r) +end + +""" + _rotl90_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} + ) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + +Get the position of `site` after counter-clockwise (left) rotation by 90 degrees. +""" +function _rotl90_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} +) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + r, c = site[1], site[2] + return CartesianIndex(1 + Ncol - c, r) +end + +""" + _rot180_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} + ) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + +Get the position of `site` after rotation by 180 degrees. +""" +function _rot180_site( + site::S, (Nrow, Ncol)::NTuple{2,Int} +) where {S<:Union{CartesianIndex{2},NTuple{2,Int}}} + r, c = site[1], site[2] + return CartesianIndex(1 + Nrow - r, 1 + Ncol - c) +end + +""" + mirror_antidiag(H::LocalOperator) + +Mirror a `LocalOperator` across the anti-diagonal axis of its lattice. +""" +function mirror_antidiag(H::LocalOperator) + lattice2 = mirror_antidiag(H.lattice) + terms2 = ( + (Tuple(_mirror_antidiag_site(site, size(H.lattice)) for site in sites) => op) for + (sites, op) in H.terms + ) + return LocalOperator(lattice2, terms2...) +end + +function Base.rotr90(H::LocalOperator) + lattice2 = rotr90(H.lattice) + terms2 = ( + (Tuple(_rotr90_site(site, size(H.lattice)) for site in sites) => op) for + (sites, op) in H.terms + ) + return LocalOperator(lattice2, terms2...) +end + +function Base.rotl90(H::LocalOperator) + lattice2 = rotl90(H.lattice) + terms2 = ( + (Tuple(_rotl90_site(site, size(H.lattice)) for site in sites) => op) for + (sites, op) in H.terms + ) + return LocalOperator(lattice2, terms2...) +end + +function Base.rot180(H::LocalOperator) + lattice2 = rot180(H.lattice) + terms2 = ( + (Tuple(_rot180_site(site, size(H.lattice)) for site in sites) => op) for + (sites, op) in H.terms + ) + return LocalOperator(lattice2, terms2...) +end diff --git a/src/operators/models.jl b/src/operators/models.jl index 324397ff..1ca82cd3 100644 --- a/src/operators/models.jl +++ b/src/operators/models.jl @@ -1,5 +1,13 @@ ## Model Hamiltonians # ------------------- +""" + nearest_neighbour_hamiltonian( + lattice::Matrix{S}, h::AbstractTensorMap{S,2,2} + ) where {S} + +Create a nearest neighbor `LocalOperator` by specifying the 2-site interaction term `h` +which acts both in horizontal and vertical direction. +""" function nearest_neighbour_hamiltonian( lattice::Matrix{S}, h::AbstractTensorMap{S,2,2} ) where {S} @@ -54,7 +62,7 @@ end """ j1_j2([elt::Type{T}], [symm::Type{S}], [lattice::InfiniteSquare]; - J1=1.0, J2=1.0, spin=1//2, sublattice=true) + J1=1.0, J2=1.0, spin=1//2, sublattice=true) Square lattice J₁-J₂ model. The `sublattice` kwarg enables a single site unit cell via a sublattice rotation. @@ -122,3 +130,52 @@ function pwave_superconductor( (neighbor => hy for neighbor in y_neighbors)..., ) end + +function MPSKitModels.hubbard_model( + T::Type{<:Number}, + particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}, + lattice::InfiniteSquare; + t=1.0, + U=1.0, + mu=0.0, + n::Integer=0, +) + @assert n == 0 "Currently no support for imposing a fixed particle number" + N = MPSKitModels.e_number(T, particle_symmetry, spin_symmetry) + pspace = space(N, 1) + unit = TensorKit.id(pspace) + hopping = + MPSKitModels.e⁺e⁻(T, particle_symmetry, spin_symmetry) + + MPSKitModels.e⁻e⁺(T, particle_symmetry, spin_symmetry) + interaction_term = MPSKitModels.nꜛnꜜ(T, particle_symmetry, spin_symmetry) + site_term = U * interaction_term - mu * N + h = (-t) * hopping + (1 / 4) * (site_term ⊗ unit + unit ⊗ site_term) + return nearest_neighbour_hamiltonian(fill(pspace, size(lattice)), h) +end + +function MPSKitModels.tj_model( + T::Type{<:Number}, + particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}, + lattice::InfiniteSquare; + t=2.5, + J=1.0, + mu=0.0, + slave_fermion::Bool=false, +) + hopping = + TJOperators.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion) + + TJOperators.e_minplus(particle_symmetry, spin_symmetry; slave_fermion) + num = TJOperators.e_number(particle_symmetry, spin_symmetry; slave_fermion) + heis = + TJOperators.S_exchange(particle_symmetry, spin_symmetry; slave_fermion) - + (1 / 4) * (num ⊗ num) + pspace = space(num, 1) + unit = TensorKit.id(pspace) + h = (-t) * hopping + J * heis - (mu / 4) * (num ⊗ unit + unit ⊗ num) + if T <: Real + h = real(h) + end + return nearest_neighbour_hamiltonian(fill(pspace, size(lattice)), h) +end diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 51b353d4..fea4ab1f 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -113,12 +113,13 @@ end Base.size(T::InfinitePEPS) = size(T.A) Base.size(T::InfinitePEPS, i) = size(T.A, i) Base.length(T::InfinitePEPS) = length(T.A) -Base.eltype(T::InfinitePEPS) = eltype(T.A) -VectorInterface.scalartype(T::InfinitePEPS) = scalartype(T.A) +Base.eltype(T::InfinitePEPS) = eltype(typeof(T)) +Base.eltype(::Type{<:InfinitePEPS{T}}) where {T} = T +VectorInterface.scalartype(::Type{T}) where {T<:InfinitePEPS} = scalartype(eltype(T)) ## Copy Base.copy(T::InfinitePEPS) = InfinitePEPS(copy(T.A)) -Base.similar(T::InfinitePEPS, args...) = InfinitePEPS(similar.(T.A, args...)) +Base.similar(T::InfinitePEPS, args...) = InfinitePEPS(similar(T.A, args...)) Base.repeat(T::InfinitePEPS, counts...) = InfinitePEPS(repeat(T.A, counts...)) Base.getindex(T::InfinitePEPS, args...) = Base.getindex(T.A, args...) diff --git a/src/states/infiniteweightpeps.jl b/src/states/infiniteweightpeps.jl new file mode 100644 index 00000000..91303bed --- /dev/null +++ b/src/states/infiniteweightpeps.jl @@ -0,0 +1,223 @@ + +""" + const PEPSWeight{S} + +Default type for PEPS bond weights with 2 virtual indices, conventionally ordered as: ``wt : WS ← EN``. +`WS`, `EN` denote the west/south, east/north spaces for x/y-weights on the square lattice, respectively. +""" +const PEPSWeight{S} = AbstractTensorMap{S,1,1} where {S<:ElementarySpace} + +""" + struct SUWeight{E<:PEPSWeight} + +Schmidt bond weights used in simple/cluster update. +Weight elements are always real. +""" +struct SUWeight{E<:PEPSWeight} + data::Array{E,3} + + function SUWeight(data::Array{E,3}) where {E<:PEPSWeight} + @assert eltype(data[1]) <: Real + return new{E}(data) + end +end + +function SUWeight(wts_mats::AbstractMatrix{E}...) where {E<:PEPSWeight} + n_mat = length(wts_mats) + Nr, Nc = size(wts_mats[1]) + @assert all((Nr, Nc) == size(wts_mat) for wts_mat in wts_mats) + weights = collect( + wts_mats[d][r, c] for (d, r, c) in Iterators.product(1:n_mat, 1:Nr, 1:Nc) + ) + return SUWeight(weights) +end + +## Shape and size +Base.size(W::SUWeight) = size(W.data) +Base.size(W::SUWeight, i) = size(W.data, i) +Base.length(W::SUWeight) = length(W.data) +Base.eltype(W::SUWeight) = eltype(typeof(W)) +Base.eltype(::Type{SUWeight{E}}) where {E} = E +VectorInterface.scalartype(::Type{T}) where {T<:SUWeight} = scalartype(eltype(T)) + +Base.getindex(W::SUWeight, args...) = Base.getindex(W.data, args...) +Base.setindex!(W::SUWeight, args...) = (Base.setindex!(W.data, args...); W) +Base.axes(W::SUWeight, args...) = axes(W.data, args...) + +function compare_weights(wts1::SUWeight, wts2::SUWeight) + @assert size(wts1) == size(wts2) + return sum(_singular_value_distance, zip(wts1.data, wts2.data)) / length(wts1) +end + +""" + struct InfiniteWeightPEPS{T<:PEPSTensor,E<:PEPSWeight} <: AbstractPEPS + +Represents an infinite projected entangled-pair state on a 2D square lattice +consisting of vertex tensors and bond weights. +""" +struct InfiniteWeightPEPS{T<:PEPSTensor,E<:PEPSWeight} <: AbstractPEPS + vertices::Matrix{T} + weights::SUWeight{E} + + function InfiniteWeightPEPS( + vertices::Matrix{T}, weights::SUWeight{E} + ) where {T<:PEPSTensor,E<:PEPSWeight} + @assert size(vertices) == size(weights)[2:end] + Nr, Nc = size(vertices) + for (r, c) in Iterators.product(1:Nr, 1:Nc) + space(weights[2, r, c], 1)' == space(vertices[r, c], 2) || throw( + SpaceMismatch("South space of bond weight y$((r, c)) does not match.") + ) + space(weights[2, r, c], 2)' == space(vertices[_prev(r, Nr), c], 4) || throw( + SpaceMismatch("North space of bond weight y$((r, c)) does not match.") + ) + space(weights[1, r, c], 1)' == space(vertices[r, c], 3) || + throw(SpaceMismatch("West space of bond weight x$((r, c)) does not match.")) + space(weights[1, r, c], 2)' == space(vertices[r, _next(c, Nc)], 5) || + throw(SpaceMismatch("East space of bond weight x$((r, c)) does not match.")) + end + return new{T,E}(vertices, weights) + end +end + +""" + InfiniteWeightPEPS( + vertices::Matrix{T}, weight_mats::Matrix{E}... + ) where {T<:PEPSTensor,E<:PEPSWeight} + +Create an InfiniteWeightPEPS from matrices of vertex tensors, +and separate matrices of weights on each type of bond at all locations in the unit cell. +""" +function InfiniteWeightPEPS( + vertices::Matrix{T}, weight_mats::Matrix{E}... +) where {T<:PEPSTensor,E<:PEPSWeight} + return InfiniteWeightPEPS(vertices, SUWeight(weight_mats...)) +end + +""" + InfiniteWeightPEPS( + f, T, Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) + ) where {S<:ElementarySpace} + +Create an InfiniteWeightPEPS by specifying its physical, north and east spaces (as `ElementarySpace`s) and unit cell size. +Use `T` to specify the element type of the vertex tensors. +Bond weights are initialized as identity matrices of element type `Float64`. +""" +function InfiniteWeightPEPS( + f, T, Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) +) where {S<:ElementarySpace} + vertices = InfinitePEPS(f, T, Pspace, Nspace, Espace; unitcell=unitcell).A + Nr, Nc = unitcell + weights = collect( + id(d == 1 ? Espace : Nspace) for (d, r, c) in Iterators.product(1:2, 1:Nr, 1:Nc) + ) + return InfiniteWeightPEPS(vertices, SUWeight(weights)) +end + +""" + absorb_weight(t::T, row::Int, col::Int, ax::Int, weights::SUWeight; + sqrtwt::Bool=false, invwt::Bool=false) where {T<:PEPSTensor} + +Absorb or remove environment weight on axis `ax` of vertex tensor `t` +known to be located at position (`row`, `col`) in the unit cell. +Weights around the tensor at `(row, col)` are +``` + ↓ + [2,r,c] + ↓ + ← [1,r,c-1] ← T[r,c] ← [1,r,c] ← + ↓ + [1,r+1,c] + ↓ +``` + +# Arguments +- `t::T`: The vertex tensor to which the weight will be absorbed. The first axis of `t` should be the physical axis. +- `row::Int`: The row index specifying the position in the tensor network. +- `col::Int`: The column index specifying the position in the tensor network. +- `ax::Int`: The axis along which the weight is absorbed. +- `weights::SUWeight`: The weight object to absorb into the tensor. +- `sqrtwt::Bool=false` (optional): If `true`, the square root of the weight is absorbed. +- `invwt::Bool=false` (optional): If `true`, the inverse of the weight is absorbed. + +# Details +The optional kwargs `sqrtwt` and `invwt` allow taking the square root or the inverse of the weight before absorption. + +# Examples +```julia +# Absorb the weight into the 2nd axis of tensor at position (2, 3) +absorb_weight(t, 2, 3, 2, weights) + +# Absorb the square root of the weight into the tensor +absorb_weight(t, 2, 3, 2, weights; sqrtwt=true) + +# Absorb the inverse of (i.e. remove) the weight into the tensor +absorb_weight(t, 2, 3, 2, weights; invwt=true) +``` +""" +function absorb_weight( + t::T, + row::Int, + col::Int, + ax::Int, + weights::SUWeight; + sqrtwt::Bool=false, + invwt::Bool=false, +) where {T<:PEPSTensor} + Nr, Nc = size(weights)[2:end] + @assert 1 <= row <= Nr && 1 <= col <= Nc + @assert 2 <= ax <= 5 + pow = (sqrtwt ? 1 / 2 : 1) * (invwt ? -1 : 1) + if ax == 2 # north + wt = weights[2, row, col] + elseif ax == 3 # east + wt = weights[1, row, col] + elseif ax == 4 # south + wt = weights[2, _next(row, Nr), col] + else # west + wt = weights[1, row, _prev(col, Nc)] + end + wt2 = sdiag_pow(wt, pow) + indices_t = collect(-1:-1:-5) + indices_t[ax] = 1 + indices_wt = (ax in (2, 3) ? [1, -ax] : [-ax, 1]) + t2 = permute(ncon((t, wt2), (indices_t, indices_wt)), ((1,), Tuple(2:5))) + return t2 +end + +""" + InfinitePEPS(peps::InfiniteWeightPEPS) + +Create `InfinitePEPS` from `InfiniteWeightPEPS` by absorbing bond weights into vertex tensors. +""" +function InfinitePEPS(peps::InfiniteWeightPEPS) + vertices = deepcopy(peps.vertices) + N1, N2 = size(vertices) + for (r, c) in Iterators.product(1:N1, 1:N2) + for ax in 2:5 + vertices[r, c] = absorb_weight( + vertices[r, c], r, c, ax, peps.weights; sqrtwt=true + ) + end + end + return InfinitePEPS(vertices) +end + +function Base.size(peps::InfiniteWeightPEPS) + return size(peps.vertices) +end + +""" + mirror_antidiag(peps::InfiniteWeightPEPS) + +Mirror the unit cell of an iPEPS with weights by its anti-diagonal line. +""" +function mirror_antidiag(peps::InfiniteWeightPEPS) + vertices2 = mirror_antidiag(peps.vertices) + for (i, t) in enumerate(vertices2) + vertices2[i] = permute(t, ((1,), (3, 2, 5, 4))) + end + weights2_x = mirror_antidiag(peps.weights[2, :, :]) + weights2_y = mirror_antidiag(peps.weights[1, :, :]) + return InfiniteWeightPEPS(vertices2, weights2_x, weights2_y) +end diff --git a/src/utility/mirror.jl b/src/utility/mirror.jl new file mode 100644 index 00000000..40626213 --- /dev/null +++ b/src/utility/mirror.jl @@ -0,0 +1,12 @@ +""" + mirror_antidiag(arr::AbstractMatrix) + +Mirror a matrix by its anti-diagonal line (the 45 degree line through the lower-left corner). + +The element originally at [r, c] is moved [Nc-c+1, Nr-r+1], i.e. the element now at [r, c] +was originally at [Nr-c+1, Nc-r+1] +""" +function mirror_antidiag(arr::AbstractMatrix) + Nr, Nc = size(arr) + return collect(arr[Nr - c + 1, Nc - r + 1] for (r, c) in Iterators.product(1:Nc, 1:Nr)) +end diff --git a/src/utility/svd.jl b/src/utility/svd.jl index a590c87a..cfd9a1f9 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -71,7 +71,7 @@ the iterative SVD didn't converge, the algorithm falls back to a dense SVD. end function random_start_vector(t::Matrix) - return randn(eltype(t), size(t, 1)) + return randn(scalartype(t), size(t, 1)) end # Compute SVD data block-wise using KrylovKit algorithm diff --git a/src/utility/util.jl b/src/utility/util.jl index 0b9fd4c4..45182f1d 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -2,7 +2,6 @@ _next(i, total) = mod1(i + 1, total) _prev(i, total) = mod1(i - 1, total) -# iterator over each coordinates """ eachcoordinate(x, dirs=1:4) @@ -21,38 +20,51 @@ function _elementwise_mult(a::AbstractTensorMap, b::AbstractTensorMap) return dst end -# Compute √S⁻¹ for diagonal TensorMaps -_safe_inv(a, tol) = abs(a) < tol ? zero(a) : inv(a) -function sdiag_inv_sqrt(S::AbstractTensorMap; tol::Real=eps(eltype(S))^(3 / 4)) - tol *= norm(S, Inf) # Relative tol w.r.t. largest singular value (use norm(∘, Inf) to make differentiable) - invsq = similar(S) +_safe_pow(a, pow, tol) = (pow < 0 && abs(a) < tol) ? zero(a) : a^pow + +""" + sdiag_pow(S::AbstractTensorMap, pow::Real; tol::Real=eps(scalartype(S))^(3 / 4)) - if sectortype(S) == Trivial +Compute `S^pow` for diagonal matrices `S`. +""" +function sdiag_pow(S::AbstractTensorMap, pow::Real; tol::Real=eps(scalartype(S))^(3 / 4)) + tol *= norm(S, Inf) # Relative tol w.r.t. largest singular value (use norm(∘, Inf) to make differentiable) + Spow = similar(S) + for (k, b) in blocks(S) copyto!( - invsq.data, - LinearAlgebra.diagm(_safe_inv.(LinearAlgebra.diag(S.data), tol) .^ (1 / 2)), + blocks(Spow)[k], + LinearAlgebra.diagm(_safe_pow.(LinearAlgebra.diag(b), pow, tol)), ) - else - for (k, b) in blocks(S) - copyto!( - blocks(invsq)[k], - LinearAlgebra.diagm(_safe_inv.(LinearAlgebra.diag(b), tol) .^ (1 / 2)), - ) - end end + return Spow +end - return invsq +""" + absorb_s(u::AbstractTensorMap, s::AbstractTensorMap, vh::AbstractTensorMap) + +Given `tsvd` result `u`, `s` and `vh`, absorb singular values `s` into `u` and `vh` by: +``` + u -> u * sqrt(s), vh -> sqrt(s) * vh +``` +""" +function absorb_s(u::AbstractTensorMap, s::AbstractTensorMap, vh::AbstractTensorMap) + sqrt_s = sdiag_pow(s, 0.5) + return u * sqrt_s, sqrt_s * vh end function ChainRulesCore.rrule( - ::typeof(sdiag_inv_sqrt), S::AbstractTensorMap; tol::Real=eps(eltype(S))^(3 / 4) + ::typeof(sdiag_pow), + S::AbstractTensorMap, + pow::Real; + tol::Real=eps(scalartype(S))^(3 / 4), ) tol *= norm(S, Inf) - invsq = sdiag_inv_sqrt(S; tol) - function sdiag_inv_sqrt_pullback(c̄) - return (ChainRulesCore.NoTangent(), -1 / 2 * _elementwise_mult(c̄, invsq'^3)) + spow = sdiag_pow(S, pow; tol) + spow_minus1_conj = scale!(sdiag_pow(S', pow - 1; tol), pow) + function sdiag_pow_pullback(c̄) + return (ChainRulesCore.NoTangent(), _elementwise_mult(c̄, spow_minus1_conj)) end - return invsq, sdiag_inv_sqrt_pullback + return spow, sdiag_pow_pullback end # Check whether diagonals contain degenerate values up to absolute or relative tolerance diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 1c05b084..a8bbc426 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -1,40 +1,104 @@ using Test using Random +using Accessors using PEPSKit using TensorKit using KrylovKit using OptimKit # initialize parameters -χbond = 2 +Dbond = 2 χenv = 16 ctm_alg = CTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) +# compare against Juraj Hasik's data: +# https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.0/state_1s_A1_j20.0_D2_chi_opt48.dat +E_ref = -0.6602310934799577 -# initialize states -Random.seed!(91283219347) -H = heisenberg_XYZ(InfiniteSquare()) -psi_init = InfinitePEPS(2, χbond) -env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg) - -# find fixedpoint -result = fixedpoint(psi_init, H, opt_alg, env_init) -ξ_h, ξ_v, = correlation_length(result.peps, result.env) - -@test result.E ≈ -0.6694421 atol = 1e-2 -@test all(@. ξ_h > 0 && ξ_v > 0) - -# same test but for 1x2 unit cell -unitcell = (1, 2) -H_1x2 = heisenberg_XYZ(InfiniteSquare(unitcell...)) -psi_init_1x2 = InfinitePEPS(2, χbond; unitcell) -env_init_1x2 = leading_boundary( - CTMRGEnv(psi_init_1x2, ComplexSpace(χenv)), psi_init_1x2, ctm_alg -) -result_1x2 = fixedpoint(psi_init_1x2, H_1x2, opt_alg, env_init_1x2) -ξ_h_1x2, ξ_v_1x2, = correlation_length(result_1x2.peps, result_1x2.env) +@testset "(1, 1) unit cell AD optimization" begin + # initialize states + Random.seed!(123) + H = heisenberg_XYZ(InfiniteSquare()) + psi_init = InfinitePEPS(2, Dbond) + env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg) + + # optimize energy and compute correlation lengths + result = fixedpoint(psi_init, H, opt_alg, env_init) + ξ_h, ξ_v, = correlation_length(result.peps, result.env) + + @test result.E ≈ E_ref atol = 1e-2 + @test all(@. ξ_h > 0 && ξ_v > 0) +end + +@testset "(1, 2) unit cell AD optimization" begin + # initialize states + Random.seed!(456) + unitcell = (1, 2) + H_1x2 = heisenberg_XYZ(InfiniteSquare(unitcell...)) + psi_init_1x2 = InfinitePEPS(2, Dbond; unitcell) + env_init_1x2 = leading_boundary( + CTMRGEnv(psi_init_1x2, ComplexSpace(χenv)), psi_init_1x2, ctm_alg + ) + + # optimize energy and compute correlation lengths + result_1x2 = fixedpoint(psi_init_1x2, H_1x2, opt_alg, env_init_1x2) + ξ_h_1x2, ξ_v_1x2, = correlation_length(result_1x2.peps, result_1x2.env) + + @test result_1x2.E ≈ 2 * E_ref atol = 1e-2 + @test all(@. ξ_h_1x2 > 0 && ξ_v_1x2 > 0) +end + +@testset "Simple update into AD optimization" begin + # random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers) + Random.seed!(789) + N1, N2 = 2, 2 + Pspace = ℂ^2 + Vspace = ℂ^Dbond + Espace = ℂ^χenv + 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 + # Heisenberg model Hamiltonian (already only includes nearest neighbor terms) + ham = heisenberg_XYZ(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)...) + + # simple update + dts = [1e-2, 1e-3, 4e-4, 1e-4] + tols = [1e-7, 1e-8, 1e-8, 1e-8] + maxiter = 5000 + for (n, (dt, tol)) in enumerate(zip(dts, tols)) + Dbond2 = (n == 1) ? Dbond + 2 : Dbond + trscheme = truncerr(1e-10) & truncdim(Dbond2) + alg = SimpleUpdate(dt, tol, maxiter, trscheme) + result = simpleupdate(peps, ham, alg; bipartite=false) + peps = result[1] + end + + # absorb weight into site tensors and CTMRG + peps = InfinitePEPS(peps) + envs = CTMRGEnv(rand, Float64, peps, Espace) + trscheme = truncerr(1e-9) & truncdim(χenv) + envs = leading_boundary(envs, peps, CTMRG(; trscheme, ctmrgscheme=:sequential)) + + # measure physical quantities + e_site = costfun(peps, envs, ham) / (N1 * N2) + @info "Simple update energy = $e_site" + # benchmark data from Phys. Rev. B 94, 035133 (2016) + @test isapprox(e_site, -0.6594; atol=1e-3) -@test result_1x2.E ≈ 2 * result.E atol = 1e-2 -@test all(@. ξ_h_1x2 > 0 && ξ_v_1x2 > 0) + # continue with auto differentiation + svd_alg_gmres = SVDAdjoint(; rrule_alg=GMRES(; tol=1e-8)) + opt_alg_gmres = @set opt_alg.boundary_alg.projector_alg.svd_alg = svd_alg_gmres + result = fixedpoint(peps, ham, opt_alg_gmres, envs) # sensitivity warnings and degeneracies due to SU(2)? + ξ_h, ξ_v, = correlation_length(result.peps, result.env) + e_site2 = result.E / (N1 * N2) + @info "Auto diff energy = $e_site2" + @test e_site2 ≈ E_ref atol = 1e-2 + @test all(@. ξ_h > 0 && ξ_v > 0) +end diff --git a/test/runtests.jl b/test/runtests.jl index 07a4cd2f..b15ad1c6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ end @time @safetestset "Heisenberg model" begin include("heisenberg.jl") end - @time @safetestset "Heisenberg model" begin + @time @safetestset "J1-J2 model" begin include("j1j2_model.jl") end @time @safetestset "P-wave superconductor" begin