diff --git a/README.md b/README.md index ebb55861..25f9a924 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/docs/src/index.md b/docs/src/index.md index 465b150f..709f4222 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -26,7 +26,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # sublattice rotation t # configuring the parameters D = 2 chi = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=1, trscheme=truncdim(chi)) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi)) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl index eb699b21..51af6dea 100644 --- a/examples/boundary_mps.jl +++ b/examples/boundary_mps.jl @@ -32,7 +32,9 @@ mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) N = abs(prod(expectation_value(mps, T))) # This can be compared to the result obtained using the CTMRG algorithm -ctm = leading_boundary(peps, CTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))) +ctm = leading_boundary( + peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20)) +) N´ = abs(norm(peps, ctm)) @show abs(N - N´) / N @@ -53,7 +55,9 @@ mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2)) mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS()) N2 = abs(prod(expectation_value(mps2, T2))) -ctm2 = leading_boundary(peps2, CTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))) +ctm2 = leading_boundary( + peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20)) +) N2´ = abs(norm(peps2, ctm2)) @show abs(N2 - N2´) / N2 diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index 25b30c02..34baed18 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -11,7 +11,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1) # Parameters χbond = 2 χenv = 20 -ctm_alg = CTMRG(; tol=1e-10, miniter=4, maxiter=100, verbosity=2) +ctm_alg = SimultaneousCTMRG(; tol=1e-10, verbosity=2) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2), diff --git a/examples/hubbard_su.jl b/examples/hubbard_su.jl index 4bd84846..6478a28e 100644 --- a/examples/hubbard_su.jl +++ b/examples/hubbard_su.jl @@ -1,5 +1,4 @@ using Test -using Printf using Random using PEPSKit using TensorKit @@ -14,15 +13,16 @@ if symm == Trivial 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) +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] @@ -31,32 +31,29 @@ 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] + peps, = simpleupdate(peps, ham, alg; bipartite=false) end -# absort weight into site tensors +# 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] - 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) + 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. -""" -# measure physical quantities -E = costfun(peps, envs, ham) / (N1 * N2) +# 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 -@info @sprintf("Energy = %.8f\n", E) -@info @sprintf("Benchmark energy = %.8f\n", E_exact) -@test isapprox(E, E_exact; atol=1e-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) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index a69e5331..74d68e03 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -3,7 +3,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Base: @kwdef using Compat -using Accessors +using Accessors: @set using VectorInterface using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations using ChainRulesCore, Zygote @@ -43,6 +43,9 @@ include("algorithms/contractions/ctmrg_contractions.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") +include("algorithms/ctmrg/projectors.jl") +include("algorithms/ctmrg/simultaneous.jl") +include("algorithms/ctmrg/sequential.jl") include("algorithms/ctmrg/gaugefix.jl") include("algorithms/time_evolution/gatetools.jl") @@ -61,12 +64,16 @@ include("utility/symmetrization.jl") const ctmrg_tol = 1e-8 const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 - const ctmrgscheme = :simultaneous const reuse_env = true const trscheme = FixedSpaceTruncation() 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 projector_alg_type = HalfInfiniteProjector + const projector_alg = projector_alg_type(svd_alg, trscheme, 2) + const ctmrg_alg = SimultaneousCTMRG( + ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg + ) const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -83,12 +90,14 @@ Module containing default values that represent typical algorithm parameters. - `ctmrg_tol`: Tolerance checking singular value and norm convergence - `fpgrad_maxiter`: Maximal number of iterations for computing the CTMRG fixed-point gradient - `fpgrad_tol`: Convergence tolerance for the fixed-point gradient iteration -- `ctmrgscheme`: Scheme for growing, projecting and renormalizing CTMRG environments - `reuse_env`: If `true`, the current optimization step is initialized on the previous environment - `trscheme`: Truncation scheme for SVDs and other decompositions - `fwd_alg`: SVD algorithm that is used in the forward pass - `rrule_alg`: Reverse-rule for differentiating that SVD - `svd_alg`: Combination of `fwd_alg` and `rrule_alg` +- `projector_alg_type`: Default type of projector algorithm +- `projector_alg`: Algorithm to compute CTMRG projectors +- `ctmrg_alg`: Algorithm for performing CTMRG runs - `optimizer`: Optimization algorithm for PEPS ground-state optimization - `gradient_linsolver`: Default linear solver for the `LinSolver` gradient algorithm - `iterscheme`: Scheme for differentiating one CTMRG iteration @@ -97,19 +106,28 @@ Module containing default values that represent typical algorithm parameters. """ module Defaults using TensorKit, KrylovKit, OptimKit, OhMyThreads - using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint + using PEPSKit: + LinSolver, + FixedSpaceTruncation, + SVDAdjoint, + HalfInfiniteProjector, + SimultaneousCTMRG + const ctmrg_tol = 1e-8 const ctmrg_maxiter = 100 const ctmrg_miniter = 4 - const ctmrg_tol = 1e-8 const fpgrad_maxiter = 30 const fpgrad_tol = 1e-6 - const ctmrgscheme = :simultaneous const sparse = false const reuse_env = true const trscheme = FixedSpaceTruncation() 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 projector_alg_type = HalfInfiniteProjector + const projector_alg = projector_alg_type(svd_alg, trscheme, 2) + const ctmrg_alg = SimultaneousCTMRG( + ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg + ) const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2) const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol @@ -165,9 +183,10 @@ end using .Defaults: set_scheduler! export set_scheduler! export SVDAdjoint, IterSVD, NonTruncSVDAdjoint -export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length +export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG +export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector export LocalOperator -export expectation_value, costfun, product_peps +export expectation_value, costfun, product_peps, correlation_length export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver export fixedpoint diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl index cb920d76..2fc43bdd 100644 --- a/src/algorithms/contractions/ctmrg_contractions.jl +++ b/src/algorithms/contractions/ctmrg_contractions.jl @@ -210,13 +210,40 @@ function right_projector(E_1, C, E_2, U, isqS, ket::PEPSTensor, bra::PEPSTensor= end """ - halfinfinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) - halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} - halfinfinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} - halfinfinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, - ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + contract_projectors(U, S, V, Q, Q_next) + +Compute projectors based on a SVD of `Q * Q_next`, where the inverse square root +`isqS` of the singular values is computed. + +Left projector: +``` + -- |~~~~~~| -- |~~| + |Q_next| |V'| -- isqS -- + == |~~~~~~| == |~~| +``` + +Right projector: +``` + |~~| -- |~~~| -- + -- isqS -- |U'| | Q | + |~~| == |~~~| == +``` +""" +function contract_projectors(U, S, V, Q, Q_next) + 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 +end + +""" + half_infinite_environment(quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3}) + half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + half_infinite_environment(C_1, C_2, E_1, E_2, E_3, E_4, x, + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} + half_infinite_environment(x, C_1, C_2, E_1, E_2, E_3, E_4, + ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2) where {P<:PEPSTensor} Contract two quadrants (enlarged corners) to form a half-infinite environment. @@ -236,7 +263,7 @@ The environment can also be contracted directly from all its constituent tensors | || || | ``` -Alternatively, contract environment with a vector `x` acting on it +Alternatively, contract the environment with a vector `x` acting on it ``` C_1 -- E_2 -- E_3 -- C_2 @@ -244,19 +271,18 @@ Alternatively, contract environment with a vector `x` acting on it E_1 == ket_bra_1 == ket_bra_2 == E_4 | || || | [~~~~~~x~~~~~~] - || | ``` or contract the adjoint environment with `x`, e.g. as needed for iterative solvers. """ -function halfinfinite_environment( +function half_infinite_environment( quadrant1::AbstractTensorMap{S,3,3}, quadrant2::AbstractTensorMap{S,3,3} ) where {S} return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := quadrant1[χ_in D_inabove D_inbelow; χ D1 D2] * quadrant2[χ D1 D2; χ_out D_outabove D_outbelow] end -function halfinfinite_environment( +function half_infinite_environment( C_1, C_2, E_1, E_2, E_3, E_4, ket_1::P, ket_2::P, bra_1::P=ket_1, bra_2::P=ket_2 ) where {P<:PEPSTensor} return @autoopt @tensor env[χ_in D_inabove D_inbelow; χ_out D_outabove D_outbelow] := @@ -271,7 +297,7 @@ function halfinfinite_environment( C_2[χ4; χ5] * E_4[χ5 D7 D8; χ_out] end -function halfinfinite_environment( +function half_infinite_environment( C_1, C_2, E_1, @@ -297,7 +323,7 @@ function halfinfinite_environment( E_4[χ5 D7 D8; χ6] * x[χ6 D11 D12] end -function halfinfinite_environment( +function half_infinite_environment( x::AbstractTensor{S,3}, C_1, C_2, @@ -310,7 +336,7 @@ function halfinfinite_environment( bra_1::P=ket_1, bra_2::P=ket_2, ) where {S,P<:PEPSTensor} - return @autoopt @tensor env_x[χ_in D_inabove D_inbelow] := + return @autoopt @tensor x_env[χ_in D_inabove D_inbelow] := x[χ1 D1 D2] * conj(E_1[χ1 D3 D4; χ2]) * conj(C_1[χ2; χ3]) * diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index 317a7231..ac5444f2 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -1,109 +1,37 @@ -""" - FixedSpaceTruncation <: TensorKit.TruncationScheme -CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD -is performed fixed. Since different environment directions and unit cell entries might -have different spaces, this truncation style is different from `TruncationSpace`. """ -struct FixedSpaceTruncation <: TensorKit.TruncationScheme end + CTMRGAlgorithm +Abstract super type for the corner transfer matrix renormalization group (CTMRG) algorithm +for contracting infinite PEPS. """ - struct ProjectorAlg{S}(; svd_alg=TensorKit.SVD(), trscheme=TensorKit.notrunc(), - fixedspace=false, verbosity=0) +abstract type CTMRGAlgorithm end -Algorithm struct collecting all projector related parameters. The truncation scheme has to be -a `TensorKit.TruncationScheme`, and some SVD algorithms might have further restrictions on what -kind of truncation scheme can be used. If `fixedspace` is true, the truncation scheme is set to -`truncspace(V)` where `V` is the environment bond space, adjusted to the corresponding -environment direction/unit cell entry. """ -@kwdef struct ProjectorAlg{S<:SVDAdjoint,T} - svd_alg::S = Defaults.svd_alg - trscheme::T = Defaults.trscheme - verbosity::Int = 0 -end -# TODO: add option for different projector styles (half-infinite, full-infinite, etc.) - -function truncation_scheme(alg::ProjectorAlg, edge) - if alg.trscheme isa FixedSpaceTruncation - return truncspace(space(edge, 1)) - else - return alg.trscheme - end -end - -function svd_algorithm(alg::ProjectorAlg, (dir, r, c)) - if alg.svd_alg isa SVDAdjoint{<:FixedSVD} - fwd_alg = alg.svd_alg.fwd_alg - fix_svd = FixedSVD(fwd_alg.U[dir, r, c], fwd_alg.S[dir, r, c], fwd_alg.V[dir, r, c]) - return SVDAdjoint(; fwd_alg=fix_svd, rrule_alg=alg.svd_alg.rrule_alg) - else - return alg.svd_alg - end -end + ctmrg_iteration(state, env, alg::CTMRGAlgorithm) -> env′, info +Perform a single CTMRG iteration in which all directions are being grown and renormalized. """ - CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, verbosity=0, - svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation(), - ctmrgscheme=Defaults.ctmrgscheme) - -Algorithm struct that represents the CTMRG algorithm for contracting infinite PEPS. -Each CTMRG run is converged up to `tol` where the singular value convergence of the -corners as well as the norm is checked. The maximal and minimal number of CTMRG iterations -is set with `maxiter` and `miniter`. Different levels of output information are printed -depending on `verbosity`, where `0` suppresses all output, `1` only prints warnings, `2` -gives information at the start and end, and `3` prints information every iteration. - -The projectors are computed from `svd_alg` SVDs where the truncation scheme is set via -`trscheme`. +function ctmrg_iteration(state, env, alg::CTMRGAlgorithm) end -In general, two different schemes can be selected with `ctmrgscheme` which determine how -CTMRG is implemented. It can either be `:sequential`, where the projectors are succesively -computed on the western side, and then applied and rotated. Or with `:simultaneous` all projectors -are computed and applied simultaneously on all sides, where in particular the corners get -contracted with two projectors at the same time. """ -struct CTMRG{S} - tol::Float64 - maxiter::Int - miniter::Int - verbosity::Int - projector_alg::ProjectorAlg -end -function CTMRG(; - tol=Defaults.ctmrg_tol, - maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, - verbosity=2, - svd_alg=Defaults.svd_alg, - trscheme=Defaults.trscheme, - ctmrgscheme::Symbol=Defaults.ctmrgscheme, -) - return CTMRG{ctmrgscheme}( - tol, maxiter, miniter, verbosity, ProjectorAlg(; svd_alg, trscheme, verbosity) - ) -end + MPSKit.leading_boundary([envinit], state, alg::CTMRGAlgorithm) -ctmrgscheme(::CTMRG{S}) where {S} = S +Contract `state` using CTMRG and return the CTM environment. Per default, a random +initial environment is used. -# aliases for the different CTMRG schemes -const SequentialCTMRG = CTMRG{:sequential} -const SimultaneousCTMRG = CTMRG{:simultaneous} +Each CTMRG run is converged up to `alg.tol` where the singular value convergence +of the corners and edges is checked. The maximal and minimal number of CTMRG +iterations is set with `alg.maxiter` and `alg.miniter`. -# supply correct constructor for Accessors.@set -Accessors.constructorof(::Type{CTMRG{S}}) where {S} = CTMRG{S} - -""" - MPSKit.leading_boundary([envinit], state, alg::CTMRG) - -Contract `state` using CTMRG and return the CTM environment. -Per default, a random initial environment is used. +Different levels of output information are printed depending on `alg.verbosity`, where `0` +suppresses all output, `1` only prints warnings, `2` gives information at the start and +end, and `3` prints information every iteration. """ -function MPSKit.leading_boundary(state, alg::CTMRG) +function MPSKit.leading_boundary(state, alg::CTMRGAlgorithm) return MPSKit.leading_boundary(CTMRGEnv(state, oneunit(spacetype(state))), state, alg) end -function MPSKit.leading_boundary(envinit, state, alg::CTMRG) +function MPSKit.leading_boundary(envinit, state, alg::CTMRGAlgorithm) CS = map(x -> tsvd(x)[2], envinit.corners) TS = map(x -> tsvd(x)[2], envinit.edges) @@ -115,7 +43,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) return LoggingExtras.withlevel(; alg.verbosity) do ctmrg_loginit!(log, η, N) for iter in 1:(alg.maxiter) - env, = ctmrg_iter(state, env, alg) # Grow and renormalize in all 4 directions + env, = ctmrg_iteration(state, env, alg) # Grow and renormalize in all 4 directions η, CS, TS = calc_convergence(env, CS, TS) N = norm(state, env) @@ -133,51 +61,7 @@ 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 - -Perform one iteration of CTMRG that maps the `state` and `envs` to a new environment, -and also returns the `info` `NamedTuple`. -""" -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 - 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) - enlarged_envs = ctmrg_expand(eachcoordinate(state, 1:4), state, envs) - projectors, info = ctmrg_projectors(enlarged_envs, envs, alg) - envs′ = ctmrg_renormalize(enlarged_envs, projectors, state, envs, alg) - return envs′, info -end - +# custom CTMRG logging ctmrg_loginit!(log, η, N) = @infov 2 loginit!(log, η, N) ctmrg_logiter!(log, iter, η, N) = @infov 3 logiter!(log, iter, η, N) ctmrg_logfinish!(log, iter, η, N) = @infov 2 logfinish!(log, iter, η, N) @@ -188,214 +72,42 @@ ctmrg_logcancel!(log, iter, η, N) = @warnv 1 logcancel!(log, iter, η, N) @non_differentiable ctmrg_logfinish!(args...) @non_differentiable ctmrg_logcancel!(args...) -# ======================================================================================== # -# Expansion step -# ======================================================================================== # - -""" - ctmrg_expand(coordinates, state, envs) - -Expand the environment by absorbing a new PEPS tensor on the given coordinates. -""" -function ctmrg_expand(coordinates, state, envs::CTMRGEnv) - return dtmap(idx -> TensorMap(EnlargedCorner(state, envs, idx), idx[1]), coordinates) -end - -# ======================================================================================== # -# Projector step -# ======================================================================================== # - -""" - ctmrg_projectors(col::Int, enlarged_envs, env, alg::CTMRG{:sequential}) - ctmrg_projectors(enlarged_envs, env, alg::CTMRG{:simultaneous}) - -Compute the CTMRG projectors based on enlarged environments. -In the `:sequential` mode the projectors are computed for the column `col`, whereas -in the `:simultaneous` mode, all projectors (and corresponding SVDs) are computed in parallel. -""" -function ctmrg_projectors( - col::Int, state::InfinitePEPS, envs::CTMRGEnv{C,E}, alg::SequentialCTMRG -) where {C,E} - projector_alg = alg.projector_alg - ϵ = zero(real(scalartype(envs))) - - # SVD half-infinite environment - coordinates = eachcoordinate(envs)[:, col] - projectors = dtmap(coordinates) do (r, c) - r′ = _prev(r, size(envs.corners, 2)) - Q1 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r, c)), SOUTHWEST) - Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST) - QQ = halfinfinite_environment(Q1, Q2) - trscheme = truncation_scheme(projector_alg, envs.edges[WEST, r′, c]) - svd_alg = svd_algorithm(projector_alg, (WEST, r, c)) - U, S, V, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - ϵ = max(ϵ, ϵ_local / norm(S)) - - # Compute SVD truncation error and check for degenerate singular values - Zygote.isderiving() && ignore_derivatives() do - if alg.verbosity > 0 && is_degenerate_spectrum(S) - svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S)) - @warn("degenerate singular values detected: ", svals) - end - end - - # Compute projectors - return build_projectors(U, S, V, Q1, Q2) - end - return (map(first, projectors), map(last, projectors)), (; err=ϵ) -end -function ctmrg_projectors( - enlarged_envs, envs::CTMRGEnv{C,E}, alg::SimultaneousCTMRG -) where {C,E} - projector_alg = alg.projector_alg - # pre-allocation - U, V = Zygote.Buffer.(projector_type(envs.edges)) - # Corner type but with real numbers - S = Zygote.Buffer(U.data, tensormaptype(spacetype(C), 1, 1, real(scalartype(E)))) - - ϵ = zero(real(scalartype(envs))) - coordinates = eachcoordinate(envs, 1:4) - projectors = dtmap(coordinates) do (dir, r, c) - # Row-column index of next enlarged corner - next_rc = if dir == 1 - (r, _next(c, size(envs.corners, 3))) - elseif dir == 2 - (_next(r, size(envs.corners, 2)), c) - elseif dir == 3 - (r, _prev(c, size(envs.corners, 3))) - elseif dir == 4 - (_prev(r, size(envs.corners, 2)), c) - end - - # SVD half-infinite environment - QQ = halfinfinite_environment( - enlarged_envs[dir, r, c], enlarged_envs[_next(dir, 4), next_rc...] - ) - - trscheme = truncation_scheme(projector_alg, envs.edges[dir, next_rc...]) - svd_alg = svd_algorithm(projector_alg, (dir, r, c)) - U_local, S_local, V_local, ϵ_local = PEPSKit.tsvd!(QQ, svd_alg; trunc=trscheme) - U[dir, r, c] = U_local - S[dir, r, c] = S_local - V[dir, r, c] = V_local - ϵ = max(ϵ, ϵ_local / norm(S_local)) - - # Compute SVD truncation error and check for degenerate singular values - Zygote.isderiving() && ignore_derivatives() do - if alg.verbosity > 0 && is_degenerate_spectrum(S_local) - svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S_local)) - @warn("degenerate singular values detected: ", svals) - end - end - - # Compute projectors - return build_projectors( - U_local, - S_local, - V_local, - enlarged_envs[dir, r, c], - enlarged_envs[_next(dir, 4), next_rc...], - ) +#= +In order to compute an error measure, we compare the singular values of the current iteration with the previous one. +However, when the virtual spaces change, this comparison is not directly possible. +Instead, we project both tensors into the smaller space and then compare the difference. + +TODO: we might want to consider embedding the smaller tensor into the larger space and then compute the difference +=# +function _singular_value_distance((S₁, S₂)) + V₁ = space(S₁, 1) + V₂ = space(S₂, 1) + if V₁ == V₂ + return norm(S₁ - S₂) + else + V = infimum(V₁, V₂) + e1 = isometry(V₁, V) + e2 = isometry(V₂, V) + return norm(e1' * S₁ * e1 - e2' * S₂ * e2) end - - P_left = map(first, projectors) - P_right = map(last, projectors) - return (P_left, P_right), (; err=ϵ, U=copy(U), S=copy(S), V=copy(V)) -end - -""" - build_projectors(U::AbstractTensorMap{E,3,1}, S::AbstractTensorMap{E,1,1}, V::AbstractTensorMap{E,1,3}, - Q::AbstractTensorMap{E,3,3}, Q_next::AbstractTensorMap{E,3,3}) where {E<:ElementarySpace} - build_projectors(U::AbstractTensorMap{E,3,1}, S::AbstractTensorMap{E,1,1}, V::AbstractTensorMap{E,1,3}, - Q::EnlargedCorner, Q_next::EnlargedCorner) where {E<:ElementarySpace} - -Construct left and right projectors where the higher-dimensional is facing left and right, respectively. -""" -function build_projectors( - U::AbstractTensorMap{E,3,1}, - S::AbstractTensorMap{E,1,1}, - V::AbstractTensorMap{E,1,3}, - Q::AbstractTensorMap{E,3,3}, - Q_next::AbstractTensorMap{E,3,3}, -) where {E<:ElementarySpace} - isqS = sdiag_pow(S, -0.5) - P_left = Q_next * V' * isqS - P_right = isqS * U' * Q - return P_left, P_right -end -function build_projectors( - U::AbstractTensorMap{E,3,1}, - S::AbstractTensorMap{E,1,1}, - V::AbstractTensorMap{E,1,3}, - Q::EnlargedCorner, - Q_next::EnlargedCorner, -) where {E<:ElementarySpace} - 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 - ) - return P_left, P_right end -# ======================================================================================== # -# Renormalization step -# ======================================================================================== # - """ - ctmrg_renormalize(col::Int, projectors, state, envs, ::CTMRG{:sequential}) - ctmrg_renormalize(enlarged_envs, projectors, state, envs, ::CTMRG{:simultaneous}) + calc_convergence(envs, CSold, TSold) -Apply projectors to renormalize corners and edges. -The `:sequential` mode renormalizes the environment on the column `col`, where as the -`:simultaneous` mode renormalizes all environment tensors simultaneously. +Given a new environment `envs` and the singular values of previous corners and edges +`CSold` and `TSold`, compute the maximal singular value distance. """ -function ctmrg_renormalize(col::Int, projectors, state, envs, ::SequentialCTMRG) - corners = Zygote.Buffer(envs.corners) - edges = Zygote.Buffer(envs.edges) +function calc_convergence(envs, CSold, TSold) + CSnew = map(x -> tsvd(x)[2], envs.corners) + ΔCS = maximum(_singular_value_distance, zip(CSold, CSnew)) - for (dir, r, c) in eachcoordinate(state, 1:4) - (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue - corners[dir, r, c] = envs.corners[dir, r, c] - end - for (dir, r, c) in eachcoordinate(state, 1:4) - (c == col && dir == WEST) && continue - edges[dir, r, c] = envs.edges[dir, r, c] - end + TSnew = map(x -> tsvd(x)[2], envs.edges) + ΔTS = maximum(_singular_value_distance, zip(TSold, TSnew)) - # Apply projectors to renormalize corners and edge - for row in axes(envs.corners, 2) - C_southwest = renormalize_bottom_corner((row, col), envs, projectors) - corners[SOUTHWEST, row, col] = C_southwest / norm(C_southwest) + @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS" - C_northwest = renormalize_top_corner((row, col), envs, projectors) - corners[NORTHWEST, row, col] = C_northwest / norm(C_northwest) - - E_west = renormalize_west_edge((row, col), envs, projectors, state) - edges[WEST, row, col] = E_west / norm(E_west) - end - - return CTMRGEnv(copy(corners), copy(edges)) + return max(ΔCS, ΔTS), CSnew, TSnew end -function ctmrg_renormalize(enlarged_envs, projectors, state, envs, ::SimultaneousCTMRG) - P_left, P_right = projectors - coordinates = eachcoordinate(envs, 1:4) - corners_edges = dtmap(coordinates) do (dir, r, c) - if dir == NORTH - corner = renormalize_northwest_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_north_edge((r, c), envs, P_left, P_right, state) - elseif dir == EAST - corner = renormalize_northeast_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_east_edge((r, c), envs, P_left, P_right, state) - elseif dir == SOUTH - corner = renormalize_southeast_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_south_edge((r, c), envs, P_left, P_right, state) - elseif dir == WEST - corner = renormalize_southwest_corner((r, c), enlarged_envs, P_left, P_right) - edge = renormalize_west_edge((r, c), envs, P_left, P_right, state) - end - return corner / norm(corner), edge / norm(edge) - end - return CTMRGEnv(map(first, corners_edges), map(last, corners_edges)) -end +@non_differentiable calc_convergence(args...) diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index 338395b2..bedc784d 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -141,51 +141,6 @@ function fix_global_phases(envprev::CTMRGEnv, envfix::CTMRGEnv) return CTMRGEnv(cornersgfix, edgesgfix) end -#= -In order to compute an error measure, we compare the singular values of the current iteration with the previous one. -However, when the virtual spaces change, this comparison is not directly possible. -Instead, we project both tensors into the smaller space and then compare the difference. - -TODO: we might want to consider embedding the smaller tensor into the larger space and then compute the difference -=# -function _singular_value_distance((S₁, S₂)) - V₁ = space(S₁, 1) - V₂ = space(S₂, 1) - if V₁ == V₂ - return norm(S₁ - S₂) - else - V = infimum(V₁, V₂) - e1 = isometry(V₁, V) - e2 = isometry(V₂, V) - return norm(e1' * S₁ * e1 - e2' * S₂ * e2) - end -end - -function calc_convergence(envs, CSold, TSold) - CSnew = map(x -> tsvd(x)[2], envs.corners) - ΔCS = maximum(_singular_value_distance, zip(CSold, CSnew)) - - TSnew = map(x -> tsvd(x)[2], envs.edges) - ΔTS = maximum(_singular_value_distance, zip(TSold, TSnew)) - - @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS" - - 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...) - """ calc_elementwise_convergence(envfinal, envfix; atol=1e-6) diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl new file mode 100644 index 00000000..7d7ca4aa --- /dev/null +++ b/src/algorithms/ctmrg/projectors.jl @@ -0,0 +1,104 @@ +""" + FixedSpaceTruncation <: TensorKit.TruncationScheme + +CTMRG specific truncation scheme for `tsvd` which keeps the bond space on which the SVD +is performed fixed. Since different environment directions and unit cell entries might +have different spaces, this truncation style is different from `TruncationSpace`. +""" +struct FixedSpaceTruncation <: TensorKit.TruncationScheme end + +""" + ProjectorAlgorithm + +Abstract super type for all CTMRG projector algorithms. +""" +abstract type ProjectorAlgorithm end + +function svd_algorithm(alg::ProjectorAlgorithm, (dir, r, c)) + if alg.svd_alg isa SVDAdjoint{<:FixedSVD} + fwd_alg = alg.svd_alg.fwd_alg + fix_svd = FixedSVD(fwd_alg.U[dir, r, c], fwd_alg.S[dir, r, c], fwd_alg.V[dir, r, c]) + return SVDAdjoint(; fwd_alg=fix_svd, rrule_alg=alg.svd_alg.rrule_alg) + else + return alg.svd_alg + end +end + +function truncation_scheme(alg::ProjectorAlgorithm, edge) + if alg.trscheme isa FixedSpaceTruncation + return truncspace(space(edge, 1)) + else + return alg.trscheme + end +end + +""" + struct HalfInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the half-infinite CTMRG environment. +""" +@kwdef struct HalfInfiniteProjector{S<:SVDAdjoint,T} <: ProjectorAlgorithm + svd_alg::S = Defaults.svd_alg + trscheme::T = Defaults.trscheme + verbosity::Int = 0 +end + +""" + struct FullInfiniteProjector{S,T}(; svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, verbosity=0) + +Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG environment. +""" +@kwdef struct FullInfiniteProjector{S<:SVDAdjoint,T} <: ProjectorAlgorithm + svd_alg::S = Defaults.svd_alg + trscheme::T = Defaults.trscheme + verbosity::Int = 0 +end + +""" + compute_projector(enlarged_corners, coordinate, alg::ProjectorAlgorithm) + +Determine left and right projectors at the bond given determined by the enlarged corners +and the given coordinate using the specified `alg`. +""" +function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjector) + # SVD half-infinite environment + halfinf = half_infinite_environment(enlarged_corners...) + svd_alg = svd_algorithm(alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=alg.trscheme) + + # Compute SVD truncation error and check for degenerate singular values + Zygote.isderiving() && ignore_derivatives() do + if alg.verbosity > 0 && is_degenerate_spectrum(S) + svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S)) + @warn("degenerate singular values detected: ", svals) + end + end + + P_left, P_right = contract_projectors(U, S, V, enlarged_corners...) + return (P_left, P_right), (; err, U, S, V) +end +function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector) + # QR left and right half-infinite environments + halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2]) + halfinf_right = half_infinite_environment(enlarged_corners[3], enlarged_corners[4]) + _, R_left = leftorth!(halfinf_left) + L_right, _ = rightorth!(halfinf_right) + + # SVD product of QRs + fullinf = R_left * L_right + svd_alg = svd_algorithm(alg, coordinate) + U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=alg.trscheme) + + # Compute SVD truncation error and check for degenerate singular values + Zygote.isderiving() && ignore_derivatives() do + if alg.verbosity > 0 && is_degenerate_spectrum(S) + svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S)) + @warn("degenerate singular values detected: ", svals) + end + end + + P_left, P_right = contract_projectors(U, S, V, R_left, L_right) + return (P_left, P_right), (; err, U, S, V) +end diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl new file mode 100644 index 00000000..94f4846e --- /dev/null +++ b/src/algorithms/ctmrg/sequential.jl @@ -0,0 +1,134 @@ +""" + SequentialCTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, verbosity=0, + projector_alg=typeof(Defaults.projector_alg), + svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) + +CTMRG algorithm where the expansions and renormalization is performed sequentially +column-wise. This is implemented as a growing and projecting step to the left, followed by +a clockwise rotation (performed four times). The projectors are computed using +`projector_alg` from `svd_alg` SVDs where the truncation scheme is set via `trscheme`. +""" +struct SequentialCTMRG <: CTMRGAlgorithm + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::ProjectorAlgorithm +end +function SequentialCTMRG(; + tol=Defaults.ctmrg_tol, + maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, + verbosity=2, + projector_alg=Defaults.projector_alg_type, + svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, +) + return SequentialCTMRG( + tol, maxiter, miniter, verbosity, projector_alg(; svd_alg, trscheme, verbosity) + ) +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) + ϵ = max(ϵ, info.err) + end + state = rotate_north(state, EAST) + envs = rotate_north(envs, EAST) + end + + return envs, (; err=ϵ) +end + +""" + sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm) + sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm) + +Compute CTMRG projectors in the `:sequential` scheme either for an entire column `col` or +for a specific `coordinate` (where `dir=WEST` is already implied in the `:sequential` scheme). +""" +function sequential_projectors( + col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm +) + # SVD half-infinite environment column-wise + ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs, 2))) + coordinates = eachcoordinate(envs)[:, col] + projectors = dtmap(coordinates) do (r, c) + trscheme = truncation_scheme(alg, envs.edges[WEST, _prev(r, size(envs, 2)), c]) + proj, info = sequential_projectors( + (WEST, r, c), state, envs, @set(alg.trscheme = trscheme) + ) + ϵ[r] = info.err / norm(info.S) + return proj + end + + return (map(first, projectors), map(last, projectors)), (; err=maximum(copy(ϵ))) +end +function sequential_projectors( + coordinate::NTuple{3,Int}, + state::InfinitePEPS, + envs::CTMRGEnv, + alg::HalfInfiniteProjector, +) + _, r, c = coordinate + r′ = _prev(r, size(envs, 2)) + Q1 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r, c)), SOUTHWEST) + Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST) + return compute_projector((Q1, Q2), coordinate, alg) +end +function sequential_projectors( + coordinate::NTuple{3,Int}, + state::InfinitePEPS, + envs::CTMRGEnv, + alg::FullInfiniteProjector, +) + rowsize, colsize = size(envs)[2:3] + coordinate_nw = _next_coordinate(coordinate, rowsize, colsize) + coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize) + coordinate_se = _next_coordinate(coordinate_ne, rowsize, colsize) + ec = ( + TensorMap(EnlargedCorner(state, envs, coordinate_se), SOUTHEAST), + TensorMap(EnlargedCorner(state, envs, coordinate), SOUTHWEST), + TensorMap(EnlargedCorner(state, envs, coordinate_nw), NORTHWEST), + TensorMap(EnlargedCorner(state, envs, coordinate_ne), NORTHEAST), + ) + return compute_projector(ec, coordinate, alg) +end + +""" + renormalize_sequentially(col::Int, projectors, state, envs) + +Renormalize one column of the CTMRG environment. +""" +function renormalize_sequentially(col::Int, projectors, state, envs) + corners = Zygote.Buffer(envs.corners) + edges = Zygote.Buffer(envs.edges) + + for (dir, r, c) in eachcoordinate(state, 1:4) + (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue + corners[dir, r, c] = envs.corners[dir, r, c] + end + for (dir, r, c) in eachcoordinate(state, 1:4) + (c == col && dir == WEST) && continue + edges[dir, r, c] = envs.edges[dir, r, c] + end + + # Apply projectors to renormalize corners and edge + for row in axes(envs.corners, 2) + C_southwest = renormalize_bottom_corner((row, col), envs, projectors) + corners[SOUTHWEST, row, col] = C_southwest / norm(C_southwest) + + C_northwest = renormalize_top_corner((row, col), envs, projectors) + corners[NORTHWEST, row, col] = C_northwest / norm(C_northwest) + + E_west = renormalize_west_edge((row, col), envs, projectors, state) + edges[WEST, row, col] = E_west / norm(E_west) + end + + return CTMRGEnv(copy(corners), copy(edges)) +end diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl new file mode 100644 index 00000000..fcae797b --- /dev/null +++ b/src/algorithms/ctmrg/simultaneous.jl @@ -0,0 +1,138 @@ +""" + SimultaneousCTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, verbosity=0, + projector_alg=Defaults.projector_alg, + svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation()) + +CTMRG algorithm where all sides are grown and renormalized at the same time. In particular, +the projectors are applied to the corners from two sides simultaneously. The projectors are +computed using `projector_alg` from `svd_alg` SVDs where the truncation scheme is set via +`trscheme`. +""" +struct SimultaneousCTMRG <: CTMRGAlgorithm + tol::Float64 + maxiter::Int + miniter::Int + verbosity::Int + projector_alg::ProjectorAlgorithm +end +function SimultaneousCTMRG(; + tol=Defaults.ctmrg_tol, + maxiter=Defaults.ctmrg_maxiter, + miniter=Defaults.ctmrg_miniter, + verbosity=2, + projector_alg=Defaults.projector_alg_type, + svd_alg=Defaults.svd_alg, + trscheme=Defaults.trscheme, +) + return SimultaneousCTMRG( + tol, maxiter, miniter, verbosity, projector_alg(; svd_alg, trscheme, verbosity) + ) +end + +function ctmrg_iteration(state, envs::CTMRGEnv, alg::SimultaneousCTMRG) + enlarged_corners = dtmap(eachcoordinate(state, 1:4)) do idx + return TensorMap(EnlargedCorner(state, envs, idx), idx[1]) + end # expand environment + projectors, info = simultaneous_projectors(enlarged_corners, envs, alg.projector_alg) # compute projectors on all coordinates + envs′ = renormalize_simultaneously(enlarged_corners, projectors, state, envs) # renormalize enlarged corners + return envs′, info +end + +# Pre-allocate U, S, and V tensor as Zygote buffers to make it differentiable +function _prealloc_svd(edges::Array{E,N}, ::HalfInfiniteProjector) where {E,N} + Sc = scalartype(E) + U = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, space(e)), edges)) + V = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, domain(e), codomain(e)), edges)) + S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc))) # Corner type but with real numbers + return U, S, V +end +function _prealloc_svd(edges::Array{E,N}, ::FullInfiniteProjector) where {E,N} + Sc = scalartype(E) + Rspace(x) = fuse(codomain(x)) + U = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, Rspace(e), domain(e)), edges)) + V = Zygote.Buffer(map(e -> TensorMap(zeros, Sc, domain(e), Rspace(e)), edges)) + S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc))) # Corner type but with real numbers + return U, S, V +end + +""" + simultaneous_projectors(enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm) + simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgorithm) + +Compute CTMRG projectors in the `:simultaneous` scheme either for all provided +enlarged corners or on a specific `coordinate`. +""" +function simultaneous_projectors( + enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm +) where {E} + U, S, V = _prealloc_svd(envs.edges, alg) + ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs))) + + projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate + coordinate′ = _next_coordinate(coordinate, size(envs)[2:3]...) + trscheme = truncation_scheme(alg, envs.edges[coordinate[1], coordinate′[2:3]...]) + proj, info = simultaneous_projectors( + coordinate, enlarged_corners, @set(alg.trscheme = trscheme) + ) + U[coordinate...] = info.U + S[coordinate...] = info.S + V[coordinate...] = info.V + ϵ[coordinate...] = info.err / norm(info.S) + return proj + end + + P_left = map(first, projectors) + P_right = map(last, projectors) + return (P_left, P_right), (; err=maximum(copy(ϵ)), U=copy(U), S=copy(S), V=copy(V)) +end +function simultaneous_projectors( + coordinate, enlarged_corners::Array{E,3}, alg::HalfInfiniteProjector +) where {E} + coordinate′ = _next_coordinate(coordinate, size(enlarged_corners)[2:3]...) + ec = (enlarged_corners[coordinate...], enlarged_corners[coordinate′...]) + return compute_projector(ec, coordinate, alg) +end +function simultaneous_projectors( + coordinate, enlarged_corners::Array{E,3}, alg::FullInfiniteProjector +) where {E} + rowsize, colsize = size(enlarged_corners)[2:3] + coordinate2 = _next_coordinate(coordinate, rowsize, colsize) + coordinate3 = _next_coordinate(coordinate2, rowsize, colsize) + coordinate4 = _next_coordinate(coordinate3, rowsize, colsize) + ec = ( + enlarged_corners[coordinate4...], + enlarged_corners[coordinate...], + enlarged_corners[coordinate2...], + enlarged_corners[coordinate3...], + ) + return compute_projector(ec, coordinate, alg) +end + +""" + renormalize_simultaneously(enlarged_corners, projectors, state, envs) + +Renormalize all enlarged corners and edges simultaneously. +""" +function renormalize_simultaneously(enlarged_corners, projectors, state, envs) + P_left, P_right = projectors + coordinates = eachcoordinate(envs, 1:4) + corners_edges = dtmap(coordinates) do (dir, r, c) + if dir == NORTH + corner = renormalize_northwest_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_north_edge((r, c), envs, P_left, P_right, state) + elseif dir == EAST + corner = renormalize_northeast_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_east_edge((r, c), envs, P_left, P_right, state) + elseif dir == SOUTH + corner = renormalize_southeast_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_south_edge((r, c), envs, P_left, P_right, state) + elseif dir == WEST + corner = renormalize_southwest_corner((r, c), enlarged_corners, P_left, P_right) + edge = renormalize_west_edge((r, c), envs, P_left, P_right, state) + end + return corner / norm(corner), edge / norm(edge) + end + + return CTMRGEnv(map(first, corners_edges), map(last, corners_edges)) +end diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 9243067e..3e8d388e 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -97,9 +97,24 @@ function renormalize_southwest_corner(ec::EnlargedCorner, P_left, P_right) ) end -# ------------------ -# Sparse environment -# ------------------ +# Wrapper around half_infinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) +function half_infinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) + return HalfInfiniteEnv(ec_1, ec_2) +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_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 + ) + return P_left, P_right +end + +# -------------------------------- +# Sparse half-infinite environment +# -------------------------------- """ struct HalfInfiniteEnv{C,E,A,A′} @@ -114,8 +129,8 @@ struct HalfInfiniteEnv{C,E,A,A′} # TODO: subtype as AbstractTensorMap once Te E_3::E E_4::E ket_1::A - bra_1::A′ ket_2::A + bra_1::A′ bra_2::A′ end @@ -128,8 +143,10 @@ function HalfInfiniteEnv(quadrant1::EnlargedCorner, quadrant2::EnlargedCorner) quadrant1.E_2, quadrant2.E_1, quadrant2.E_2, - quadrant1.ket_bra, - quadrant2.ket_bra, + quadrant1.ket, + quadrant2.ket, + quadrant1.bra, + quadrant2.bra, ) end @@ -139,7 +156,7 @@ end Instantiate half-infinite environment as `TensorMap` explicitly. """ function TensorKit.TensorMap(env::HalfInfiniteEnv) # Dense operator - return halfinfinite_environment( + return half_infinite_environment( env.C_1, env.C_2, env.E_1, @@ -161,7 +178,7 @@ Contract half-infinite environment with a vector `x`, such that the environment linear map or adjoint linear map on `x` if `Val(true)` or `Val(false)` is passed, respectively. """ function (env::HalfInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x - return halfinfinite_environment( + return half_infinite_environment( env.C_1, env.C_2, env.E_1, @@ -176,7 +193,7 @@ function (env::HalfInfiniteEnv)(x, ::Val{false}) # Linear map: env() * x ) end function (env::HalfInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * x - return halfinfinite_environment( + return half_infinite_environment( x, env.C_1, env.C_2, @@ -191,26 +208,7 @@ function (env::HalfInfiniteEnv)(x, ::Val{true}) # Adjoint linear map: env()' * ) end -# Wrapper around halfinfinite_environment contraction using EnlargedCorners (used in ctmrg_projectors) -function halfinfinite_environment(ec_1::EnlargedCorner, ec_2::EnlargedCorner) - return HalfInfiniteEnv( - ec_1.C, - ec_2.C, - ec_1.E_1, - ec_1.E_2, - ec_2.E_1, - ec_2.E_2, - ec_1.ket, - ec_2.ket, - ec_1.bra, - ec_2.bra, - ) -end - -# ----------------------------------------------------- # AbstractTensorMap subtyping and IterSVD compatibility -# ----------------------------------------------------- - function TensorKit.domain(env::HalfInfiniteEnv) return domain(env.E_4) * domain(env.ket_2)[3] * domain(env.bra_2)[3]' end diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 8b927f93..7be93c92 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -77,31 +77,31 @@ function LinSolver(; end """ - PEPSOptimize{G}(; boundary_alg=CTMRG(), optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer - reuse_env::Bool=true, gradient_alg::G=LinSolver()) + PEPSOptimize{G}(; boundary_alg=Defaults.ctmrg_alg, optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer + reuse_env::Bool=true, gradient_alg::G=Defaults.gradient_alg) 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 step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm. """ struct PEPSOptimize{G} - boundary_alg::CTMRG + boundary_alg::CTMRGAlgorithm optimizer::OptimKit.OptimizationAlgorithm reuse_env::Bool gradient_alg::G function PEPSOptimize( # Inner constructor to prohibit illegal setting combinations - boundary_alg::CTMRG{S}, + boundary_alg::CTMRGAlgorithm, optimizer, reuse_env, gradient_alg::G, - ) where {S,G} + ) where {G} if gradient_alg isa GradMode - if S === :sequential && iterscheme(gradient_alg) === :fixed + if boundary_alg isa SequentialCTMRG && iterscheme(gradient_alg) === :fixed throw(ArgumentError(":sequential and :fixed are not compatible")) end end @@ -109,7 +109,7 @@ struct PEPSOptimize{G} end end function PEPSOptimize(; - boundary_alg=CTMRG(), + boundary_alg=Defaults.ctmrg_alg, optimizer=Defaults.optimizer, reuse_env=Defaults.reuse_env, gradient_alg=Defaults.gradient_alg, @@ -216,26 +216,24 @@ Evaluating the gradient of the cost function for CTMRG: function _rrule( gradmode::GradMode{:diffgauge}, - ::RuleConfig, + config::RuleConfig, ::typeof(MPSKit.leading_boundary), envinit, state, - alg::CTMRG, + alg::CTMRGAlgorithm, ) - envs = leading_boundary(envinit, state, alg) #TODO: fixed space for unit cells + envs = leading_boundary(envinit, state, alg) function leading_boundary_diffgauge_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) # find partial gradients of gauge_fixed single CTMRG iteration - # TODO: make this rrule_via_ad so it's zygote-agnostic - _, env_vjp = pullback(state, envs) do A, x - return gauge_fix(x, ctmrg_iter(A, x, alg)[1])[1] - end + f(A, x) = gauge_fix(x, ctmrg_iteration(A, x, alg)[1])[1] + _, env_vjp = rrule_via_ad(config, f, state, envs) # evaluate the geometric sum - ∂f∂A(x)::typeof(state) = env_vjp(x)[1] - ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂f∂A(x)::typeof(state) = env_vjp(x)[2] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[3] ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() @@ -247,44 +245,40 @@ end # Here f is differentiated from an pre-computed SVD with fixed U, S and V function _rrule( gradmode::GradMode{:fixed}, - ::RuleConfig, + config::RuleConfig, ::typeof(MPSKit.leading_boundary), envinit, state, - alg::CTMRG{C}, -) where {C} - @assert C === :simultaneous + alg::SimultaneousCTMRG, +) @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg) envs = leading_boundary(envinit, state, alg) - envsconv, info = ctmrg_iter(state, envs, alg) - envsfix, signs = gauge_fix(envs, envsconv) + envsconv, info = ctmrg_iteration(state, envs, alg) + envs_fixed, signs = gauge_fix(envs, envsconv) # Fix SVD Ufix, Vfix = fix_relative_phases(info.U, info.V, signs) svd_alg_fixed = SVDAdjoint(; fwd_alg=FixedSVD(Ufix, info.S, Vfix), rrule_alg=alg.projector_alg.svd_alg.rrule_alg ) - alg_fixed = CTMRG(; - svd_alg=svd_alg_fixed, trscheme=notrunc(), ctmrgscheme=:simultaneous - ) + alg_fixed = @set alg.projector_alg.svd_alg = svd_alg_fixed + alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc() function leading_boundary_fixed_pullback(Δenvs′) Δenvs = unthunk(Δenvs′) - _, env_vjp = pullback(state, envsfix) do A, x - e, = ctmrg_iter(A, x, alg_fixed) - return fix_global_phases(x, e) - end + f(A, x) = fix_global_phases(x, ctmrg_iteration(A, x, alg_fixed)[1]) + _, env_vjp = rrule_via_ad(config, f, state, envs_fixed) # evaluate the geometric sum - ∂f∂A(x)::typeof(state) = env_vjp(x)[1] - ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂f∂A(x)::typeof(state) = env_vjp(x)[2] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[3] ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() end - return envsfix, leading_boundary_fixed_pullback + return envs_fixed, leading_boundary_fixed_pullback end @doc """ diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 7bca4b90..97d890de 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -355,6 +355,24 @@ function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnv, name::Symbol) end end +Base.size(env::CTMRGEnv, args...) = size(env.corners, args...) +Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) +function eachcoordinate(x::CTMRGEnv) + return collect(Iterators.product(axes(x, 2), axes(x, 3))) +end +function eachcoordinate(x::CTMRGEnv, dirs) + return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) +end +Base.real(env::CTMRGEnv) = CTMRGEnv(real.(env.corners), real.(env.edges)) +Base.complex(env::CTMRGEnv) = CTMRGEnv(complex.(env.corners), complex.(env.edges)) + +# In-place update of environment +function update!(env::CTMRGEnv{C,T}, env´::CTMRGEnv{C,T}) where {C,T} + env.corners .= env´.corners + env.edges .= env´.edges + return env +end + # Rotate corners & edges counter-clockwise function Base.rotl90(env::CTMRGEnv{C,T}) where {C,T} # Initialize rotated corners & edges with rotated sizes @@ -400,24 +418,6 @@ function Base.rot180(env::CTMRGEnv{C,T}) where {C,T} return CTMRGEnv(copy(corners′), copy(edges′)) end -Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) -function eachcoordinate(x::CTMRGEnv) - return collect(Iterators.product(axes(x, 2), axes(x, 3))) -end -function eachcoordinate(x::CTMRGEnv, dirs) - return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) -end - -Base.real(env::CTMRGEnv) = CTMRGEnv(real.(env.corners), real.(env.edges)) -Base.complex(env::CTMRGEnv) = CTMRGEnv(complex.(env.corners), complex.(env.edges)) - -# In-place update of environment -function update!(env::CTMRGEnv{C,T}, env´::CTMRGEnv{C,T}) where {C,T} - env.corners .= env´.corners - env.edges .= env´.edges - return env -end - # Functions used for FP differentiation and by KrylovKit.linsolve function Base.:+(e₁::CTMRGEnv, e₂::CTMRGEnv) return CTMRGEnv(e₁.corners + e₂.corners, e₁.edges + e₂.edges) diff --git a/src/utility/util.jl b/src/utility/util.jl index 45182f1d..075f39ea 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -1,7 +1,32 @@ -# Get next and previous directional CTM enviroment index, respecting periodicity +# Get next and previous directional CTMRG environment index, respecting periodicity _next(i, total) = mod1(i + 1, total) _prev(i, total) = mod1(i - 1, total) +# Get next and previous coordinate (direction, row, column), given a direction and going around the environment clockwise +function _next_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_next(dir, 4), row, _next(col, colsize)) + elseif dir == 2 + return (_next(dir, 4), _next(row, rowsize), col) + elseif dir == 3 + return (_next(dir, 4), row, _prev(col, colsize)) + elseif dir == 4 + return (_next(dir, 4), _prev(row, rowsize), col) + end +end +function _prev_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_prev(dir, 4), _next(row, rowsize), col) + elseif dir == 2 + return (_prev(dir, 4), row, _prev(col, colsize)) + elseif dir == 3 + return (_prev(dir, 4), _prev(row, rowsize), col) + elseif dir == 4 + return (_prev(dir, 4), row, _next(col, colsize)) + end +end + +# iterator over each coordinates """ eachcoordinate(x, dirs=1:4) @@ -80,28 +105,6 @@ function is_degenerate_spectrum( return false end -""" - projector_type(T::DataType, size) - projector_type(edges::Array{<:AbstractTensorMap}) - -Create two arrays of specified `size` that contain undefined tensors representing -left and right acting projectors, respectively. The projector types are inferred -from the TensorMap type `T` which avoids having to recompute transpose tensors. -Alternatively, supply an array of edge tensors from which left and right projectors -are intialized explicitly with zeros. -""" -function projector_type(T::DataType, size) - P_left = Array{T,length(size)}(undef, size) - Prtype = tensormaptype(spacetype(T), numin(T), numout(T), storagetype(T)) - P_right = Array{Prtype,length(size)}(undef, size) - return P_left, P_right -end -function projector_type(edges::Array{<:AbstractTensorMap}) - P_left = map(e -> TensorMap(zeros, scalartype(e), space(e)), edges) - P_right = map(e -> TensorMap(zeros, scalartype(e), domain(e), codomain(e)), edges) - return P_left, P_right -end - # There are no rrules for rotl90 and rotr90 in ChainRules.jl function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix) function rotl90_pullback(x) diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl index a161c073..2ec6265c 100644 --- a/test/boundarymps/vumps.jl +++ b/test/boundarymps/vumps.jl @@ -17,7 +17,9 @@ const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitia mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(sum(expectation_value(mps, T))) - ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) + ctm = leading_boundary( + CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1) + ) N´ = abs(norm(psi, ctm)) @test N ≈ N´ atol = 1e-3 @@ -31,7 +33,9 @@ end mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(prod(expectation_value(mps, T))) - ctm = leading_boundary(CTMRGEnv(psi, ComplexSpace(20)), psi, CTMRG(; verbosity=1)) + ctm = leading_boundary( + CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1) + ) N´ = abs(norm(psi, ctm)) @test N ≈ N´ rtol = 1e-2 diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl index 3f005986..208fdece 100644 --- a/test/ctmrg/fixed_iterscheme.jl +++ b/test/ctmrg/fixed_iterscheme.jl @@ -5,7 +5,7 @@ using TensorKit, KrylovKit using PEPSKit using PEPSKit: FixedSVD, - ctmrg_iter, + ctmrg_iteration, gauge_fix, fix_relative_phases, fix_global_phases, @@ -15,14 +15,17 @@ using PEPSKit: χbond = 2 χenv = 16 svd_algs = [SVDAdjoint(; fwd_alg=TensorKit.SDD()), SVDAdjoint(; fwd_alg=IterSVD())] +projector_algs = [HalfInfiniteProjector] #, FullInfiniteProjector] unitcells = [(1, 1), (3, 4)] +atol = 1e-5 # test for element-wise convergence after application of fixed step -@testset "$unitcell unit cell with $(typeof(svd_alg.fwd_alg))" for (unitcell, svd_alg) in - Iterators.product( - unitcells, svd_algs +@testset "$unitcell unit cell with $(typeof(svd_alg.fwd_alg)) and $projector_alg" for ( + unitcell, svd_alg, projector_alg +) in Iterators.product( + unitcells, svd_algs, projector_algs ) - ctm_alg = CTMRG(; svd_alg) + ctm_alg = SimultaneousCTMRG(; svd_alg, projector_alg) # initialize states Random.seed!(2394823842) @@ -30,27 +33,29 @@ unitcells = [(1, 1), (3, 4)] env_conv1 = leading_boundary(CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg) # do extra iteration to get SVD - env_conv2, info = ctmrg_iter(psi, env_conv1, ctm_alg) + env_conv2, info = ctmrg_iteration(psi, env_conv1, ctm_alg) env_fix, signs = gauge_fix(env_conv1, env_conv2) - @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = atol # fix gauge of SVD U_fix, V_fix = fix_relative_phases(info.U, info.V, signs) svd_alg_fix = SVDAdjoint(; fwd_alg=FixedSVD(U_fix, info.S, V_fix)) - ctm_alg_fix = CTMRG(; svd_alg=svd_alg_fix, trscheme=notrunc()) + ctm_alg_fix = SimultaneousCTMRG(; + projector_alg, svd_alg=svd_alg_fix, trscheme=notrunc() + ) # do iteration with FixedSVD - env_fixedsvd, = ctmrg_iter(psi, env_conv1, ctm_alg_fix) + env_fixedsvd, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix) env_fixedsvd = fix_global_phases(env_conv1, env_fixedsvd) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = atol end @testset "Element-wise consistency of TensorKit.SDD and IterSVD" begin - ctm_alg_iter = CTMRG(; + ctm_alg_iter = SimultaneousCTMRG(; maxiter=200, svd_alg=SVDAdjoint(; fwd_alg=IterSVD(; alg=GKL(; tol=1e-14, krylovdim=χenv + 10))), ) - ctm_alg_full = CTMRG(; svd_alg=SVDAdjoint(; fwd_alg=TensorKit.SDD())) + ctm_alg_full = SimultaneousCTMRG(; svd_alg=SVDAdjoint(; fwd_alg=TensorKit.SDD())) # initialize states Random.seed!(91283219347) @@ -59,34 +64,33 @@ end env_conv1 = leading_boundary(env_init, psi, ctm_alg_iter) # do extra iteration to get SVD - env_conv2_iter, info_iter = ctmrg_iter(psi, env_conv1, ctm_alg_iter) + env_conv2_iter, info_iter = ctmrg_iteration(psi, env_conv1, ctm_alg_iter) env_fix_iter, signs_iter = gauge_fix(env_conv1, env_conv2_iter) - @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = atol - env_conv2_full, info_full = ctmrg_iter(psi, env_conv1, ctm_alg_full) + env_conv2_full, info_full = ctmrg_iteration(psi, env_conv1, ctm_alg_full) env_fix_full, signs_full = gauge_fix(env_conv1, env_conv2_full) - @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = atol # fix gauge of SVD U_fix_iter, V_fix_iter = fix_relative_phases(info_iter.U, info_iter.V, signs_iter) svd_alg_fix_iter = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_iter, info_iter.S, V_fix_iter)) - ctm_alg_fix_iter = CTMRG(; svd_alg=svd_alg_fix_iter, trscheme=notrunc()) + ctm_alg_fix_iter = SimultaneousCTMRG(; svd_alg=svd_alg_fix_iter, trscheme=notrunc()) U_fix_full, V_fix_full = fix_relative_phases(info_full.U, info_full.V, signs_full) svd_alg_fix_full = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_full, info_full.S, V_fix_full)) - ctm_alg_fix_full = CTMRG(; svd_alg=svd_alg_fix_full, trscheme=notrunc()) + ctm_alg_fix_full = SimultaneousCTMRG(; svd_alg=svd_alg_fix_full, trscheme=notrunc()) # do iteration with FixedSVD - env_fixedsvd_iter, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_iter) + env_fixedsvd_iter, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_iter) env_fixedsvd_iter = fix_global_phases(env_conv1, env_fixedsvd_iter) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = 1e-6 # This doesn't work for x₀ = rand(size(b, 1))? + @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = atol # This doesn't work for x₀ = rand(size(b, 1))? - env_fixedsvd_full, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_full) + env_fixedsvd_full, = ctmrg_iteration(psi, env_conv1, ctm_alg_fix_full) env_fixedsvd_full = fix_global_phases(env_conv1, env_fixedsvd_full) - @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = 1e-6 + @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = atol # check matching decompositions - atol = 1e-12 decomposition_check = all( zip(info_iter.U, info_iter.S, info_iter.V, info_full.U, info_full.S, info_full.V), ) do (U_iter, S_iter, V_iter, U_full, S_full, V_full) diff --git a/test/ctmrg/ctmrgschemes.jl b/test/ctmrg/flavors.jl similarity index 76% rename from test/ctmrg/ctmrgschemes.jl rename to test/ctmrg/flavors.jl index 406d4c1f..5d8b440c 100644 --- a/test/ctmrg/ctmrgschemes.jl +++ b/test/ctmrg/flavors.jl @@ -7,11 +7,15 @@ using PEPSKit # initialize parameters χbond = 2 χenv = 16 -ctm_alg_sequential = CTMRG(; ctmrgscheme=:sequential) -ctm_alg_simultaneous = CTMRG(; ctmrgscheme=:simultaneous) +ctm_alg_sequential = SequentialCTMRG() +ctm_alg_simultaneous = SimultaneousCTMRG() unitcells = [(1, 1), (3, 4)] +projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] -@testset "$(unitcell) unit cell" for unitcell in unitcells +@testset "$(unitcell) unit cell with $projector_alg" for (unitcell, projector_alg) in + Iterators.product( + unitcells, projector_algs +) # compute environments Random.seed!(32350283290358) psi = InfinitePEPS(2, χbond; unitcell) @@ -52,13 +56,13 @@ unitcells = [(1, 1), (3, 4)] end # test fixedspace actually fixes space -@testset "Fixedspace truncation ($scheme)" for scheme in [:sequential, :simultaneous] - ctm_alg = CTMRG(; - tol=1e-6, - maxiter=1, - verbosity=0, - ctmrgscheme=scheme, - trscheme=FixedSpaceTruncation(), +@testset "Fixedspace truncation using $ctmrg_alg and $projector_alg" for ( + ctmrg_alg, projector_alg +) in Iterators.product( + [SequentialCTMRG, SimultaneousCTMRG], projector_algs +) + ctm_alg = ctmrg_alg(; + tol=1e-6, maxiter=1, verbosity=0, trscheme=FixedSpaceTruncation(), projector_alg ) Ds = fill(2, 3, 3) χs = [16 17 18; 15 20 21; 14 19 22] diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl index ef36943c..048d875d 100644 --- a/test/ctmrg/gaugefix.jl +++ b/test/ctmrg/gaugefix.jl @@ -3,81 +3,49 @@ using Random using PEPSKit using TensorKit -using PEPSKit: ctmrg_iter, gauge_fix, calc_elementwise_convergence +using PEPSKit: ctmrg_iteration, gauge_fix, calc_elementwise_convergence +spacetypes = [ComplexSpace, Z2Space] scalartypes = [Float64, ComplexF64] unitcells = [(1, 1), (2, 2), (3, 2)] -maxiter = 200 -schemes = [:simultaneous, :sequential] +ctmrg_algs = [SequentialCTMRG, SimultaneousCTMRG] +projector_algs = [HalfInfiniteProjector, FullInfiniteProjector] χ = 6 atol = 1e-4 -function _make_symmetric!(psi) - if ==(size(psi)...) - return symmetrize!(psi, RotateReflect()) - else - return symmetrize!(symmetrize!(psi, ReflectDepth()), ReflectWidth()) - end +function _pre_converge_env( + ::Type{T}, physical_space, peps_space, ctm_space, unitcell; seed=12345 +) where {T} + Random.seed!(seed) # Seed RNG to make random environment consistent + psi = InfinitePEPS(rand, T, physical_space, peps_space; unitcell) + env₀ = CTMRGEnv(psi, ctm_space) + env_conv = leading_boundary(env₀, psi, SequentialCTMRG()) + return env_conv, psi end -# If I can't make the rng seed behave, I'll just randomly define a peps somehow -function _semi_random_peps!(psi::InfinitePEPS) - i = 0 - A′ = map(psi.A) do a - for (_, b) in blocks(a) - l = length(b) - b .= reshape(collect((1:l) .+ i), size(b)) - i += l - end - return a +# pre-converge CTMRG environments with given spacetype, scalartype and unit cell +preconv = Dict() +for (S, T, unitcell) in Iterators.product(spacetypes, scalartypes, unitcells) + if S == ComplexSpace + result = _pre_converge_env(T, S(2), S(2), S(χ), unitcell) + elseif S == Z2Space + result = _pre_converge_env( + T, S(0 => 1, 1 => 1), S(0 => 1, 1 => 1), S(0 => χ ÷ 2, 1 => χ ÷ 2), unitcell + ) end - return InfinitePEPS(A′) + push!(preconv, (S, T, unitcell) => result) end -@testset "Trivial symmetry ($T) - ($unitcell) - ($ctmrgscheme)" for ( - T, unitcell, ctmrgscheme +@testset "($S) - ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg)" for ( + S, T, unitcell, ctmrg_alg, projector_alg ) in Iterators.product( - scalartypes, unitcells, schemes -) - physical_space = ComplexSpace(2) - peps_space = ComplexSpace(2) - ctm_space = ComplexSpace(χ) - - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) - - Random.seed!(987654321) # Seed RNG to make random environment consistent - ctm = CTMRGEnv(psi, ctm_space) - - alg = CTMRG(; maxiter, ctmrgscheme) - - ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iter(psi, ctm, alg) - ctm_fixed, = gauge_fix(ctm, ctm2) - @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol -end - -@testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrgscheme)" for (T, unitcell, ctmrgscheme) in - Iterators.product( - scalartypes, unitcells, schemes + spacetypes, scalartypes, unitcells, ctmrg_algs, projector_algs ) - physical_space = Z2Space(0 => 1, 1 => 1) - peps_space = Z2Space(0 => 1, 1 => 1) - ctm_space = Z2Space(0 => χ ÷ 2, 1 => χ ÷ 2) - - psi = InfinitePEPS(undef, T, physical_space, peps_space; unitcell) - _semi_random_peps!(psi) - _make_symmetric!(psi) - - Random.seed!(987654321) # Seed RNG to make random environment consistent - psi = InfinitePEPS(physical_space, peps_space; unitcell) - ctm = CTMRGEnv(psi, ctm_space) - - alg = CTMRG(; maxiter, ctmrgscheme) - - ctm = leading_boundary(ctm, psi, alg) - ctm2, = ctmrg_iter(psi, ctm, alg) - ctm_fixed, = gauge_fix(ctm, ctm2) - @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = atol + alg = ctmrg_alg(; projector_alg) + env_pre, psi = preconv[(S, T, unitcell)] + env_pre + env = leading_boundary(env_pre, psi, alg) + env′, = ctmrg_iteration(psi, env, alg) + env_fixed, = gauge_fix(env, env′) + @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol end diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 2c18ef90..648d48bd 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -17,9 +17,12 @@ models = [heisenberg_XYZ(InfiniteSquare()), pwave_superconductor(InfiniteSquare( names = ["Heisenberg", "p-wave superconductor"] gradtol = 1e-4 -boundary_algs = [ - CTMRG(; verbosity=0, ctmrgscheme=:simultaneous), - CTMRG(; verbosity=0, ctmrgscheme=:sequential), +ctmrg_algs = [ + [ + SimultaneousCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; verbosity=0, projector_alg=FullInfiniteProjector), + ], + [SequentialCTMRG(; verbosity=0, projector_alg=HalfInfiniteProjector)], ] gradmodes = [ [ @@ -50,14 +53,15 @@ steps = -0.01:0.005:0.01 Vspace = Pspaces[i] Espace = Espaces[i] gms = gradmodes[i] - boundary_alg = boundary_algs[i] + calgs = ctmrg_algs[i] psi_init = InfinitePEPS(Pspace, Vspace, Vspace) - @testset "$alg_rrule" for alg_rrule in gms - @info "optimtest of $alg_rrule on $(names[i])" + @testset "$ctmrg_alg and $alg_rrule" for (ctmrg_alg, alg_rrule) in + Iterators.product(calgs, gms) + @info "optimtest of $ctmrg_alg and $alg_rrule on $(names[i])" Random.seed!(42039482030) dir = InfinitePEPS(Pspace, Vspace, Vspace) psi = InfinitePEPS(Pspace, Vspace, Vspace) - env = leading_boundary(CTMRGEnv(psi, Espace), psi, boundary_alg) + env = leading_boundary(CTMRGEnv(psi, Espace), psi, ctmrg_alg) alphas, fs, dfs1, dfs2 = OptimKit.optimtest( (psi, env), dir; @@ -66,9 +70,7 @@ steps = -0.01:0.005:0.01 inner=PEPSKit.real_inner, ) do (peps, envs) E, g = Zygote.withgradient(peps) do psi - envs2 = PEPSKit.hook_pullback( - leading_boundary, envs, psi, boundary_alg; alg_rrule - ) + envs2 = PEPSKit.hook_pullback(leading_boundary, envs, psi, ctmrg_alg; alg_rrule) return costfun(psi, envs2, models[i]) end diff --git a/test/ctmrg/gradparts.jl b/test/ctmrg/gradparts.jl deleted file mode 100644 index 846ea90a..00000000 --- a/test/ctmrg/gradparts.jl +++ /dev/null @@ -1,76 +0,0 @@ -using Test -using Random -using PEPSKit -using TensorKit -using PEPSKit: - NORTH, - SOUTH, - WEST, - EAST, - NORTHWEST, - NORTHEAST, - SOUTHEAST, - SOUTHWEST, - rotate_north, - left_move, - ctmrg_iter - -using Zygote -using ChainRulesCore -using ChainRulesTestUtils - -include(joinpath("..", "utility.jl")) - -## Test spaces, tested functions and CTMRG algorithm -# -------------------------------------------------- -χbond = 2 -χenv = 4 -Pspaces = [ComplexSpace(2), Vect[FermionParity](0 => 1, 1 => 1)] -Vspaces = [ComplexSpace(χbond), Vect[FermionParity](0 => χbond / 2, 1 => χbond / 2)] -Espaces = [ComplexSpace(χenv), Vect[FermionParity](0 => χenv / 2, 1 => χenv / 2)] -tol = 1e-10 -atol = 1e-6 -boundary_algs = [ - CTMRG(; tol, verbosity=0, ctmrgscheme=:simultaneous), - CTMRG(; tol, verbosity=0, ctmrgscheme=:sequential), -] - -## Gauge invariant function of the environment -# -------------------------------------------- -function rho(env) - @tensor ρ[-1 -2 -3 -4 -5 -6 -7 -8] := - env.edges[WEST, 1, 1][1 -1 -2; 4] * - env.corners[NORTHWEST, 1, 1][4; 5] * - env.edges[NORTH, 1, 1][5 -3 -4; 8] * - env.corners[NORTHEAST, 1, 1][8; 9] * - env.edges[EAST, 1, 1][9 -5 -6; 12] * - env.corners[SOUTHEAST, 1, 1][12; 13] * - env.edges[SOUTH, 1, 1][13 -7 -8; 16] * - env.corners[SOUTHWEST, 1, 1][16; 1] - return ρ -end - -## Tests -# ------ -@testset "Reverse rules for ctmrg_iter with spacetype $(Vspaces[i])" for i in - eachindex(Pspaces) - Random.seed!(42039482030) - psi = InfinitePEPS(Pspaces[i], Vspaces[i], Vspaces[i]) - env = CTMRGEnv(psi, Espaces[i]) - - @testset "$alg" for alg in boundary_algs - @info "$(typeof(alg)) on $(Vspaces[i])" - f(state, env) = rho(ctmrg_iter(state, env, alg)[1]) - - # use rrule testing functionality but compute rrule via Zygote - test_rrule( - Zygote.ZygoteRuleConfig(), - f, - psi, - env; - check_inferred=false, - atol, - rrule_f=rrule_via_ad, - ) - end -end diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 0ffbdc09..126126f2 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -1,22 +1,35 @@ using Test using Random using PEPSKit -using PEPSKit: _prev, _next, ctmrg_iter +using PEPSKit: _prev, _next, ctmrg_iteration using TensorKit # settings Random.seed!(91283219347) stype = ComplexF64 -ctm_alg = CTMRG() +ctm_algs = [ + SequentialCTMRG(; projector_alg=HalfInfiniteProjector), + SequentialCTMRG(; projector_alg=FullInfiniteProjector), + SimultaneousCTMRG(; projector_alg=HalfInfiniteProjector), + SimultaneousCTMRG(; projector_alg=FullInfiniteProjector), +] function test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) peps = InfinitePEPS(randn, stype, Pspaces, Nspaces, Espaces) env = CTMRGEnv(randn, stype, peps, chis_north, chis_east, chis_south, chis_west) # apply one CTMRG iteration with fixeds - env′, = ctmrg_iter(peps, env, ctm_alg) + env′, = ctmrg_iteration(peps, env, ctm_alg) # compute random expecation value to test matching bonds random_op = LocalOperator( @@ -40,7 +53,7 @@ function random_dualize!(M::AbstractMatrix{<:ElementarySpace}) return M end -@testset "Integer space specifiers" begin +@testset "Integer space specifiers with $ctm_alg" for ctm_alg in ctm_algs unitcell = (3, 3) Pspaces = rand(2:3, unitcell...) @@ -52,11 +65,19 @@ end chis_west = rand(5:10, unitcell...) test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) end -@testset "Random Cartesian spaces" begin +@testset "Random Cartesian spaces with $ctm_alg" for ctm_alg in ctm_algs unitcell = (3, 3) Pspaces = random_dualize!(ComplexSpace.(rand(2:3, unitcell...))) @@ -68,11 +89,19 @@ end chis_west = random_dualize!(ComplexSpace.(rand(5:10, unitcell...))) test_unitcell( - unitcell, Pspaces, Nspaces, Espaces, chis_north, chis_east, chis_south, chis_west + ctm_alg, + unitcell, + Pspaces, + Nspaces, + Espaces, + chis_north, + chis_east, + chis_south, + chis_west, ) end -@testset "Specific U1 spaces" begin +@testset "Specific U1 spaces with $ctm_alg" for ctm_alg in ctm_algs unitcell = (2, 2) PA = U1Space(-1 => 1, 0 => 1) @@ -84,5 +113,5 @@ end Nspaces = [Vpeps Vpeps'; Vpeps' Vpeps] chis = [Venv Venv; Venv Venv] - test_unitcell(unitcell, Pspaces, Nspaces, Nspaces, chis, chis, chis, chis) + test_unitcell(ctm_alg, unitcell, Pspaces, Nspaces, Nspaces, chis, chis, chis, chis) end diff --git a/test/heisenberg.jl b/test/heisenberg.jl index a8bbc426..5865ba46 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -9,7 +9,7 @@ using OptimKit # initialize parameters Dbond = 2 χenv = 16 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) ) @@ -82,9 +82,8 @@ 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)) + envs₀ = CTMRGEnv(rand, Float64, peps, Espace) + envs = leading_boundary(envs₀, peps, SimultaneousCTMRG()) # measure physical quantities e_site = costfun(peps, envs, ham) / (N1 * N2) diff --git a/test/j1j2_model.jl b/test/j1j2_model.jl index 45acd5b1..8dd88397 100644 --- a/test/j1j2_model.jl +++ b/test/j1j2_model.jl @@ -8,7 +8,7 @@ using OptimKit # initialize parameters χbond = 2 χenv = 12 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2), diff --git a/test/pwave.jl b/test/pwave.jl index d663ec0d..dcdecf7e 100644 --- a/test/pwave.jl +++ b/test/pwave.jl @@ -10,7 +10,7 @@ unitcell = (2, 2) H = pwave_superconductor(InfiniteSquare(unitcell...)) χbond = 2 χenv = 16 -ctm_alg = CTMRG(; maxiter=150) +ctm_alg = SimultaneousCTMRG(; maxiter=150) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=2) ) diff --git a/test/runtests.jl b/test/runtests.jl index b15ad1c6..c80a4b43 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,8 +26,8 @@ end @time @safetestset "Unit cells" begin include("ctmrg/unitcell.jl") end - @time @safetestset "CTMRG schemes" begin - include("ctmrg/ctmrgschemes.jl") + @time @safetestset "Flavors" begin + include("ctmrg/flavors.jl") end end if GROUP == "ALL" || GROUP == "BOUNDARYMPS" diff --git a/test/tf_ising.jl b/test/tf_ising.jl index 0b4d72f5..8b84696f 100644 --- a/test/tf_ising.jl +++ b/test/tf_ising.jl @@ -19,7 +19,7 @@ mˣ = 0.91 # initialize parameters χbond = 2 χenv = 16 -ctm_alg = CTMRG() +ctm_alg = SimultaneousCTMRG() opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2) )