diff --git a/Project.toml b/Project.toml index aa1fc857..903111cd 100644 --- a/Project.toml +++ b/Project.toml @@ -7,12 +7,14 @@ version = "0.2.1" Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502" MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab" OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" +Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 5a4fdb0a..4fe6c037 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -2,6 +2,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Base: @kwdef +using Parameters using Compat using Accessors using VectorInterface @@ -11,81 +12,44 @@ using LoggingExtras using MPSKit: loginit!, logiter!, logfinish!, logcancel! using MPSKitModels + include("utility/util.jl") include("utility/svd.jl") include("utility/rotations.jl") include("utility/diffset.jl") include("utility/hook_pullback.jl") include("utility/autoopt.jl") +include("utility/loginfo.jl") include("states/abstractpeps.jl") include("states/infinitepeps.jl") -include("operators/transferpeps.jl") -include("operators/infinitepepo.jl") -include("operators/transferpepo.jl") -include("operators/derivatives.jl") include("operators/localoperator.jl") include("operators/lattices/squarelattice.jl") include("operators/models.jl") include("environments/ctmrg_environments.jl") -include("environments/transferpeps_environments.jl") -include("environments/transferpepo_environments.jl") +include("environments/vumps_environments.jl") include("algorithms/contractions/localoperator.jl") include("algorithms/contractions/ctmrg_contractions.jl") +include("algorithms/contractions/vumps_contractions.jl") include("algorithms/ctmrg/ctmrg.jl") include("algorithms/ctmrg/gaugefix.jl") +include("algorithms/vumps/vumps.jl") + include("algorithms/toolbox.jl") include("algorithms/peps_opt.jl") include("utility/symmetrization.jl") -""" - module Defaults - const ctmrg_maxiter = 100 - const ctmrg_miniter = 4 - const ctmrg_tol = 1e-12 - const fpgrad_maxiter = 100 - const fpgrad_tol = 1e-6 - end - -Module containing default values that represent typical algorithm parameters. - -- `ctmrg_maxiter = 100`: Maximal number of CTMRG iterations per run -- `ctmrg_miniter = 4`: Minimal number of CTMRG carried out -- `ctmrg_tol = 1e-12`: Tolerance checking singular value and norm convergence -- `fpgrad_maxiter = 100`: Maximal number of iterations for computing the CTMRG fixed-point gradient -- `fpgrad_tol = 1e-6`: Convergence tolerance for the fixed-point gradient iteration -""" -module Defaults - using TensorKit, KrylovKit, OptimKit - using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint - 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 reuse_env = true - const trscheme = FixedSpaceTruncation() - const iterscheme = :fixed - const fwd_alg = TensorKit.SVD() - const rrule_alg = GMRES(; tol=1e1ctmrg_tol) - const svd_alg = SVDAdjoint(; fwd_alg, rrule_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 - ) - const gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme) -end +include("utility/defaults.jl") export SVDAdjoint, IterSVD, NonTruncSVDAdjoint -export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length +export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, VUMPS, VUMPSRuntime, VUMPSEnv, correlation_length export LocalOperator export expectation_value, costfun, product_peps export leading_boundary @@ -93,7 +57,6 @@ export PEPSOptimize, GeomSum, ManualIter, LinSolver export fixedpoint export InfinitePEPS, InfiniteTransferPEPS -export InfinitePEPO, InfiniteTransferPEPO export initializeMPS, initializePEPS export ReflectDepth, ReflectWidth, Rotate, RotateReflect export symmetrize!, symmetrize_retract_and_finalize! diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl new file mode 100644 index 00000000..ecab0335 --- /dev/null +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -0,0 +1,234 @@ + +""" + ρ = ρmap(ρ::Matrix{<:AbstractTensorMap}, A::Matrix{<:AbstractTensorMap}) +```` + ┌─ Aᵢⱼ─ ┌─ + ρᵢⱼ │ = ρⱼ₊₁ + └─ Aᵢⱼ─ └─ +```` +""" +function ρmap(ρ::Matrix{<:AbstractTensorMap}, A::Matrix{<:AbstractTensorMap}) + Ni, Nj = size(ρ) + ρ = deepcopy(ρ) + @inbounds for j in 1:Nj, i in 1:Ni + jr = mod1(j + 1, Nj) + @tensor ρ[i,jr][-1; -2] = ρ[i,j][4; 1] * A[i,j][1 2 3; -2] * conj(A[i,j][4 2 3; -1]) + end + return ρ +end + +""" + C = LRtoC(L::Matrix{<:AbstractTensorMap}, R::Matrix{<:AbstractTensorMap}) + +``` + ── Cᵢⱼ ── = ── Lᵢⱼ ── Rᵢⱼ₊₁ ── +``` +""" +function LRtoC(L::Matrix{<:AbstractTensorMap}, R::Matrix{<:AbstractTensorMap}) + Rijr = circshift(R, (0, -1)) + return L .* Rijr +end + +""" + AC = ALCtoAC(AL::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap}) + +``` + ── ACᵢⱼ ── = ── ALᵢⱼ ── Cᵢⱼ ── + | | +``` +""" +function ALCtoAC(AL::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap}) + return AL .* C +end + +""" + FLm = FLmap(FLi::Vector{<:AbstractTensorMap}, + ALui::Vector{<:AbstractTensorMap}, + ALdir::Vector{<:AbstractTensorMap}, + Ati::Vector{<:AbstractTensorMap}, + Abi::Vector{<:AbstractTensorMap}) + +``` + ┌── ┌── ALuᵢⱼ ── + │ │ │ +FLᵢⱼ₊₁ = FLᵢⱼ ─ Oᵢⱼ ── + │ │ │ + └── └── ALdᵢᵣⱼ ─ +``` +""" +function FLmap(FLi::Vector{<:AbstractTensorMap}, + ALui::Vector{<:AbstractTensorMap}, + ALdir::Vector{<:AbstractTensorMap}, + Ati::Vector{<:AbstractTensorMap}, + Abi::Vector{<:AbstractTensorMap}) + FLm = [@tensoropt FL[-1 -2 -3; -4] := FL[6 5 4; 1] * ALu[1 2 3; -4] * At[9; 2 -2 8 5] * + Ab[3 -3 7 4; 9] * ALd[-1; 6 8 7] for (FL, ALu, ALd, At, Ab) in zip(FLi, ALui, ALdir, Ati, Abi)] + + return circshift(FLm, 1) +end + +""" + ``` + ┌── ALuᵢⱼ ── ┌── + Lᵢⱼ | = Lᵢⱼ₊₁ + └── ALdᵢᵣⱼ ── └── + ``` +""" +function Lmap(Li::Vector{<:AbstractTensorMap}, + ALui::Vector{<:AbstractTensorMap}, + ALdir::Vector{<:AbstractTensorMap}) + Lm = [@tensoropt L[-6; -4] := ALu[1 2 3; -4] * L[5; 1] * ALd[-6; 5 2 3] for (L, ALu, ALd) in zip(Li, ALui, ALdir)] + + return circshift(Lm, 1) +end + +""" + FRm = FRmap(FRi::Vector{<:AbstractTensorMap}, + ARui::Vector{<:AbstractTensorMap}, + ARdir::Vector{<:AbstractTensorMap}, + Ati::Vector{<:AbstractTensorMap}, + Abi::Vector{<:AbstractTensorMap}) + +``` + ── ARuᵢⱼ ──┐ ──┐ + │ │ │ + ── Oᵢⱼ ──FRᵢⱼ = ──FRᵢⱼ₋₁ + │ │ │ + ── ARdᵢᵣⱼ ──┘ ──┘ +``` +""" +function FRmap(FRi::Vector{<:AbstractTensorMap}, + ARui::Vector{<:AbstractTensorMap}, + ARdir::Vector{<:AbstractTensorMap}, + Ati::Vector{<:AbstractTensorMap}, + Abi::Vector{<:AbstractTensorMap}) + FRm = [@tensoropt FR[-1 -2 -3; -4] := ARu[-1 1 2; 3] * FR[3 4 5; 8] * At[9; 1 4 7 -2] * + Ab[2 5 6 -3; 9] * ARd[8; -4 7 6] for (FR, ARu, ARd, At, Ab) in zip(FRi, ARui, ARdir, Ati, Abi)] + + return circshift(FRm, -1) +end + +""" + Rm = Rmap(FRi::Vector{<:AbstractTensorMap}, + ARui::Vector{<:AbstractTensorMap}, + ARdir::Vector{<:AbstractTensorMap}, + ) + +``` + ── ARuᵢⱼ ──┐ ──┐ + │ Rᵢⱼ = Rᵢⱼ₋₁ + ── ARdᵢᵣⱼ ──┘ ──┘ +``` +""" +function Rmap(Ri::Vector{<:AbstractTensorMap}, + ARui::Vector{<:AbstractTensorMap}, + ARdir::Vector{<:AbstractTensorMap}) + Rm = [@tensoropt R[-1; -5] := ARu[-1 2 3; 4] * R[4; 6] * ARd[6; -5 2 3] for (R, ARu, ARd) in zip(Ri, ARui, ARdir)] + + return circshift(Rm, -1) +end + +""" + ACm = ACmap(ACj::Vector{<:AbstractTensorMap}, + FLj::Vector{<:AbstractTensorMap}, + FRj::Vector{<:AbstractTensorMap}, + Atj::Vector{<:AbstractTensorMap}, + Abj::Vector{<:AbstractTensorMap}) + +``` + ┌─────── ACᵢⱼ ─────┐ +┌───── ACᵢ₊₁ⱼ ─────┐ │ │ │ +│ │ │ = FLᵢⱼ ─── Oᵢⱼ ───── FRᵢⱼ + │ │ │ + +``` +""" +function ACmap(ACj::Vector{<:AbstractTensorMap}, + FLj::Vector{<:AbstractTensorMap}, + FRj::Vector{<:AbstractTensorMap}, + Atj::Vector{<:AbstractTensorMap}, + Abj::Vector{<:AbstractTensorMap}) + ACm = [@tensoropt AC[-1 -2 -3; -4] := AC[1 2 3; 4] * FL[-1 6 5; 1]* At[9; 2 7 -2 6] * + Ab[3 8 -3 5; 9] * FR[4 7 8; -4] for (AC, FL, FR, At, Ab) in zip(ACj, FLj, FRj, Atj, Abj)] + + return circshift(ACm, 1) +end + +""" + Cmap(Cij, FLjp, FRj, II) + +``` + ┌────Cᵢⱼ ───┐ +┌── Cᵢ₊₁ⱼ ──┐ │ │ +│ │ = FLᵢⱼ₊₁ ──── FRᵢⱼ + │ │ + +``` +""" +function Cmap(Cj::Vector{<:AbstractTensorMap}, + FLjr::Vector{<:AbstractTensorMap}, + FRj::Vector{<:AbstractTensorMap}) + Cm = [@tensoropt C[-1; -2] := C[1; 2] * FL[-1 3 4; 1] * FR[2 3 4; -2] for (C, FL, FR) in zip(Cj, FLjr, FRj)] + + return circshift(Cm, 1) +end + +""" + nearest_neighbour_energy(ipeps::InfinitePEPS, Hh, Hv, env::VUMPSEnv) + + Compute the energy of the nearest neighbour Hamiltonian for an infinite PEPS. + +``` + ┌────── ACuᵢⱼ ────ARuᵢᵣ ──────┐ + │ │ │ │ + FLoᵢⱼ ── Oᵢⱼ ────── Oᵢᵣ ──── FRoᵢᵣ ir = Ni + 1 - i + │ │ │ │ jr = j + 1 + └───── ACdᵢᵣⱼ ─────ARdᵢᵣᵣ─────┘ + + ┌─────── ACuᵢⱼ ─────┐ + │ │ │ + FLuᵢⱼ ─── Oᵢⱼ ───── FRuᵢⱼ ir = i + 1 + │ │ │ irr = Ni - i + FLoᵢᵣⱼ ── Oᵢᵣⱼ ─── FRoᵢᵣⱼ + │ │ │ + └────── ACdᵢᵣᵣⱼ ───┘ +``` +""" +function nearest_neighbour_energy(ipeps::InfinitePEPS, Hh, Hv, env::VUMPSEnv) + @unpack ACu, ARu, ACd, ARd, FLu, FRu, FLo, FRo = env + Ni, Nj = size(ipeps) + + energy_tol = 0 + for j in 1:Nj, i in 1:Ni + # horizontal contraction + ir = Ni + 1 - i + jr = mod1(j + 1, Nj) + @tensoropt oph[-1 -2; -3 -4] := FLo[i,j][18 12 15; 5] * ACu[i,j][5 6 7; 8] * ipeps.A[i,j][-1; 6 13 19 12] * + conj(ipeps.A[i,j][-3; 7 16 20 15]) * conj(ACd[ir,j][18 19 20; 21]) * + ARu[i,jr][8 9 10; 11] * ipeps.A[i,jr][-2; 9 14 22 13] * + conj(ipeps.A[i,jr][-4; 10 17 23 16]) * conj(ARd[ir,jr][21 22 23; 24]) * FRo[i,jr][11 14 17; 24] + + @tensor eh = oph[1 2; 3 4] * Hh[3 4; 1 2] + @tensor nh = oph[1 2; 1 2] + energy_tol += eh / nh + # @show eh / nh eh nh + + # vertical contraction + ir = mod1(i + 1, Ni) + irr = mod1(Ni - i, Ni) + @tensoropt opv[-1 -2; -3 -4] := FLu[i,j][21 19 20; 18] * ACu[i,j][18 12 15 5] * ipeps.A[i,j][-1; 12 6 13 19] * + conj(ipeps.A[i,j][-3; 15 7 16 20]) * FRu[i,j][5 6 7; 8] * FLo[ir,j][24 22 23; 21] * + ipeps.A[ir,j][-2; 13 9 14 22] * conj(ipeps.A[ir,j][-4; 16 10 17 23]) * + FRo[ir,j][8 9 10; 11] * conj(ACd[irr,j][24 14 17; 11]) + + @tensor ev = opv[1 2; 3 4] * Hv[3 4; 1 2] + @tensor nv = opv[1 2; 1 2] + energy_tol += ev / nv + # @show ev / nv ev nv + + # penalty term + # energy_tol += 0.1 * abs(eh / nh - eh / nh) + end + + return energy_tol/Ni/Nj +end \ No newline at end of file diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl index bba37f25..2cfacdf1 100644 --- a/src/algorithms/ctmrg/ctmrg.jl +++ b/src/algorithms/ctmrg/ctmrg.jl @@ -43,8 +43,8 @@ function svd_algorithm(alg::ProjectorAlg, (dir, r, c)) end """ - CTMRG(; tol=Defaults.ctmrg_tol, maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, verbosity=0, + CTMRG(; tol=Defaults.contr_tol, maxiter=Defaults.contr_maxiter, + miniter=Defaults.contr_miniter, verbosity=0, svd_alg=SVDAdjoint(), trscheme=FixedSpaceTruncation(), ctmrgscheme=Defaults.ctmrgscheme) @@ -72,9 +72,9 @@ struct CTMRG{S} projector_alg::ProjectorAlg end function CTMRG(; - tol=Defaults.ctmrg_tol, - maxiter=Defaults.ctmrg_maxiter, - miniter=Defaults.ctmrg_miniter, + tol=Defaults.contr_tol, + maxiter=Defaults.contr_maxiter, + miniter=Defaults.contr_miniter, verbosity=2, svd_alg=Defaults.svd_alg, trscheme=Defaults.trscheme, @@ -92,15 +92,15 @@ const SequentialCTMRG = CTMRG{:sequential} const SimultaneousCTMRG = CTMRG{:simultaneous} """ - MPSKit.leading_boundary([envinit], state, alg::CTMRG) + leading_boundary([envinit], state, alg::CTMRG) Contract `state` using CTMRG and return the CTM environment. Per default, a random initial environment is used. """ -function MPSKit.leading_boundary(state, alg::CTMRG) - return MPSKit.leading_boundary(CTMRGEnv(state, oneunit(spacetype(state))), state, alg) +function leading_boundary(state, alg::CTMRG) + return leading_boundary(CTMRGEnv(state, oneunit(spacetype(state))), state, alg) end -function MPSKit.leading_boundary(envinit, state, alg::CTMRG) +function leading_boundary(envinit, state, alg::CTMRG) CS = map(x -> tsvd(x; alg=TensorKit.SVD())[2], envinit.corners) TS = map(x -> tsvd(x; alg=TensorKit.SVD())[2], envinit.edges) @@ -110,22 +110,22 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRG) log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG")) return LoggingExtras.withlevel(; alg.verbosity) do - ctmrg_loginit!(log, η, N) + contr_loginit!(log, η, N) for iter in 1:(alg.maxiter) env, = ctmrg_iter(state, env, alg) # Grow and renormalize in all 4 directions η, CS, TS = calc_convergence(env, CS, TS) N = norm(state, env) - ctmrg_logiter!(log, iter, η, N) + contr_logiter!(log, iter, η, N) if η ≤ alg.tol - ctmrg_logfinish!(log, iter, η, N) + contr_logfinish!(log, iter, η, N) break end if iter == alg.maxiter - ctmrg_logcancel!(log, iter, η, N) + contr_logcancel!(log, iter, η, N) else - ctmrg_logiter!(log, iter, η, N) + contr_logiter!(log, iter, η, N) end end return env @@ -161,16 +161,6 @@ function ctmrg_iter(state, envs::CTMRGEnv, alg::SimultaneousCTMRG) return envs′, info end -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) -ctmrg_logcancel!(log, iter, η, N) = @warnv 1 logcancel!(log, iter, η, N) - -@non_differentiable ctmrg_loginit!(args...) -@non_differentiable ctmrg_logiter!(args...) -@non_differentiable ctmrg_logfinish!(args...) -@non_differentiable ctmrg_logcancel!(args...) - # ======================================================================================== # # Expansion step # ======================================================================================== # diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl index 3b569e6b..09fcf3f3 100644 --- a/src/algorithms/peps_opt.jl +++ b/src/algorithms/peps_opt.jl @@ -89,13 +89,13 @@ step by setting `reuse_env` to true. Otherwise a random environment is used at e step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm. """ struct PEPSOptimize{G} - boundary_alg::CTMRG + boundary_alg::Union{CTMRG, VUMPS} optimizer::OptimKit.OptimizationAlgorithm reuse_env::Bool gradient_alg::G function PEPSOptimize( # Inner constructor to prohibit illegal setting combinations - boundary_alg::CTMRG{S}, + boundary_alg::Union{CTMRG{S}, VUMPS}, optimizer, reuse_env, gradient_alg::G, @@ -106,6 +106,8 @@ struct PEPSOptimize{G} elseif boundary_alg.projector_alg.svd_alg.fwd_alg isa IterSVD && iterscheme(gradient_alg) === :fixed throw(ArgumentError("IterSVD and :fixed are currently not compatible")) + elseif boundary_alg isa VUMPS && iterscheme(gradient_alg) === :fixed + throw(ArgumentError("VUMPS and :fixed are currently not compatible")) end end return new{G}(boundary_alg, optimizer, reuse_env, gradient_alg) @@ -148,7 +150,7 @@ function fixedpoint( ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, - env₀::CTMRGEnv=CTMRGEnv(ψ₀, field(T)^20); + env₀=CTMRGEnv(ψ₀, field(T)^20); (finalize!)=OptimKit._finalize!, symmetrization=nothing, ) where {T} @@ -177,6 +179,7 @@ function fixedpoint( return costfun(ψ, envs´, H) end g = only(gs) # `withgradient` returns tuple of gradients `gs` + envs isa Union{VUMPSRuntime, Tuple{VUMPSRuntime, VUMPSRuntime}} && (g = InfinitePEPS(g.A)) # KrylovKit patch return E, g end @@ -213,11 +216,12 @@ Evaluating the gradient of the cost function for CTMRG: function _rrule( gradmode::GradMode{:diffgauge}, ::RuleConfig, - ::typeof(MPSKit.leading_boundary), + ::typeof(leading_boundary), envinit, state, alg::CTMRG, ) + envs = leading_boundary(envinit, state, alg) #TODO: fixed space for unit cells function leading_boundary_diffgauge_pullback(Δenvs′) @@ -240,15 +244,47 @@ function _rrule( return envs, leading_boundary_diffgauge_pullback end +function _rrule( + gradmode::GradMode{:diffgauge}, + ::RuleConfig, + ::typeof(leading_boundary), + envinit, + state, + alg::VUMPS, +) + + 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 + + # evaluate the geometric sum + ∂f∂A(x)::typeof(state) = env_vjp(x)[1] + ∂f∂x(x)::typeof(envs) = env_vjp(x)[2] + ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode) + + return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent() + end + + return envs, leading_boundary_diffgauge_pullback +end + # Here f is differentiated from an pre-computed SVD with fixed U, S and V function _rrule( gradmode::GradMode{:fixed}, ::RuleConfig, - ::typeof(MPSKit.leading_boundary), + ::typeof(leading_boundary), envinit, state, alg::CTMRG{C}, ) where {C} + @assert C === :simultaneous @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg) envs = leading_boundary(envinit, state, alg) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 4f48cdf4..787db82a 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -1,4 +1,4 @@ -function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CTMRGEnv) +function expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CTMRGEnv) checklattice(peps, O) return sum(O.terms) do (inds, operator) # TODO: parallelize this map, especially for the backwards pass contract_localoperator(inds, operator, peps, peps, envs) / @@ -6,8 +6,16 @@ function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CT end end -function costfun(peps::InfinitePEPS, envs::CTMRGEnv, O::LocalOperator) - E = MPSKit.expectation_value(peps, O, envs) +function expectation_value(peps::InfinitePEPS, O::LocalOperator, rt::Union{VUMPSRuntime, Tuple{VUMPSRuntime, VUMPSRuntime}}) + # checklattice(peps, O) + env = VUMPSEnv(rt, peps) + Hh = O.terms[1].second + Hv = O.terms[2].second + return nearest_neighbour_energy(peps, Hh, Hv, env) +end + +function costfun(peps::InfinitePEPS, envs, O::LocalOperator) + E = expectation_value(peps, O, envs) ignore_derivatives() do isapprox(imag(E), 0; atol=sqrt(eps(real(E)))) || @warn "Expectation value is not real: $E." @@ -56,6 +64,37 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) return total end +function LinearAlgebra.norm(ipeps::InfinitePEPS, env::VUMPSEnv) + @unpack ACu, ARu, ACd, ARd, FLu, FRu, FLo, FRo = env + Ni, Nj = size(ipeps) + + λFLo, _ = rightenv(ARu, adjoint.(ARd), ipeps; ifobs=true) + λC, _ = rightCenv(ARu, adjoint.(ARd); ifobs=true) + + return prod(λFLo./λC)^(1/Ni) + + # test for 1x1 unitcell + # _, FLo = leftenv(ARu, adjoint.(ARd), itp; ifobs=true) + # _, FRo = rightenv(ARu, adjoint.(ARd), itp; ifobs=true) + # _, Lo = leftCenv(ARu, adjoint.(ARd); ifobs=true) + # _, Ro = rightCenv(ARu, adjoint.(ARd); ifobs=true) + + # λFLo = Zygote.Buffer(zeros(ComplexF64, Ni, Nj)) + # λCo = Zygote.Buffer(zeros(ComplexF64, Ni, Nj)) + # for i in 1:Ni, j in 1:Nj + # ir = Ni + 1 - i + # irr = mod1(Ni + 2 - i, Ni) + # @tensoropt FLm[-1 -2 -3; -4] := FLo[i,j][6 5 4; 1] * ARu[i,j][1 2 3; -4] * itp.top[i,j][9; 2 -2 8 5] * + # itp.bot[i,j][3 -3 7 4; 9] * adjoint(ARd[ir,j])[-1; 6 8 7] + # @tensoropt Lm[-6; -4] := ARu[i,j][1 2 3; -4] * Lo[i,j][5; 1] * adjoint(ARd[irr,j])[-6; 5 2 3] + # λFLo[i, j] = (@tensor FLm[1 2 3; 4] * FRo[i, j][4 2 3; 1]) / + # (@tensor FLo[i, j][1 2 3; 4] * FRo[i, j][4 2 3; 1]) + # λCo[i, j] = (@tensor Lm[1; 2] * Ro[i, j][2; 1]) / + # (@tensor Lo[i, j][1; 2] * Ro[i, j][2; 1]) + # end + # return prod(copy(λFLo)./copy(λCo))^(1/Ni/Nj) +end + """ correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2) diff --git a/src/algorithms/vumps/vumps.jl b/src/algorithms/vumps/vumps.jl new file mode 100644 index 00000000..59c5589e --- /dev/null +++ b/src/algorithms/vumps/vumps.jl @@ -0,0 +1,120 @@ +@kwdef mutable struct VUMPS + ifupdown::Bool = true + tol::Float64 = Defaults.contr_tol + maxiter::Int = Defaults.contr_maxiter + miniter::Int = Defaults.contr_miniter + verbosity::Int = Defaults.verbosity +end + +function VUMPSRuntime(ipeps::InfinitePEPS, χ::VectorSpace) + A = initial_A(ipeps, χ) + AL, L, _ = left_canonical(A) + R, AR, _ = right_canonical(AL) + + _, FL = leftenv(AL, adjoint.(AL), ipeps) + _, FR = rightenv(AR, adjoint.(AR), ipeps) + C = LRtoC(L, R) + return VUMPSRuntime(AL, AR, C, FL, FR) +end +@non_differentiable VUMPSRuntime(ipeps::InfinitePEPS, χ::VectorSpace) + +function _down_ipeps(ipeps::InfinitePEPS) + Ni, Nj = size(ipeps) + ipepsd = Zygote.Buffer(ipeps.A) + for j in 1:Nj, i in 1:Ni + ir = Ni + 1 - i + Ad = permute(ipeps[ir,j]', ((5,), (3,2,1,4))) + ipepsd[i, j] = _fit_spaces(Ad, ipeps[ir,j]) + end + + return InfinitePEPS(copy(ipepsd)) +end + +VUMPSRuntime(ipeps::InfinitePEPS, χ::Int, alg::VUMPS) = VUMPSRuntime(ipeps, ℂ^χ, alg) +function VUMPSRuntime(ipeps::InfinitePEPS, χ::VectorSpace, alg::VUMPS) + Ni, Nj = size(ipeps) + + rtup = VUMPSRuntime(ipeps, χ) + alg.verbosity >= 2 && Zygote.@ignore @info "VUMPS init: cell=($(Ni)×$(Nj)) χ = $(χ) up(↑) environment" + + if alg.ifupdown + ipepsd = _down_ipeps(ipeps) + rtdown = VUMPSRuntime(ipepsd, χ) + alg.verbosity >= 2 && Zygote.@ignore @info "VUMPS init: cell=($(Ni)×$(Nj)) χ = $(χ) down(↓) environment" + + return rtup, rtdown + else + return rtup + end +end +@non_differentiable VUMPSRuntime(ipeps::InfinitePEPS, χ::VectorSpace, alg::VUMPS) + +function vumps_itr(rt::VUMPSRuntime, ipeps::InfinitePEPS, alg::VUMPS) + t = Zygote.@ignore time() + for i in 1:alg.maxiter + rt, err = vumps_step(rt, ipeps, alg) + alg.verbosity >= 3 && Zygote.@ignore @info @sprintf("VUMPS@step: %4d\terr = %.3e\ttime = %.3f sec", i, err, time()-t) + if err < alg.tol && i >= alg.miniter + alg.verbosity >= 2 && Zygote.@ignore @info @sprintf("VUMPS conv@step: %4d\terr = %.3e\ttime = %.3f sec", i, err, time()-t) + break + end + if i == alg.maxiter + alg.verbosity >= 2 && Zygote.@ignore @warn @sprintf("VUMPS cancel@step: %4d\terr = %.3e\ttime = %.3f sec", i, err, time()-t) + end + end + return rt +end + +function leading_boundary(rt::VUMPSRuntime, ipeps::InfinitePEPS, alg::VUMPS) + return vumps_itr(rt, ipeps, alg) +end + +function VUMPSEnv(rt::VUMPSRuntime, ::InfinitePEPS) + @unpack AL, AR, C, FL, FR = rt + AC = ALCtoAC(AL, C) + return VUMPSEnv(AC, AR, AC, AR, FL, FR, FL, FR) +end + +function leading_boundary(rt::Tuple{VUMPSRuntime, VUMPSRuntime}, ipeps::InfinitePEPS, alg::VUMPS) + rtup, rtdown = rt + + rtup = vumps_itr(rtup, ipeps, alg) + + ipepsd = _down_ipeps(ipeps) + rtdown = vumps_itr(rtdown, ipepsd, alg) + + return rtup, rtdown +end + +function VUMPSEnv(rt::Tuple{VUMPSRuntime, VUMPSRuntime}, ipeps::InfinitePEPS) + rtup, rtdown = rt + + ALu, ARu, Cu, FLu, FRu = rtup.AL, rtup.AR, rtup.C, rtup.FL, rtup.FR + ACu = ALCtoAC(ALu, Cu) + + ALd, ARd, Cd = rtdown.AL, rtdown.AR, rtdown.C + ACd = ALCtoAC(ALd, Cd) + + _, FLo = leftenv(ALu, adjoint.(ALd), ipeps, FLu; ifobs = true) + _, FRo = rightenv(ARu, adjoint.(ARd), ipeps, FRu; ifobs = true) + + return VUMPSEnv(ACu, ARu, ACd, ARd, FLu, FRu, FLo, FRo) +end + +function vumps_step(rt::VUMPSRuntime, ipeps::InfinitePEPS, alg::VUMPS) + verbosity = alg.verbosity + @unpack AL, C, AR, FL, FR = rt + AC = Zygote.@ignore ALCtoAC(AL,C) + _, ACp = ACenv(AC, FL, FR, ipeps; verbosity) + _, Cp = Cenv( C, FL, FR; verbosity) + ALp, ARp, _, _ = ACCtoALAR(ACp, Cp) + _, FL = leftenv(AL, adjoint.(ALp), ipeps, FL; verbosity) + _, FR = rightenv(AR, adjoint.(ARp), ipeps, FR; verbosity) + _, ACp = ACenv(ACp, FL, FR, ipeps; verbosity) + _, Cp = Cenv( Cp, FL, FR; verbosity) + ALp, ARp, errL, errR = ACCtoALAR(ACp, Cp) + err = errL + errR + alg.verbosity >= 4 && err > 1e-8 && println("errL=$errL, errR=$errR") + + return VUMPSRuntime(ALp, ARp, Cp, FL, FR), err +end \ No newline at end of file diff --git a/src/environments/transferpepo_environments.jl b/src/environments/transferpepo_environments.jl deleted file mode 100644 index f79d641b..00000000 --- a/src/environments/transferpepo_environments.jl +++ /dev/null @@ -1,103 +0,0 @@ - -function MPSKit.environments(state::InfiniteMPS, O::InfiniteTransferPEPO; kwargs...) - return environments( - convert(MPSMultiline, state), convert(TransferPEPOMultiline, O); kwargs... - ) -end - -function MPSKit.environments( - state::MPSMultiline, O::TransferPEPOMultiline; solver=MPSKit.Defaults.eigsolver -) - (lw, rw) = MPSKit.mixed_fixpoints(state, O, state; solver) - return MPSKit.PerMPOInfEnv(nothing, O, state, solver, lw, rw, ReentrantLock()) -end - -function MPSKit.mixed_fixpoints( - above::MPSMultiline, - O::TransferPEPOMultiline, - below::MPSMultiline, - init=gen_init_fps(above, O, below); - solver=MPSKit.Defaults.eigsolver, -) - T = eltype(above) - - (numrows, numcols) = size(above) - @assert size(above) == size(O) - @assert size(below) == size(O) - - envtype = eltype(init[1]) - lefties = PeriodicArray{envtype,2}(undef, numrows, numcols) - righties = PeriodicArray{envtype,2}(undef, numrows, numcols) - - @threads for cr in 1:numrows - c_above = above[cr] # TODO: Update index convention to above[cr - 1] - c_below = below[cr + 1] - - (L0, R0) = init[cr] - - @sync begin - Threads.@spawn begin - E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL) - (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver) - convhist.converged < 1 && - @info "left eigenvalue failed to converge $(convhist.normres)" - L0 = first(Ls) - end - - Threads.@spawn begin - E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR) - (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver) - convhist.converged < 1 && - @info "right eigenvalue failed to converge $(convhist.normres)" - R0 = first(Rs) - end - end - - lefties[cr, 1] = L0 - for loc in 2:numcols - lefties[cr, loc] = - lefties[cr, loc - 1] * - TransferMatrix(c_above.AL[loc - 1], O[cr, loc - 1], c_below.AL[loc - 1]) - end - - renormfact::scalartype(T) = dot(c_below.CR[0], PEPO_∂∂C(L0, R0) * c_above.CR[0]) - - righties[cr, end] = R0 / sqrt(renormfact) - lefties[cr, 1] /= sqrt(renormfact) - - for loc in (numcols - 1):-1:1 - righties[cr, loc] = - TransferMatrix(c_above.AR[loc + 1], O[cr, loc + 1], c_below.AR[loc + 1]) * - righties[cr, loc + 1] - - renormfact = dot( - c_below.CR[loc], - PEPO_∂∂C(lefties[cr, loc + 1], righties[cr, loc]) * c_above.CR[loc], - ) - righties[cr, loc] /= sqrt(renormfact) - lefties[cr, loc + 1] /= sqrt(renormfact) - end - end - - return (lefties, righties) -end - -function gen_init_fps(above::MPSMultiline, O::TransferPEPOMultiline, below::MPSMultiline) - T = eltype(above) - - map(1:size(O, 1)) do cr - L0::T = TensorMap( - rand, - scalartype(T), - left_virtualspace(below, cr + 1, 0) * prod(adjoint.(west_spaces(O[cr], 1))), - left_virtualspace(above, cr, 0), # TODO: Update index convention to above[cr - 1] - ) - R0::T = TensorMap( - rand, - scalartype(T), - right_virtualspace(above, cr, 0) * prod(adjoint.(east_spaces(O[cr], 1))), - right_virtualspace(below, cr + 1, 0), - ) - (L0, R0) - end -end diff --git a/src/environments/transferpeps_environments.jl b/src/environments/transferpeps_environments.jl deleted file mode 100644 index 4332c232..00000000 --- a/src/environments/transferpeps_environments.jl +++ /dev/null @@ -1,137 +0,0 @@ - -function MPSKit.environments(state::InfiniteMPS, O::InfiniteTransferPEPS; kwargs...) - return environments( - convert(MPSMultiline, state), convert(TransferPEPSMultiline, O); kwargs... - ) -end - -import MPSKit.MPSMultiline - -function MPSKit.environments( - state::MPSMultiline, O::TransferPEPSMultiline; solver=MPSKit.Defaults.eigsolver -) - (lw, rw) = MPSKit.mixed_fixpoints(state, O, state; solver) - return MPSKit.PerMPOInfEnv(nothing, O, state, solver, lw, rw, ReentrantLock()) -end - -function MPSKit.mixed_fixpoints( - above::MPSMultiline, - O::TransferPEPSMultiline, - below::MPSMultiline, - init=gen_init_fps(above, O, below); - solver=MPSKit.Defaults.eigsolver, -) - T = eltype(above) - - (numrows, numcols) = size(above) - @assert size(above) == size(O) - @assert size(below) == size(O) - - envtype = eltype(init[1]) - lefties = PeriodicArray{envtype,2}(undef, numrows, numcols) - righties = PeriodicArray{envtype,2}(undef, numrows, numcols) - - @threads for cr in 1:numrows - c_above = above[cr] # TODO: Update index convention to above[cr - 1] - c_below = below[cr + 1] - - (L0, R0) = init[cr] - - @sync begin - Threads.@spawn begin - E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL) - (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver) - convhist.converged < 1 && - @info "left eigenvalue failed to converge $(convhist.normres)" - L0 = first(Ls) - end - - Threads.@spawn begin - E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR) - (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver) - convhist.converged < 1 && - @info "right eigenvalue failed to converge $(convhist.normres)" - R0 = first(Rs) - end - end - - lefties[cr, 1] = L0 - for loc in 2:numcols - lefties[cr, loc] = - lefties[cr, loc - 1] * - TransferMatrix(c_above.AL[loc - 1], O[cr, loc - 1], c_below.AL[loc - 1]) - end - - renormfact::scalartype(T) = dot(c_below.CR[0], PEPS_∂∂C(L0, R0) * c_above.CR[0]) - - righties[cr, end] = R0 / sqrt(renormfact) - lefties[cr, 1] /= sqrt(renormfact) - - for loc in (numcols - 1):-1:1 - righties[cr, loc] = - TransferMatrix(c_above.AR[loc + 1], O[cr, loc + 1], c_below.AR[loc + 1]) * - righties[cr, loc + 1] - - renormfact = dot( - c_below.CR[loc], - PEPS_∂∂C(lefties[cr, loc + 1], righties[cr, loc]) * c_above.CR[loc], - ) - righties[cr, loc] /= sqrt(renormfact) - lefties[cr, loc + 1] /= sqrt(renormfact) - end - end - - return (lefties, righties) -end - -function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSMultiline) - T = eltype(above) - - map(1:size(O, 1)) do cr - L0::T = TensorMap( - rand, - scalartype(T), - left_virtualspace(below, cr + 1, 0) * - space(O[cr].top[1], 5)' * - space(O[cr].bot[1], 5), - left_virtualspace(above, cr, 0), # TODO: Update index convention to above[cr - 1] - ) - R0::T = TensorMap( - rand, - scalartype(T), - right_virtualspace(above, cr, 0) * - space(O[cr].top[1], 3)' * - space(O[cr].bot[1], 3), - right_virtualspace(below, cr + 1, 0), - ) - (L0, R0) - end -end - -function MPSKit.transfer_spectrum( - above::MPSMultiline, - O::TransferPEPSMultiline, - below::MPSMultiline, - init=gen_init_fps(above, O, below); - num_vals=2, - solver=MPSKit.Defaults.eigsolver, -) - @assert size(above) == size(O) - @assert size(below) == size(O) - - numrows = size(above, 1) - envtype = eltype(init[1]) - eigenvals = Vector{Vector{scalartype(envtype)}}(undef, numrows) - - @threads for cr in 1:numrows - L0, = init[cr] - - E_LL = TransferMatrix(above[cr - 1].AL, O[cr], below[cr + 1].AL) # Note that this index convention is different from above! - λ, _, convhist = eigsolve(flip(E_LL), L0, num_vals, :LM, solver) - convhist.converged < num_vals && - @warn "correlation length failed to converge: normres = $(convhist.normres)" - eigenvals[cr] = λ - end - - return eigenvals -end diff --git a/src/environments/vumps_environments.jl b/src/environments/vumps_environments.jl new file mode 100644 index 00000000..33d28650 --- /dev/null +++ b/src/environments/vumps_environments.jl @@ -0,0 +1,516 @@ +""" + VUMPSEnv{T<:Number, S<:IndexSpace, + OT<:AbstractTensorMap{S, 2, 2}, + ET<:AbstractTensorMap{S, 2, 1}, + CT<:AbstractTensorMap{S, 1, 1}} + +A struct that contains the environment of the VUMPS algorithm for calculate observables. + +For a `Ni` x `Nj` unitcell, each is a Matrix, containing + +- `AC`: The mixed canonical environment tensor. +- `AR`: The right canonical environment tensor. +- `Lu`: The left upper environment tensor. +- `Ru`: The right upper environment tensor. +- `Lo`: The left mixed environment tensor. +- `Ro`: The right mixed environment tensor. +""" +struct VUMPSEnv{T<:Number, S<:IndexSpace, + ET<:AbstractTensorMap{S, 3, 1}} + ACu::Matrix{ET} + ARu::Matrix{ET} + ACd::Matrix{ET} + ARd::Matrix{ET} + FLu::Matrix{ET} + FRu::Matrix{ET} + FLo::Matrix{ET} + FRo::Matrix{ET} + function VUMPSEnv(ACu::Matrix{ET}, + ARu::Matrix{ET}, + ACd::Matrix{ET}, + ARd::Matrix{ET}, + FLu::Matrix{ET}, + FRu::Matrix{ET}, + FLo::Matrix{ET}, + FRo::Matrix{ET}) where {ET} + T = eltype(ACu[1]) + S = spacetype(ACu[1]) + new{T, S, ET}(ACu, ARu, ACd, ARd, FLu, FRu, FLo, FRo) + end +end + +""" + VUMPSRuntime{T<:Number, S<:IndexSpace, + OT<:AbstractTensorMap{S, 2, 2}, + ET<:AbstractTensorMap{S, 2, 1}, + CT<:AbstractTensorMap{S, 1, 1}} + +A struct that contains the environment of the VUMPS algorithm for runtime calculations. + +For a `Ni` x `Nj` unitcell, each is a Matrix, containing + +- `O`: The center transfer matrix PEPO tensor. +- `AL`: The left canonical environment tensor. +- `AR`: The right canonical environment tensor. +- `C`: The canonical environment tensor. +- `L`: The left environment tensor. +- `R`: The right environment tensor. +""" +struct VUMPSRuntime{T<:Number, S<:IndexSpace, + ET<:AbstractTensorMap{S, 3, 1}, + CT<:AbstractTensorMap{S, 1, 1}} + AL::Matrix{ET} + AR::Matrix{ET} + C::Matrix{CT} + FL::Matrix{ET} + FR::Matrix{ET} + function VUMPSRuntime(AL::Matrix{ET}, + AR::Matrix{ET}, + C::Matrix{CT}, + FL::Matrix{ET}, + FR::Matrix{ET}) where {ET, CT} + T = eltype(AL[1]) + S = spacetype(AL[1]) + new{T, S, ET, CT}(AL, AR, C, FL, FR) + end +end + +# In-place update of environment +function update!(env::VUMPSRuntime, env´::VUMPSRuntime) + env.AL .= env´.AL + env.AR .= env´.AR + env.C .= env´.C + env.FL .= env´.FL + env.FR .= env´.FR + return env +end + +function update!(env::Tuple{VUMPSRuntime, VUMPSRuntime}, env´::Tuple{VUMPSRuntime, VUMPSRuntime}) + update!(env[1], env´[1]) + update!(env[2], env´[2]) + return env +end + +""" + +```` + + l ←------- r + / \ + / \ + t d +```` +Initalize a boundary MPS for the transfer operator `O` by specifying an array of virtual +spaces consistent with the unit cell. +""" +function initial_A(ipeps::InfinitePEPS, χ::VectorSpace) + T = eltype(ipeps[1]) + Ni, Nj = size(ipeps) + A = [(D = space(ipeps[i, j], 2)'; + TensorMap(rand, T, χ * D * D', χ)) for i in 1:Ni, j in 1:Nj] + return A +end + +function initial_C(A::Matrix{<:AbstractTensorMap}) + T = eltype(A[1]) + Ni, Nj = size(A) + C = [(χ = space(A[i, j], 1); + isomorphism(Matrix{T}, χ, χ)) for i in 1:Ni, j in 1:Nj] + return C +end + +# KrylovKit patch +TensorKit.inner(x::AbstractArray{<:AbstractTensorMap}, y::AbstractArray{<:AbstractTensorMap}) = sum(map(TensorKit.inner, x, y)) +TensorKit.add!!(x::AbstractArray{<:AbstractTensorMap}, y::AbstractArray{<:AbstractTensorMap}, a::Number, b::Number) = map((x, y) -> TensorKit.add!!(x, y, a, b), x, y) +TensorKit.scale!!(x::AbstractArray{<:AbstractTensorMap}, a::Number) = map(x -> TensorKit.scale!!(x, a), x) + +""" + λs[1], Fs[1] = selectpos(λs, Fs) + +Select the max positive one of λs and corresponding Fs. +""" +function selectpos(λs, Fs, N) + if length(λs) > 1 && norm(abs(λs[1]) - abs(λs[2])) < 1e-12 + # @show "selectpos: λs are degeneracy" + N = min(N, length(λs)) + p = argmax(real(λs[1:N])) + # @show λs p abs.(λs) + return λs[1:N][p], Fs[1:N][p] + else + return λs[1], Fs[1] + end +end + +""" + L = getL!(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap}; verbosity = Defaults.verbosity, kwargs...) + +```` + ┌─ Aᵢⱼ ─ Aᵢⱼ₊₁─ ┌─ L ─ + ρᵢⱼ │ │ = ρᵢⱼ = │ + └─ Aᵢⱼ─ Aᵢⱼ₊₁─ └─ L'─ +```` + +ρ=L'*L, return L, where `L`is guaranteed to have positive diagonal elements. + +""" +function getL!(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap}; verbosity = Defaults.verbosity, kwargs...) + Ni, Nj = size(A) + λs, ρs, info = eigsolve(ρ -> ρmap(ρ, A), L, 1, :LM; ishermitian = false, maxiter = 1, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "getL not converged" + _, ρs1 = selectpos(λs, ρs, Nj) + @inbounds for j = 1:Nj, i = 1:Ni + ρ = ρs1[i,j] + ρs1[i,j]' + ρ /= tr(ρ) + _, S, Vt = tsvd!(ρ) + # Lo = lmul!(Diagonal(sqrt.(F.S)), F.Vt) + Lo = sqrt(S) * Vt + _, R = leftorth!(Lo) + L[i,j] = R + end + return L +end + +function _to_front(t::AbstractTensorMap) # make TensorMap{S,N₁+N₂-1,1} + I1 = TensorKit.codomainind(t) + I2 = TensorKit.domainind(t) + return transpose(t, ((I1..., reverse(Base.tail(I2))...), (I2[1],))) +end + +function _to_tail(t::AbstractTensorMap) # make TensorMap{S,1,N₁+N₂-1} + I1 = TensorKit.codomainind(t) + I2 = TensorKit.domainind(t) + return transpose(t, ((I1[1],), (I2..., reverse(Base.tail(I1))...))) +end + +""" + AL, Le, λ = getAL(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap}) + +Given an MPS tensor `A` and `L` ，return a left-canonical MPS tensor `AL`, a gauge transform `R` and +a scalar factor `λ` such that ``λ * AL * Le = L * A`` +""" +function getAL(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap}) + Ni, Nj = size(A) + AL = similar(A) + Le = similar(L) + λ = zeros(Ni, Nj) + @inbounds for j = 1:Nj, i = 1:Ni + Q, R = leftorth!(_to_front(L[i,j] * _to_tail(A[i,j]))) + AL[i,j] = Q + λ[i,j] = norm(R) + Le[i,j] = rmul!(R, 1/λ[i,j]) + end + return AL, Le, λ +end + +function getLsped(Le::Matrix{<:AbstractTensorMap}, A::Matrix{<:AbstractTensorMap}, AL::Matrix{<:AbstractTensorMap}; verbosity = Defaults.verbosity, kwargs...) + Ni, Nj = size(A) + L = similar(Le) + @inbounds for j = 1:Nj, i = 1:Ni + λs, Ls, info = eigsolve(L -> (@tensor Ln[-1; -2] := L[4; 1] * A[i,j][1 2 3; -2] * conj(AL[i,j][4 2 3; -1])), Le[i,j], 1, :LM; ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "getLsped not converged" + _, Ls1 = selectpos(λs, Ls, Nj) + _, R = leftorth!(Ls1) + L[i,j] = R + end + return L +end + +""" + AL, L, λ = left_canonical(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap} = initial_C(A); kwargs...) + +Given an MPS tensor `A`, return a left-canonical MPS tensor `AL`, a gauge transform `L` and +a scalar factor `λ` such that ``λ*AL*L = L*A``, where an initial guess for `L` can be +provided. +""" +function left_canonical(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap} = initial_C(A); tol = 1e-12, maxiter = 100, kwargs...) + L = getL!(A, L; kwargs...) + AL, Le, λ = getAL(A, L;kwargs...) + numiter = 1 + while norm(L.-Le) > tol && numiter < maxiter + L = getLsped(Le, A, AL; kwargs...) + AL, Le, λ = getAL(A, L; kwargs...) + numiter += 1 + end + L = Le + return AL, L, λ +end + +""" + R, AR, λ = right_canonical(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap} = initial_C(A); tol = 1e-12, maxiter = 100, kwargs...) + +Given an MPS tensor `A`, return a gauge transform R, a right-canonical MPS tensor `AR`, and +a scalar factor `λ` such that ``λ * R * AR = A * R``, where an initial guess for `R` can be +provided. +""" +function right_canonical(A::Matrix{<:AbstractTensorMap}, L::Matrix{<:AbstractTensorMap} = initial_C(A); tol = 1e-12, maxiter = 100, kwargs...) + Ar = [permute(A, ((4,2,3), (1,))) for A in A] + Lr = [permute(L, ((2, ), (1,))) for L in L] + + AL, L, λ = left_canonical(Ar, Lr; tol, maxiter, kwargs...) + + R = [permute( L, ((2,), (1, ))) for L in L] + AR = [permute(AL, ((4,2,3), (1,))) for AL in AL] + return R, AR, λ +end + +function initial_FL(AL::Matrix{<:AbstractTensorMap}, ipeps::InfinitePEPS) + T = eltype(ipeps[1]) + FL = [(D = space(ipeps, 5)'; + χ = space(AL, 4)'; + TensorMap(rand, T, χ * D * D', χ)) for (ipeps, AL) in zip(ipeps.A, AL)] + return FL +end + +function initial_FR(AR::Matrix{<:AbstractTensorMap}, ipeps::InfinitePEPS) + T = eltype(ipeps[1]) + FR = [(D = space(ipeps, 3)'; + χ = space(AR, 4)'; + TensorMap(rand, T, χ * D * D', χ)) for (ipeps, AR) in zip(ipeps.A, AR)] + return FR +end + +""" + λL, FL = leftenv(ALu, ALd, O, FL = initial_FL(ALu,O); kwargs...) + +Compute the left environment tensor for MPS A and MPO O, by finding the left fixed point +of ALu - O - ALd contracted along the physical dimension. +``` + ┌── ALuᵢⱼ ── ┌── + │ │ │ +FLᵢⱼ ─ Oᵢⱼ ── = λLᵢⱼ FLᵢⱼ₊₁ + │ │ │ + └── ALdᵢᵣⱼ ─ └── +``` +""" +function leftenv(ALu::Matrix{<:AbstractTensorMap}, + ALd::Matrix{<:AbstractTensorMap}, + ipeps::InfinitePEPS, + FL::Matrix{<:AbstractTensorMap} = initial_FL(ALu,ipeps); + ifobs=false, verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(ipeps) + λL = Zygote.Buffer(zeros(eltype(ipeps[1]), Ni)) + FL′ = Zygote.Buffer(FL) + for i in 1:Ni + ir = ifobs ? Ni + 1 - i : mod1(i + 1, Ni) + λLs, FL1s, info = eigsolve(FLi -> FLmap(FLi, ALu[i,:], ALd[ir,:], ipeps.A[i,:], adjoint.(ipeps.A[i,:])), + FL[i,:], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "leftenv not converged" + λL[i], FL′[i,:] = selectpos(λLs, FL1s, Nj) + end + return copy(λL), copy(FL′) +end + +""" + leftCenv(ALu::Matrix{<:AbstractTensorMap}, + ALd::Matrix{<:AbstractTensorMap}, + L::Matrix{<:AbstractTensorMap} = initial_C(ALu); + ifobs=false, verbosity = Defaults.verbosity, kwargs...) + +Compute the left environment tensor for MPS A, by finding the left fixed point +of ALu - ALd contracted along the physical dimension. +``` + ┌── ALuᵢⱼ ── ┌── + Lᵢⱼ | = λLᵢⱼ Lᵢⱼ₊₁ + └── ALdᵢᵣⱼ ── └── +``` +""" +function leftCenv(ALu::Matrix{<:AbstractTensorMap}, + ALd::Matrix{<:AbstractTensorMap}, + L::Matrix{<:AbstractTensorMap} = initial_C(ALu); + ifobs=false, verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(ALu) + λL = Zygote.Buffer(zeros(eltype(ALu[1]), Ni)) + L′ = Zygote.Buffer(L) + for i in 1:Ni + ir = ifobs ? mod1(Ni - i + 2, Ni) : i + λLs, L1s, info = eigsolve(L -> Lmap(L, ALu[i,:], ALd[ir,:]), + L[i,:], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "leftenv not converged" + λL[i], L′[i,:] = selectpos(λLs, L1s, Nj) + end + return copy(λL), copy(L′) +end + +""" + λR, FR = rightenv(ARu, ARd, M, FR = FRint(ARu,M); kwargs...) + +Compute the right environment tensor for MPS A and MPO M, by finding the left fixed point +of AR - M - conj(AR) contracted along the physical dimension. +``` + ── ARuᵢⱼ ──┐ ──┐ + │ │ │ + ── Mᵢⱼ ──FRᵢⱼ = λRᵢⱼ──FRᵢⱼ₋₁ + │ │ │ + ── ARdᵢᵣⱼ ──┘ ──┘ +``` +""" +function rightenv(ARu::Matrix{<:AbstractTensorMap}, + ARd::Matrix{<:AbstractTensorMap}, + ipeps::InfinitePEPS, + FR::Matrix{<:AbstractTensorMap} = initial_FR(ARu,ipeps); + ifobs=false, ifinline=false,verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(ipeps) + λR = Zygote.Buffer(zeros(eltype(ipeps[1]), Ni)) + FR′ = Zygote.Buffer(FR) + for i in 1:Ni + ir = ifobs ? Ni + 1 - i : mod1(i + 1, Ni) + ifinline && (ir = i) + λRs, FR1s, info = eigsolve(FR -> FRmap(FR, ARu[i,:], ARd[ir,:], ipeps.A[i,:], adjoint.(ipeps.A[i,:])), + FR[i,:], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "rightenv not converged" + λR[i], FR′[i,:] = selectpos(λRs, FR1s, Nj) + end + return copy(λR), copy(FR′) +end + +""" + λR, FR = rightCenv(ARu::Matrix{<:AbstractTensorMap}, + ARd::Matrix{<:AbstractTensorMap}, + R::Matrix{<:AbstractTensorMap} = initial_C(ARu); + kwargs...) + +Compute the right environment tensor for MPS A by finding the left fixed point +of AR - conj(AR) contracted along the physical dimension. +``` + ── ARuᵢⱼ ──┐ ──┐ + | Rᵢⱼ = λRᵢⱼ Rᵢⱼ₋₁ + ── ARdᵢᵣⱼ ──┘ ──┘ +``` +""" +function rightCenv(ARu::Matrix{<:AbstractTensorMap}, + ARd::Matrix{<:AbstractTensorMap}, + R::Matrix{<:AbstractTensorMap} = initial_C(ARu); + ifobs=false, verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(ARu) + λR = Zygote.Buffer(zeros(eltype(ARu[1]), Ni)) + R′ = Zygote.Buffer(R) + for i in 1:Ni + ir = ifobs ? mod1(Ni - i + 2, Ni) : i + λRs, R1s, info = eigsolve(R -> Rmap(R, ARu[i,:], ARd[ir,:]), + R[i,:], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "rightenv not converged" + λR[i], R′[i,:] = selectpos(λRs, R1s, Nj) + end + return copy(λR), copy(R′) +end + +""" + ACenv(AC, FL, M, FR;kwargs...) + +Compute the up environment tensor for MPS `FL`,`FR` and MPO `M`, by finding the up fixed point + of `FL - M - FR` contracted along the physical dimension. +``` +┌─────── ACᵢⱼ ─────┐ +│ │ │ = λACᵢⱼ ┌─── ACᵢ₊₁ⱼ ──┐ +FLᵢⱼ ─── Oᵢⱼ ───── FRᵢⱼ │ │ │ +│ │ │ +``` +""" +function ACenv(AC::Matrix{<:AbstractTensorMap}, + FL::Matrix{<:AbstractTensorMap}, + FR::Matrix{<:AbstractTensorMap}, + ipeps::InfinitePEPS; + verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(ipeps) + λAC = Zygote.Buffer(zeros(eltype(ipeps[1]),Nj)) + AC′ = Zygote.Buffer(AC) + for j in 1:Nj + λACs, ACs, info = eigsolve(AC -> ACmap(AC, FL[:,j], FR[:,j], ipeps.A[:,j], adjoint.(ipeps.A[:,j])), + AC[:,j], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "ACenv Not converged" + λAC[j], AC′[:,j] = selectpos(λACs, ACs, Ni) + end + return copy(λAC), copy(AC′) +end + +""" + λC, C = Cenv(C::Matrix{<:AbstractTensorMap}, + FL::Matrix{<:AbstractTensorMap}, + FR::Matrix{<:AbstractTensorMap}; + kwargs...) = Cenv!(copy(C), FL, FR; kwargs...) + +Compute the up environment tensor for MPS `FL` and `FR`, by finding the up fixed point + of `FL - FR` contracted along the physical dimension. +``` +┌────Cᵢⱼ ───┐ +│ │ = λCᵢⱼ ┌──Cᵢⱼ ─┐ +FLᵢⱼ₊₁ ──── FRᵢⱼ │ │ +│ │ +``` +""" +function Cenv(C::Matrix{<:AbstractTensorMap}, + FL::Matrix{<:AbstractTensorMap}, + FR::Matrix{<:AbstractTensorMap}; + verbosity = Defaults.verbosity, kwargs...) + + Ni, Nj = size(C) + λC = Zygote.Buffer(zeros(eltype(C[1]), Nj)) + C′ = Zygote.Buffer(C) + for j in 1:Nj + jr = mod1(j + 1, Nj) + λCs, Cs, info = eigsolve(C -> Cmap(C, FL[:,jr], FR[:,j]), + C[:,j], 1, :LM; alg_rrule=KrylovKit.GMRES(), maxiter=100, ishermitian = false, kwargs...) + verbosity >= 1 && info.converged == 0 && @warn "Cenv Not converged" + λC[j], C′[:,j] = selectpos(λCs, Cs, Ni) + end + return copy(λC), copy(C′) +end + +function env_norm(F::Matrix{<:AbstractTensorMap}) + Ni, Nj = size(F) + buf = Zygote.Buffer(F) + @inbounds for j in 1:Nj, i in 1:Ni + buf[i,j] = F[i,j]/norm(F[i,j]) + end + return copy(buf) +end + +""" + AL, AR = ACCtoALAR(AC, C) + +QR factorization to get `AL` and `AR` from `AC` and `C` + +```` +──ALᵢⱼ──Cᵢⱼ── = ──ACᵢⱼ── = ──Cᵢ₋₁ⱼ ──ARᵢⱼ── + │ │ │ +```` +""" +function ACCtoALAR(AC::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap}) + AC = env_norm(AC) + C = env_norm( C) + AL, errL = ACCtoAL(AC, C) + AR, errR = ACCtoAR(AC, C) + return AL, AR, errL, errR +end + +function ACCtoAL(AC::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap}) + Ni, Nj = size(AC) + errL = 0.0 + AL = Zygote.Buffer(AC) + @inbounds for j in 1:Nj, i in 1:Ni + QAC, RAC = TensorKit.leftorth(AC[i,j]) + QC, RC = TensorKit.leftorth( C[i,j]) + errL += norm(RAC - RC) + AL[i,j] = QAC * QC' + end + return copy(AL), errL +end + +function ACCtoAR(AC::Matrix{<:AbstractTensorMap}, C::Matrix{<:AbstractTensorMap}) + Ni, Nj = size(AC) + errR = 0.0 + AR = Zygote.Buffer(AC) + @inbounds for j in 1:Nj, i in 1:Ni + jr = mod1(j - 1, Nj) + LAC, QAC = rightorth(_to_tail(AC[i,j])) + LC, QC = rightorth(C[i,jr]) + errR += norm(LAC - LC) + AR[i,j] = _to_front(QC' * QAC) + end + return copy(AR), errR +end \ No newline at end of file diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl deleted file mode 100644 index a33c7e2c..00000000 --- a/src/operators/infinitepepo.jl +++ /dev/null @@ -1,139 +0,0 @@ -""" - struct InfinitePEPO{T<:PEPOTensor} - -Represents an infinite projected entangled-pair operator (PEPO) on a 3D cubic lattice. -""" -struct InfinitePEPO{T<:PEPOTensor} <: AbstractPEPO - A::Array{T,3} - - function InfinitePEPO(A::Array{T,3}) where {T<:PEPOTensor} - # space checks - for (d, w, h) in Tuple.(CartesianIndices(A)) - space(A[d, w, h], 1) == space(A[d, w, _next(h, end)], 2)' || - throw(SpaceMismatch("Physical space at site $((d, w, h)) does not match.")) - space(A[d, w, h], 3) == space(A[_prev(d, end), w, h], 5)' || throw( - SpaceMismatch("North virtual space at site $((d, w, h)) does not match."), - ) - space(A[d, w, h], 4) == space(A[d, _next(w, end), h], 6)' || throw( - SpaceMismatch("East virtual space at site $((d, w, h)) does not match.") - ) - end - return new{T}(A) - end -end - -## Constructors -""" - InfinitePEPO(A::AbstractArray{T, 3}) - -Allow users to pass in an array of tensors. -""" -function InfinitePEPO(A::AbstractArray{T,3}) where {T<:PEPOTensor} - return InfinitePEPO(Array(deepcopy(A))) -end - -""" - InfinitePEPO(f=randn, T=ComplexF64, Pspaces, Nspaces, Espaces) - -Allow users to pass in arrays of spaces. -""" -function InfinitePEPO( - Pspaces::A, Nspaces::A, Espaces::A=Nspaces -) where {A<:AbstractArray{<:ElementarySpace,3}} - return InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) -end -function InfinitePEPO( - f, T, Pspaces::A, Nspaces::A, Espaces::A=Nspaces -) where {A<:AbstractArray{<:ElementarySpace,3}} - size(Pspaces) == size(Nspaces) == size(Espaces) || - throw(ArgumentError("Input spaces should have equal sizes.")) - - Sspaces = adjoint.(circshift(Nspaces, (1, 0, 0))) - Wspaces = adjoint.(circshift(Espaces, (0, -1, 0))) - Ppspaces = adjoint.(circshift(Pspaces, (0, 0, -1))) - - P = map(Pspaces, Ppspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, Pp, N, E, S, W - return TensorMap(f, T, P * Pp ← N * E * S * W) - end - - return InfinitePEPO(P) -end - -function InfinitePEPO( - Pspaces::A, Nspaces::A, Espaces::A=Nspaces -) where {A<:AbstractArray{<:ElementarySpace,2}} - size(Pspaces) == size(Nspaces) == size(Espaces) || - throw(ArgumentError("Input spaces should have equal sizes.")) - - Pspaces = reshape(Pspaces, (size(Pspaces)..., 1)) - Nspaces = reshape(Pspaces, (size(Nspaces)..., 1)) - Espaces = reshape(Pspaces, (size(Espaces)..., 1)) - - return InfinitePEPO(Pspaces, Nspaces, Espaces) -end - -""" - InfinitePEPO(A; unitcell=(1, 1, 1)) - -Create an InfinitePEPO by specifying a tensor and unit cell. -""" -function InfinitePEPO(A::T; unitcell::Tuple{Int,Int,Int}=(1, 1, 1)) where {T<:PEPOTensor} - return InfinitePEPO(fill(A, unitcell)) -end - -""" - InfinitePEPO(f=randn, T=ComplexF64, Pspace, Nspace, [Espace]; unitcell=(1,1,1)) - -Create an InfinitePEPO by specifying its spaces and unit cell. -""" -function InfinitePEPO( - Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int,Int}=(1, 1, 1) -) where {S<:ElementarySpace} - return InfinitePEPO( - randn, - ComplexF64, - fill(Pspace, unitcell), - fill(Nspace, unitcell), - fill(Espace, unitcell), - ) -end -function InfinitePEPO( - f, T, Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int,Int}=(1, 1, 1) -) where {S<:ElementarySpace} - return InfinitePEPO( - f, T, fill(Pspace, unitcell), fill(Nspace, unitcell), fill(Espace, unitcell) - ) -end - -## Shape and size -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) - -## Copy -Base.copy(T::InfinitePEPO) = InfinitePEPO(copy(T.A)) -Base.similar(T::InfinitePEPO) = InfinitePEPO(similar(T.A)) -Base.repeat(T::InfinitePEPO, counts...) = InfinitePEPO(repeat(T.A, counts...)) - -Base.getindex(T::InfinitePEPO, args...) = Base.getindex(T.A, args...) -Base.axes(T::InfinitePEPO, args...) = axes(T.A, args...) -TensorKit.space(T::InfinitePEPO, i, j) = space(T[i, j, end], 1) - -function initializePEPS( - T::InfinitePEPO{<:PEPOTensor{S}}, vspace::S -) where {S<:ElementarySpace} - Pspaces = Array{S,2}(undef, size(T, 1), size(T, 2)) - for (i, j) in product(1:size(T, 1), 1:size(T, 2)) - Pspaces[i, j] = space(T, i, j) - end - Nspaces = repeat([vspace], size(T, 1), size(T, 2)) - Espaces = repeat([vspace], size(T, 1), size(T, 2)) - return InfinitePEPS(Pspaces, Nspaces, Espaces) -end - -# Rotations -Base.rotl90(T::InfinitePEPO) = InfinitePEPO(stack(rotl90, eachslice(T.A; dims=3))) -Base.rotr90(T::InfinitePEPO) = InfinitePEPO(stack(rotr90, eachslice(T.A; dims=3))) -Base.rot180(T::InfinitePEPO) = InfinitePEPO(stack(rot180, eachslice(T.A; dims=3))) diff --git a/src/operators/lattices/squarelattice.jl b/src/operators/lattices/squarelattice.jl index 097a4f7b..c5bc22a9 100644 --- a/src/operators/lattices/squarelattice.jl +++ b/src/operators/lattices/squarelattice.jl @@ -16,6 +16,19 @@ function vertices(lattice::InfiniteSquare) return CartesianIndices((1:(lattice.Nrows), 1:(lattice.Ncols))) end +""" + nearest_neighbours(lattice::InfiniteSquare) + +Return the nearest neighbors of the lattice `lattice`. + +```` + +---*---+ + | + * + | + + +```` +""" function nearest_neighbours(lattice::InfiniteSquare) neighbors = Tuple{CartesianIndex,CartesianIndex}[] for idx in vertices(lattice) @@ -25,6 +38,19 @@ function nearest_neighbours(lattice::InfiniteSquare) return neighbors end +""" + next_nearest_neighbours(lattice::InfiniteSquare) + +Return the next nearest neighbors of the lattice `lattice`. + +```` + +------+ + | \\ ╱ | + | * | + | ╱ \\| + +------+ +```` +""" function next_nearest_neighbours(lattice::InfiniteSquare) neighbors = Tuple{CartesianIndex,CartesianIndex}[] for idx in vertices(lattice) diff --git a/src/operators/transferpepo.jl b/src/operators/transferpepo.jl deleted file mode 100644 index 5215834c..00000000 --- a/src/operators/transferpepo.jl +++ /dev/null @@ -1,232 +0,0 @@ -""" - InfiniteTransferPEPO{T,O} - -Represents an infinite transfer operator corresponding to a single row of a partition -function which corresponds to the expectation value of an `InfinitePEPO` between 'ket' and -'bra' `InfinitePEPS` states. -""" -struct InfiniteTransferPEPO{T,O} - top::PeriodicArray{T,1} - mid::PeriodicArray{O,2} - bot::PeriodicArray{T,1} -end - -InfiniteTransferPEPO(top, mid) = InfiniteTransferPEPO(top, mid, top) - -""" - InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) - -Constructs a transfer operator corresponding to a single row of a partition function -representing the expectation value of `O` for the state `T`. The partition function is first -rotated such that the direction `dir` faces north, after which its `row`th row from the -north is selected. -""" -function InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) - T = rotate_north(T, dir) - O = rotate_north(O, dir) - return InfiniteTransferPEPO(PeriodicArray(T[row, :]), PeriodicArray(O[row, :, :])) -end - -Base.size(transfer::InfiniteTransferPEPO) = size(transfer.top) -Base.size(transfer::InfiniteTransferPEPO, args...) = size(transfer.top, args...) -Base.length(transfer::InfiniteTransferPEPO) = size(transfer, 1) -height(transfer::InfiniteTransferPEPO) = size(transfer.mid, 2) -Base.getindex(O::InfiniteTransferPEPO, i) = (O.top[i], O.bot[i], Tuple(O.mid[i, :])) # TODO: not too sure about this - -Base.iterate(O::InfiniteTransferPEPO, i=1) = i > length(O) ? nothing : (O[i], i + 1) - -function virtual_spaces(O::InfiniteTransferPEPO, i, dir) - return [ - space(O.top[i], dir + 1), - space.(O.mid[i, :], Ref(dir + 2))..., - space(O.bot[i], dir + 1)', - ] -end -north_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, NORTH) -east_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, EAST) -south_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, SOUTH) -west_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, WEST) - -function initializeMPS(O::InfiniteTransferPEPO, virtualspaces::AbstractArray{S,1}) where {S} - return InfiniteMPS([ - TensorMap( - rand, - MPSKit.Defaults.eltype, # should be scalartype of transfer PEPO? - virtualspaces[_prev(i, end)] * prod(adjoint.(north_spaces(O, i))), - virtualspaces[mod1(i, end)], - ) for i in 1:length(O) - ]) -end - -function initializeMPS(O::InfiniteTransferPEPO, χ::Int) - return InfiniteMPS([ - TensorMap( - rand, MPSKit.Defaults.eltype, ℂ^χ * prod(adjoint.(north_spaces(O, i))), ℂ^χ - ) for i in 1:length(O) - ]) -end - -""" - const TransferPEPOMultiline = MPSKit.Multiline{<:InfiniteTransferPEPO} - -Type that represents a multi-line transfer operator, where each line each corresponds to a -row of a partition function encoding the overlap of an `InfinitePEPO` between 'ket' and -'bra' `InfinitePEPS` states. -""" -const TransferPEPOMultiline = MPSKit.Multiline{<:InfiniteTransferPEPO} -Base.convert(::Type{TransferPEPOMultiline}, O::InfiniteTransferPEPO) = MPSKit.Multiline([O]) -Base.getindex(t::TransferPEPOMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) -Base.getindex(t::TransferPEPOMultiline, i::Int, j) = Base.getindex(t.data[i], j) - -""" - TransferPEPOMultiline(T::InfinitePEPS, O::InfinitePEPO, dir) - -Construct a multi-row transfer operator corresponding to the partition function representing -the expectation value of `O` for the state `T`. The partition function is first rotated such -that the direction `dir` faces north. -""" -function TransferPEPOMultiline(T::InfinitePEPS, O::InfinitePEPO, dir) - rowsize = size(T, mod1(dir, 2)) # depends on dir - return MPSKit.Multiline(map(cr -> InfiniteTransferPEPO(T, O, dir, cr), 1:rowsize)) -end - -# specialize simple case -function MPSKit.transfer_left( - GL::GenericMPSTensor{S,4}, - O::Tuple{T,T,Tuple{P}}, - A::GenericMPSTensor{S,4}, - Ā::GenericMPSTensor{S,4}, -) where {S,T<:PEPSTensor,P<:PEPOTensor} - @tensor GL′[-1 -2 -3 -4; -5] := - GL[10 7 4 2; 1] * - conj(Ā[10 11 12 13; -1]) * - O[1][8; 9 -2 11 7] * - O[3][1][5 8; 6 -3 12 4] * - conj(O[2][5; 3 -4 13 2]) * - A[1 9 6 3; -5] -end - -# general case -function MPSKit.transfer_left( - GL::GenericMPSTensor{S,N}, - O::Tuple{T,T,Tuple{Vararg{P,H}}}, - A::GenericMPSTensor{S,N}, - Ā::GenericMPSTensor{S,N}, -) where {S,T<:PEPSTensor,P<:PEPOTensor,N,H} - # sanity check - @assert H == N - 3 - - # collect tensors in convenient order: env, above, below, top, mid, bot - tensors = [GL, A, Ā, O[1], O[3]..., O[2]] - - # contraction order: GL, A, top, mid..., bot, Ā - - # number of contracted legs for full top-mid-bot stack - nlegs_tmb = 5 + 3 * H - - # assign and collect all contraction indices - indicesGL = [2 + nlegs_tmb, 2, ((1:3:((H + 1) * 3)) .+ 3)..., 1] - indicesA = [1, 3, ((1:3:((H + 1) * 3)) .+ 4)..., -(N + 1)] - indicesĀ = [((1:N) .+ (1 + nlegs_tmb))..., -1] - indicesTop = [6, 3, -2, 3 + nlegs_tmb, 2] - indicesBot = [1 + nlegs_tmb, nlegs_tmb, -N, 4 + H + nlegs_tmb, nlegs_tmb - 1] - indicesMid = Vector{Vector{Int}}(undef, H) - for h in 1:H - indicesMid[h] = [ - 3 + 3 * (h + 1), 3 + 3 * h, 2 + 3 * h, -(2 + h), 3 + h + nlegs_tmb, 1 + 3 * h - ] - end - indices = [indicesGL, indicesA, indicesĀ, indicesTop, indicesMid..., indicesBot] - - # record conjflags - conjlist = [false, false, true, false, repeat([false], H)..., true] - - # perform contraction, permute to restore partition - GL′ = permute(ncon(tensors, indices, conjlist), (Tuple(1:N), (N + 1,))) - - return GL′ -end - -# specialize simple case -function MPSKit.transfer_right( - GR::GenericMPSTensor{S,4}, - O::Tuple{T,T,Tuple{P}}, - A::GenericMPSTensor{S,4}, - Ā::GenericMPSTensor{S,4}, -) where {S,T<:PEPSTensor,P<:PEPOTensor} - return @tensor GR′[-1 -2 -3 -4; -5] := - GR[10 7 4 2; 1] * - conj(Ā[-5 9 6 3; 1]) * - O[1][8; 11 7 9 -2] * - O[3][1][5 8; 12 4 6 -3] * - conj(O[2][5; 13 2 3 -4]) * - A[-1 11 12 13; 10] -end - -# general case -function MPSKit.transfer_right( - GR::GenericMPSTensor{S,N}, - O::Tuple{T,T,Tuple{Vararg{P,H}}}, - A::GenericMPSTensor{S,N}, - Ā::GenericMPSTensor{S,N}, -) where {S,T<:PEPSTensor,P<:PEPOTensor,N,H} - # sanity check - @assert H == N - 3 - - # collect tensors in convenient order: env, above, below, top, mid, bot - tensors = [GR, A, Ā, O[1], O[3]..., O[2]] - - # contraction order: GR, A, top, mid..., bot, Ā - - # number of contracted legs for full top-mid-bot stack - nlegs_tmb = 5 + 3 * H - - # assign and collect all contraction indices - indicesGR = [1, 2, ((1:3:((H + 1) * 3)) .+ 3)..., 2 + nlegs_tmb] - indicesA = [-1, 3, ((1:3:((H + 1) * 3)) .+ 4)..., 1] - indicesĀ = [-(N + 1), ((2:N) .+ (1 + nlegs_tmb))..., 2 + nlegs_tmb] - indicesTop = [6, 3, 2, 3 + nlegs_tmb, -2] - indicesBot = [1 + nlegs_tmb, nlegs_tmb, nlegs_tmb - 1, 4 + H + nlegs_tmb, -N] - indicesMid = Vector{Vector{Int}}(undef, H) - for h in 1:H - indicesMid[h] = [ - 3 + 3 * (h + 1), 3 + 3 * h, 2 + 3 * h, 1 + 3 * h, 3 + h + nlegs_tmb, -(2 + h) - ] - end - indices = [indicesGR, indicesA, indicesĀ, indicesTop, indicesMid..., indicesBot] - - # record conjflags - conjlist = [false, false, true, false, repeat([false], H)..., true] - - # perform contraction, permute to restore partition - GR′ = permute(ncon(tensors, indices, conjlist), (Tuple(1:N), (N + 1,))) - - return GR′ -end - -function MPSKit.expectation_value(st::InfiniteMPS, transfer::InfiniteTransferPEPO) - return expectation_value( - convert(MPSMultiline, st), convert(TransferPEPOMultiline, transfer) - ) -end -function MPSKit.expectation_value(st::MPSMultiline, mpo::TransferPEPOMultiline) - return expectation_value(st, environments(st, mpo)) -end -function MPSKit.expectation_value( - st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A} -) where {H<:TransferPEPOMultiline,V,S,A} - return expectation_value(st, ca.opp, ca) -end -function MPSKit.expectation_value( - st::MPSMultiline, opp::TransferPEPOMultiline, ca::MPSKit.PerMPOInfEnv -) - return prod(product(1:size(st, 1), 1:size(st, 2))) do (i, j) - O_ij = opp[i, j] - N = height(opp[1]) + 4 - # just reuse left environment contraction - GL´ = transfer_left(leftenv(ca, i, j, st), O_ij, st.AC[i, j], st.AC[i + 1, j]) - return TensorOperations.tensorscalar( - ncon([GL´, rightenv(ca, i, j, st)], [[N, (2:(N - 1))..., 1], [(1:N)...]]) - ) - end -end diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl deleted file mode 100644 index 4706c293..00000000 --- a/src/operators/transferpeps.jl +++ /dev/null @@ -1,179 +0,0 @@ -""" - InfiniteTransferPEPS{T} - -Represents an infinite transfer operator corresponding to a single row of a partition -function which corresponds to the overlap between 'ket' and 'bra' `InfinitePEPS` states. -""" -struct InfiniteTransferPEPS{T} - top::PeriodicArray{T,1} - bot::PeriodicArray{T,1} -end - -InfiniteTransferPEPS(top) = InfiniteTransferPEPS(top, top) - -""" - InfiniteTransferPEPS(T::InfinitePEPS, dir, row) - -Constructs a transfer operator corresponding to a single row of a partition function -representing the norm of the state `T`. The partition function is first rotated such that -the direction `dir` faces north, after which its `row`th row from the north is selected. -""" -function InfiniteTransferPEPS(T::InfinitePEPS, dir, row) - T = rotate_north(T, dir) - return InfiniteTransferPEPS(PeriodicArray(T[row, :])) -end - -Base.size(transfer::InfiniteTransferPEPS) = size(transfer.top) -Base.size(transfer::InfiniteTransferPEPS, args...) = size(transfer.top, args...) -Base.length(transfer::InfiniteTransferPEPS) = size(transfer, 1) -Base.getindex(O::InfiniteTransferPEPS, i) = (O.top[i], O.bot[i]) - -Base.iterate(O::InfiniteTransferPEPS, i=1) = i > length(O) ? nothing : (O[i], i + 1) - -import MPSKit.GenericMPSTensor - -""" - const TransferPEPSMultiline = MPSKit.Multiline{<:InfiniteTransferPEPS} - -Type that represents a multi-line transfer operator, where each line each corresponds to a -row of a partition function encoding the overlap between 'ket' and 'bra' `InfinitePEPS` -states. -""" -const TransferPEPSMultiline = MPSKit.Multiline{<:InfiniteTransferPEPS} -Base.convert(::Type{TransferPEPSMultiline}, O::InfiniteTransferPEPS) = MPSKit.Multiline([O]) -Base.getindex(t::TransferPEPSMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) -Base.getindex(t::TransferPEPSMultiline, i::Int, j) = Base.getindex(t.data[i], j) - -""" - TransferPEPSMultiline(T::InfinitePEPS, dir) - -Construct a multi-row transfer operator corresponding to the partition function representing -the norm of the state `T`. The partition function is first rotated such -that the direction `dir` faces north. -""" -function TransferPEPSMultiline(T::InfinitePEPS, dir) - rowsize = size(T, mod1(dir, 2)) # depends on dir - return MPSKit.Multiline(map(cr -> InfiniteTransferPEPS(T, dir, cr), 1:rowsize)) -end - -""" - initializeMPS( - O::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, - virtualspaces::AbstractArray{<:ElementarySpace,1} - ) - initializeMPS( - O::Union{TransferPEPSMultiline,TransferPEPOMultiline}, - virtualspaces::AbstractArray{<:ElementarySpace,2} - ) - -Inialize a boundary MPS for the transfer operator `O` by specifying an array of virtual -spaces consistent with the unit cell. -""" -function initializeMPS(O::InfiniteTransferPEPS, virtualspaces::AbstractArray{S,1}) where {S} - return InfiniteMPS([ - TensorMap( - rand, - MPSKit.Defaults.eltype, # should be scalartype of transfer PEPS? - virtualspaces[_prev(i, end)] * space(O.top[i], 2)' * space(O.bot[i], 2), - virtualspaces[mod1(i, end)], - ) for i in 1:length(O) - ]) -end -function initializeMPS(O::InfiniteTransferPEPS, χ::Int) - return InfiniteMPS([ - TensorMap( - rand, - MPSKit.Defaults.eltype, - ℂ^χ * space(O.top[i], 2)' * space(O.bot[i], 2), - ℂ^χ, - ) for i in 1:length(O) - ]) -end -function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,2}) where {S} - mpss = map(cr -> initializeMPS(O[cr], virtualspaces[cr, :]), 1:size(O, 1)) - return MPSKit.Multiline(mpss) -end -function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,1}) where {S} - return initializeMPS(O, repeat(virtualspaces, length(O), 1)) -end -function initializeMPS(O::MPSKit.Multiline, V::ElementarySpace) - return initializeMPS(O, repeat([V], length(O), length(O[1]))) -end -function initializeMPS(O::MPSKit.Multiline, χ::Int) - return initializeMPS(O, repeat([ℂ^χ], length(O), length(O[1]))) -end - -function MPSKit.transfer_left( - GL::GenericMPSTensor{S,3}, - O::NTuple{2,PEPSTensor}, - A::GenericMPSTensor{S,3}, - Ā::GenericMPSTensor{S,3}, -) where {S} - return @tensor GL′[-1 -2 -3; -4] := - GL[1 2 4; 7] * - conj(Ā[1 3 6; -1]) * - O[1][5; 8 -2 3 2] * - conj(O[2][5; 9 -3 6 4]) * - A[7 8 9; -4] -end - -function MPSKit.transfer_right( - GR::GenericMPSTensor{S,3}, - O::NTuple{2,PEPSTensor}, - A::GenericMPSTensor{S,3}, - Ā::GenericMPSTensor{S,3}, -) where {S} - return @tensor GR′[-1 -2 -3; -4] := - GR[7 6 2; 1] * - conj(Ā[-4 4 3; 1]) * - O[1][5; 9 6 4 -2] * - conj(O[2][5; 8 2 3 -3]) * - A[-1 9 8 7] -end - -@doc """ - MPSKit.expectation_value(st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}) - MPSKit.expectation_value(st::MPSMultiline, op::Union{TransferPEPSMultiline,TransferPEPOMultiline}) - -Compute expectation value of the transfer operator `op` for the state `st` for each site in -the unit cell. -""" MPSKit.expectation_value(st, op) - -function MPSKit.expectation_value(st::InfiniteMPS, transfer::InfiniteTransferPEPS) - return expectation_value( - convert(MPSMultiline, st), convert(TransferPEPSMultiline, transfer) - ) -end -function MPSKit.expectation_value(st::MPSMultiline, mpo::TransferPEPSMultiline) - return expectation_value(st, environments(st, mpo)) -end -function MPSKit.expectation_value( - st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A} -) where {H<:TransferPEPSMultiline,V,S,A} - return expectation_value(st, ca.opp, ca) -end -function MPSKit.expectation_value( - st::MPSMultiline, opp::TransferPEPSMultiline, ca::MPSKit.PerMPOInfEnv -) - return prod(product(1:size(st, 1), 1:size(st, 2))) do (i, j) - O_ij = opp[i, j] - return @tensor leftenv(ca, i, j, st)[1 2 4; 7] * - conj(st.AC[i + 1, j][1 3 6; 13]) * - O_ij[1][5; 8 11 3 2] * - conj(O_ij[2][5; 9 12 6 4]) * - st.AC[i, j][7 8 9; 10] * - rightenv(ca, i, j, st)[10 11 12; 13] - end -end - -@doc """ - MPSKit.leading_boundary( - st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, alg, [envs] - ) - MPSKit.leading_boundary( - st::MPSMulitline, op::Union{TransferPEPSMultiline,TransferPEPOMultiline}, alg, [envs] - ) - -Approximate the leading boundary MPS eigenvector for the transfer operator `op` using `st` -as initial guess. -""" MPSKit.leading_boundary(st, op, alg) diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 93a27142..c6694c11 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -1,6 +1,14 @@ """ const PEPSTensor{S} +```` + N + ↓ +W ← + ← E + ↓ ⬊ + S P + +```` Default type for PEPS tensors with a single physical index, and 4 virtual indices, conventionally ordered as: ``T : P ← N ⊗ E ⊗ S ⊗ W``. Here, ``P``, ``N``, ``E``, ``S`` and ``W`` denote the physics, north, east, south and west spaces, respectively. diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 0a7d63aa..e6236a62 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -134,6 +134,10 @@ Base.:/(ψ::InfinitePEPS, α::Number) = InfinitePEPS(ψ.A / α) LinearAlgebra.dot(ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) = dot(ψ₁.A, ψ₂.A) LinearAlgebra.norm(ψ::InfinitePEPS) = norm(ψ.A) +## not a nice solution, but it works +Base.:+(ipeps::InfinitePEPS, x::NamedTuple) = InfinitePEPS(ipeps.A .+ x.A) +Base.:+(x::NamedTuple, ipeps::InfinitePEPS) = InfinitePEPS(ipeps.A .+ x.A) + ## (Approximate) equality function Base.:(==)(ψ₁::InfinitePEPS, ψ₂::InfinitePEPS) return all(zip(ψ₁.A, ψ₂.A)) do (p₁, p₂) diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl new file mode 100644 index 00000000..5a5b120c --- /dev/null +++ b/src/utility/defaults.jl @@ -0,0 +1,54 @@ +""" + module Defaults + +Module containing default values that represent typical algorithm parameters. + +- `contr_maxiter`: Maximum number of iterations for the contraction algorithm. +- `contr_miniter`: Minimum number of iterations for the contraction algorithm. +- `contr_tol`: Tolerance for the contraction algorithm. +- `fpgrad_maxiter`: Maximum number of iterations for the fixed-point gradient algorithm. +- `fpgrad_tol`: Tolerance for the fixed-point gradient algorithm. +- `verbosity`: Level of verbosity for the algorithm. +- `contractionscheme`: Scheme for contracting the environment. +- `reuse_env`: Whether to reuse the environment. +- `trscheme`: Truncation scheme. +- `iterscheme`: Iteration scheme. +- `fwd_alg`: Forward algorithm for the SVD. +- `rrule_alg`: Rule algorithm for the SVD. +- `svd_alg`: SVD algorithm. +- `optimizer`: Optimizer for the algorithm. +- `gradient_linsolver`: Linear solver for the gradient. +- `gradient_alg`: Algorithm for the gradient. +- `_finalize`: Function to finalize the algorithm. +""" +module Defaults + const VERBOSE_NONE = 0 + const VERBOSE_WARN = 1 + const VERBOSE_CONV = 2 + const VERBOSE_ITER = 3 + const VERBOSE_ALL = 4 + + using TensorKit, KrylovKit, OptimKit + using PEPSKit: LinSolver, FixedSpaceTruncation, SVDAdjoint + const eltype = ComplexF64 + const contr_maxiter = 100 + const contr_miniter = 4 + const contr_tol = 1e-8 + const fpgrad_maxiter = 30 + const fpgrad_tol = 1e-6 + const verbosity = VERBOSE_ITER + const ctmrgscheme = :simultaneous + const reuse_env = true + const trscheme = FixedSpaceTruncation() + const iterscheme = :fixed + const fwd_alg = TensorKit.SVD() + const rrule_alg = GMRES(; tol=1e1contr_tol) + const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg) + const optimizer = LBFGS(32; maxiter=100, gradtol=1e-8, verbosity=2) + const gradient_linsolver = KrylovKit.BiCGStab(; + maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol + ) + const gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme) + + _finalize(iter, state, opp, envs) = (state, envs) +end \ No newline at end of file diff --git a/src/utility/loginfo.jl b/src/utility/loginfo.jl new file mode 100644 index 00000000..4847e00e --- /dev/null +++ b/src/utility/loginfo.jl @@ -0,0 +1,9 @@ +contr_loginit!(log, η, N) = @infov 2 loginit!(log, η, N) +contr_logiter!(log, iter, η, N) = @infov 3 logiter!(log, iter, η, N) +contr_logfinish!(log, iter, η, N) = @infov 2 logfinish!(log, iter, η, N) +contr_logcancel!(log, iter, η, N) = @warnv 1 logcancel!(log, iter, η, N) + +@non_differentiable contr_loginit!(args...) +@non_differentiable contr_logiter!(args...) +@non_differentiable contr_logfinish!(args...) +@non_differentiable contr_logcancel!(args...) \ No newline at end of file diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl index 5bec6a96..8e3bc5f6 100644 --- a/src/utility/symmetrization.jl +++ b/src/utility/symmetrization.jl @@ -68,8 +68,18 @@ function _fit_spaces( end return y end + +function ChainRulesCore.rrule(::typeof(_fit_spaces), y::AbstractTensorMap, x::AbstractTensorMap) + function pullback(Δ) + return NoTangent(), _fit_spaces(Δ, y), NoTangent() + end + return _fit_spaces(y, x), pullback +end + _fit_spaces(y::InfinitePEPS, x::InfinitePEPS) = InfinitePEPS(map(_fit_spaces, y.A, x.A)) + + function herm_depth_inv(x::Union{PEPSTensor,PEPOTensor}) return 0.5 * (x + _fit_spaces(herm_depth(x), x)) end diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl deleted file mode 100644 index a161c073..00000000 --- a/test/boundarymps/vumps.jl +++ /dev/null @@ -1,69 +0,0 @@ -using Test -using Random -using PEPSKit -using TensorKit -using MPSKit -using LinearAlgebra - -Random.seed!(29384293742893) - -const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitian=false)) -@testset "(1, 1) PEPS" begin - psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) - - T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1) - mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) - - 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)) - N´ = abs(norm(psi, ctm)) - - @test N ≈ N´ atol = 1e-3 -end - -@testset "(2, 2) PEPS" begin - psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2)) - T = PEPSKit.TransferPEPSMultiline(psi, 1) - - mps = PEPSKit.initializeMPS(T, fill(ComplexSpace(20), 2, 2)) - 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)) - N´ = abs(norm(psi, ctm)) - - @test N ≈ N´ rtol = 1e-2 -end - -@testset "PEPO runthrough" begin - function ising_pepo(beta; unitcell=(1, 1, 1)) - t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)] - q = sqrt(t) - - O = zeros(2, 2, 2, 2, 2, 2) - O[1, 1, 1, 1, 1, 1] = 1 - O[2, 2, 2, 2, 2, 2] = 1 - @tensor o[-1 -2; -3 -4 -5 -6] := - O[1 2; 3 4 5 6] * - q[-1; 1] * - q[-2; 2] * - q[-3; 3] * - q[-4; 4] * - q[-5; 5] * - q[-6; 6] - - O = TensorMap(o, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') - - return InfinitePEPO(O; unitcell) - end - - psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) - O = ising_pepo(1) - T = InfiniteTransferPEPO(psi, O, 1, 1) - - mps = PEPSKit.initializeMPS(T, [ComplexSpace(10)]) - mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) - f = abs(prod(expectation_value(mps, T))) -end diff --git a/test/heisenberg.jl b/test/heisenberg_ctmrg.jl similarity index 91% rename from test/heisenberg.jl rename to test/heisenberg_ctmrg.jl index 4fa67362..11669eba 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg_ctmrg.jl @@ -18,7 +18,7 @@ ctm_alg = CTMRG(; ) opt_alg = PEPSOptimize(; boundary_alg=ctm_alg, - optimizer=LBFGS(4; maxiter=100, gradtol=1e-3, verbosity=2), + optimizer=LBFGS(4; maxiter=0, gradtol=1e-3, verbosity=2), gradient_alg=LinSolver(; solver=GMRES(; tol=1e-6), iterscheme=:fixed), reuse_env=true, ) @@ -30,7 +30,7 @@ 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) +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 diff --git a/test/heisenberg_vumps.jl b/test/heisenberg_vumps.jl new file mode 100644 index 00000000..c90b9f1f --- /dev/null +++ b/test/heisenberg_vumps.jl @@ -0,0 +1,123 @@ +using Test +using Random +using PEPSKit +using TensorKit +using KrylovKit +using OptimKit +using Printf + +@testset "heisenberg_XYZ 1x1 unitcell C4 symmetry" begin + Random.seed!(100) + # initialize parameters + χbond = 2 + χenv = 10 + + # initialize states + H = heisenberg_XYZ(InfiniteSquare()) + psi_init = InfinitePEPS(2, χbond; unitcell=(1, 1)) + + # find fixedpoint one-site ctmrg + boundary_alg = VUMPS( + ifupdown=false, + tol=1e-10, + miniter=3, + maxiter=10, + verbosity=1 + ) + opt_alg = PEPSOptimize(; + boundary_alg, + optimizer=LBFGS(20; maxiter=100, gradtol=1e-6, verbosity=2), + gradient_alg=nothing, + reuse_env=true + ) + + t = time() + function finalize!(x, f, g, numiter) + print(@sprintf("%.3f sec", time()-t)) + return x, f, g + end + + env_init = VUMPSRuntime(psi_init, χenv, boundary_alg) + result = fixedpoint(psi_init, H, opt_alg, env_init; + symmetrization=RotateReflect(), + finalize! + ); + @test result.E ≈ -0.66023 atol = 1e-4 +end + +@testset "heisenberg_XYZ 1x1 unitcell without C4 symmetry" begin + Random.seed!(100) + # initialize parameters + χbond = 2 + χenv = 10 + + # initialize states + H = heisenberg_XYZ(InfiniteSquare()) + psi_init = InfinitePEPS(2, χbond; unitcell=(1, 1)) + + # find fixedpoint one-site ctmrg + boundary_alg = VUMPS( + ifupdown=true, + tol=1e-10, + miniter=1, + maxiter=10, + verbosity=1 + ) + opt_alg = PEPSOptimize(; + boundary_alg, + optimizer=LBFGS(20; maxiter=100, gradtol=1e-6, verbosity=2), + gradient_alg=nothing, + reuse_env=true + ) + + t = time() + function finalize!(x, f, g, numiter) + print(@sprintf("%.3f sec", time()-t)) + return x, f, g + end + + env_init = VUMPSRuntime(psi_init, χenv, boundary_alg) + result = fixedpoint(psi_init, H, opt_alg, env_init; + finalize! + ); + @test result.E ≈ -0.66251 atol = 1e-4 +end + +@testset "heisenberg_XYZ 2x2 unitcell without C4 symmetry" begin + Random.seed!(42) + # initialize parameters + χbond = 2 + χenv = 16 + + # initialize states + H = heisenberg_XYZ(InfiniteSquare()) + psi_init = InfinitePEPS(2, χbond; unitcell=(2, 2)) + + # find fixedpoint one-site ctmrg + boundary_alg = VUMPS( + ifupdown=true, + tol=1e-10, + miniter=1, + maxiter=10, + verbosity=2 + ) + opt_alg = PEPSOptimize(; + boundary_alg, + optimizer=LBFGS(20; maxiter=100, gradtol=1e-6, verbosity=2), + gradient_alg=nothing, + reuse_env=true + ) + + t = time() + function finalize!(x, f, g, numiter) + print(@sprintf("%.3f sec", time()-t)) + return x, f, g + end + + env_init = VUMPSRuntime(psi_init, χenv, boundary_alg) + result = fixedpoint(psi_init, H, opt_alg, env_init; + finalize! + ); + @show result.E + # @test result.E ≈ -0.66251 atol = 1e-4 +end diff --git a/test/vumps/gradparts.jl b/test/vumps/gradparts.jl new file mode 100644 index 00000000..7d619a12 --- /dev/null +++ b/test/vumps/gradparts.jl @@ -0,0 +1,193 @@ +using Test +using Random +using PEPSKit +using PEPSKit: initial_A, left_canonical, right_canonical, leftenv, rightenv, ACenv, Cenv, LRtoC, ALCtoAC, ACCtoALAR +using TensorKit +using Zygote +using LinearAlgebra + +begin "test utility" + ds = [ℂ^2] + Ds = [ℂ^2] + χs = [ℂ^5] + + function num_grad(f, K::Number; δ::Real=1e-5) + if eltype(K) == ComplexF64 + (f(K + δ / 2) - f(K - δ / 2)) / δ + + (f(K + δ / 2 * 1.0im) - f(K - δ / 2 * 1.0im)) / δ * 1.0im + else + (f(K + δ / 2) - f(K - δ / 2)) / δ + end + end +end + +@testset "InfinitePEPS for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipepes = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + function f(β) + norm(β * ipepes) + end + @test Zygote.gradient(f, 1.0)[1] ≈ num_grad(f, 1.0) +end + +@testset "leftenv and rightenv for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs), ifobs in [true, false] + Random.seed!(50) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, FL = leftenv(AL, adjoint.(AL), ipeps; ifobs) + λR, FR = rightenv(AR, adjoint.(AR), ipeps; ifobs) + + S1 = TensorMap(randn, ComplexF64, χ*D'*D*χ' ← χ*D'*D*χ') + S2 = TensorMap(randn, ComplexF64, χ*D*D'*χ' ← χ*D*D'*χ') + + function foo1(β) + ipeps = β * ipeps + + _, FL = leftenv(AL, adjoint.(AL), ipeps; ifobs) + + tol = [(@tensor conj(FL[1 2 3 4]) * S1[1 2 3 4; 5 6 7 8] * FL[5 6 7 8]) / dot(FL, FL) for FL in FL] + return norm(tol) + end + @test Zygote.gradient(foo1, 1.0)[1] ≈ num_grad(foo1, 1.0) atol = 1e-8 + + function foo2(β) + ipeps = β * ipeps + + _, FR = rightenv(AR, adjoint.(AR), ipeps; ifobs) + + tol = [(@tensor conj(FR[1 2 3 4]) * S2[1 2 3 4; 5 6 7 8] * FR[5 6 7 8]) / dot(FR, FR) for FR in FR] + return norm(tol) + end + @test Zygote.gradient(foo2, 1.0)[1] ≈ num_grad(foo2, 1.0) atol = 1e-8 +end + +@testset "ACenv and Cenv for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(50) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + _, FL = leftenv(AL, adjoint.(AL), ipeps) + _, FR = rightenv(AR, adjoint.(AR), ipeps) + C = LRtoC(L, R) + AC = ALCtoAC(AL, C) + + S1 = TensorMap(randn, ComplexF64, χ*D*D'*χ' ← χ*D*D'*χ') + S2 = TensorMap(randn, ComplexF64, χ*χ' ← χ*χ') + + function foo1(β) + ipeps = β * ipeps + _, AC = ACenv(AC, FL, FR, ipeps) + + tol = [(@tensor conj(AC[1 2 3 4]) * S1[1 2 3 4; 5 6 7 8] * AC[5 6 7 8]) / dot(AC, AC) for AC in AC] + return norm(tol) + end + @test Zygote.gradient(foo1, 1.0)[1] ≈ num_grad(foo1, 1.0) atol = 1e-8 + + function foo2(β) + _, C = Cenv(C, β*FL, β*FR) + + tol = [(@tensor d = conj(C[1 2]) * S2[1 2; 3 4] * C[3 4]) / dot(C, C) for C in C] + return norm(tol) + end + @test Zygote.gradient(foo2, 1.0)[1] ≈ num_grad(foo2, 1.0) atol = 1e-8 +end + +@testset "ACCtoALAR for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, FL = leftenv(AL, adjoint.(AL), ipeps) + λR, FR = rightenv(AR, adjoint.(AR), ipeps) + + C = LRtoC(L, R) + AC = ALCtoAC(AL, C) + + λAC, AC = ACenv(AC, FL, FR, ipeps) + λC, C = Cenv( C, FL, FR) + + S = TensorMap(randn, ComplexF64, χ*D*D'*χ' ← χ*D*D'*χ') + function foo1(β) + ipeps = β * ipeps + _, AC = ACenv(AC, FL, FR, ipeps) + AL, AR, _, _ = ACCtoALAR(AC, C) + tol1 = [(@tensor conj(AL[1 2 3 4]) * S[1 2 3 4; 5 6 7 8] * AL[5 6 7 8]) / dot(AL, AL) for AL in AL] + tol2 = [(@tensor conj(AR[1 2 3 4]) * S[1 2 3 4; 5 6 7 8] * AR[5 6 7 8]) / dot(AR, AR) for AR in AR] + return norm(tol1 + tol2) + end + @test Zygote.gradient(foo1, 1.0)[1] ≈ num_grad(foo1, 1.0) atol = 1e-8 +end + +@testset "ad vumps iPEPS one side for unitcell $Ni x $Nj" for Ni in 1:1, Nj in 1:1, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(50) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + ipeps = symmetrize!(ipeps, RotateReflect()) + + alg = VUMPS(maxiter=100, verbosity=2, ifupdown=false) + rt = leading_boundary(VUMPSRuntime(ipeps, χ, alg), ipeps, alg) + + function foo1(ipeps) + rt = leading_boundary(rt, ipeps, alg) + env = VUMPSEnv(rt, ipeps) + Z = abs(norm(ipeps, env)) + return Z + end + + function foo2(ipeps) + ctm = leading_boundary(CTMRGEnv(ipeps, χ), ipeps, CTMRG(; verbosity=2)) + Z = abs(norm(ipeps, ctm))^(1/Ni/Nj) + return Z + end + + @show norm(foo1(ipeps) - foo2(ipeps)) < 1e-8 + @test norm(Zygote.gradient(foo1, ipeps)[1].A - Zygote.gradient(foo2, ipeps)[1].A) < 1e-8 +end + +@testset "_fit_spaces" for (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + A = TensorMap(randn, ComplexF64, d ← D*D*D'*D') + B = TensorMap(randn, ComplexF64, d ← D*D*D'*D) + + function f(β) + C = PEPSKit._fit_spaces(β * B, A) + return norm(C) + end + + @test Zygote.gradient(f, 1.0)[1] ≈ num_grad(f, 1.0) +end + +@testset "ad vumps iPEPS two side for unitcell $Ni x $Nj" for Ni in 2:2, Nj in 2:2, (d, D, χ) in zip(ds, Ds, [10]) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + alg = VUMPS(maxiter=100, verbosity=3, ifupdown=true) + rt = leading_boundary(VUMPSRuntime(ipeps, χ, alg), ipeps, alg) + + function foo1(ipeps) + alg = VUMPS(maxiter=2, verbosity=3, ifupdown=true) + rt = leading_boundary(rt, ipeps, alg) + env = VUMPSEnv(rt, ipeps) + Z = abs(norm(ipeps, env)) + return Z + end + Zygote.gradient(foo1, ipeps)[1] + + function foo2(ipeps) + ctm = leading_boundary(CTMRGEnv(ipeps, χ), ipeps, CTMRG(; verbosity=2)) + Z = abs(norm(ipeps, ctm))^(1/Ni/Nj) + return Z + end + + # @test norm(foo1(ipeps) - foo2(ipeps)) < 1e-6 + @show norm(Zygote.gradient(foo1, ipeps)[1].A - Zygote.gradient(foo2, ipeps)[1].A) + # @test norm(Zygote.gradient(foo1, ipeps)[1].A - Zygote.gradient(foo2, ipeps)[1].A) < 1e-5 +end \ No newline at end of file diff --git a/test/vumps/vumps.jl b/test/vumps/vumps.jl new file mode 100644 index 00000000..5bc03862 --- /dev/null +++ b/test/vumps/vumps.jl @@ -0,0 +1,78 @@ +using Test +using Random +using PEPSKit +using PEPSKit: nearest_neighbour_energy +using TensorKit +using LinearAlgebra + +begin "test utility" + ds = [ℂ^2] + Ds = [ℂ^2] + χs = [ℂ^20] +end + +@testset "initialize environments for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + rt = VUMPSRuntime(ipeps, χ) + @test rt isa VUMPSRuntime +end + +@testset "vumps one side runtime for unitcell $Ni x $Nj" for Ni in 1:1, Nj in 1:1, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(100) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + ipeps = symmetrize!(ipeps, RotateReflect()) + + alg = VUMPS(maxiter=100, verbosity=2, ifupdown=false) + + rt = leading_boundary(VUMPSRuntime(ipeps, χ, alg), ipeps, alg) + env = VUMPSEnv(rt, ipeps) + @test env isa VUMPSEnv + + Z = abs(norm(ipeps, env)) + + ctm = leading_boundary(CTMRGEnv(ipeps, χ), ipeps, CTMRG(; verbosity=2)) + Z′ = abs(norm(ipeps, ctm))^(1/Ni/Nj) + + @test Z ≈ Z′ rtol = 1e-12 +end + +@testset "vumps one side runtime energy for unitcell $Ni x $Nj" for Ni in 1:1, Nj in 1:1, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(100) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + ipeps = symmetrize!(ipeps, RotateReflect()) + + alg = VUMPS(maxiter=100, verbosity=2, ifupdown=false) + + rt = leading_boundary(VUMPSRuntime(ipeps, χ, alg), ipeps, alg) + env = VUMPSEnv(rt, ipeps) + H = heisenberg_XYZ(InfiniteSquare()) + H = H.terms[1].second + # Hh, Hv = H.terms[1] + + @show nearest_neighbour_energy(ipeps, H, H, env) + # Z = abs(norm(ipeps, env)) + + # ctm = leading_boundary(CTMRGEnv(ipeps, χ), ipeps, CTMRG(; verbosity=2)) + # Z′ = abs(norm(ipeps, ctm))^(1/Ni/Nj) + + # @test Z ≈ Z′ rtol = 1e-12 +end + +@testset "vumps two side runtime for unitcell $Ni x $Nj" for Ni in 1:3, Nj in 1:3, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + alg = PEPSKit.VUMPS(maxiter=100, verbosity=2, ifupdown=true) + + rt = leading_boundary(VUMPSRuntime(ipeps, χ, alg), ipeps, alg) + env = VUMPSEnv(rt, ipeps) + @test env isa VUMPSEnv + + Z = abs(norm(ipeps, env)) + + ctm = leading_boundary(CTMRGEnv(ipeps, χ), ipeps, CTMRG(; verbosity=2)) + Z′ = abs(norm(ipeps, ctm))^(1/Ni/Nj) + + @test Z ≈ Z′ rtol = 1e-8 +end \ No newline at end of file diff --git a/test/vumps/vumps_environment.jl b/test/vumps/vumps_environment.jl new file mode 100644 index 00000000..88539093 --- /dev/null +++ b/test/vumps/vumps_environment.jl @@ -0,0 +1,175 @@ +using Test +using Random +using PEPSKit +using PEPSKit: initial_A, initial_C, initial_FL, initial_FR +using PEPSKit: ρmap, getL!, getAL, getLsped, _to_tail, _to_front, left_canonical, right_canonical +using PEPSKit: leftenv, FLmap, rightenv, FRmap, ACenv, ACmap, Cenv, Cmap, leftCenv, Lmap, rightCenv, Rmap +using PEPSKit: LRtoC, ALCtoAC, ACCtoALAR +using TensorKit +using LinearAlgebra + +begin "test utility" + ds = [ℂ^2] + Ds = [ℂ^3] + χs = [ℂ^4] +end + +@testset "initialize A C for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + C = initial_C(A) + @test size(A) == (Ni, Nj) + @test size(C) == (Ni, Nj) + @test all(i -> space(i) == (χ * D * D' ← χ), A) + @test all(i -> space(i) == (χ ← χ), C) +end + +@testset "getL!, getAL and getLsped for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + A = initial_A(ipeps, χ) + C = initial_C(A) + C = ρmap(C, A) + + L = getL!(A, C) + + @test all(i -> all(j -> real(j) > 0, diag(i).values[]), L) + @test all(i -> all(j -> imag(j) ≈ 0, diag(i).values[]), L) + @test all(i -> space(i) == (χ ← χ), L) + + AL, Le, λ = getAL(A, L) + @test all(i -> space(i) == (χ * D * D' ← χ), AL) + @test all(i -> space(i) == (χ ← χ), Le) + @test all(map((λ, AL, Le, A, L) -> λ * AL * Le ≈ transpose(L*transpose(A, ((1,),(4,3,2))), ((1,4,3),(2,))), λ, AL, Le, A, L)) + + L = getLsped(Le, A, AL) + @test all(i -> space(i) == (χ ← χ), L) + @test all(i -> all(j -> real(j) > 0, diag(i).values[]), L) + @test all(i -> all(j -> imag(j) ≈ 0, diag(i).values[]), L) +end + +@testset "canonical form for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + A = initial_A(ipeps, χ) + + AL, L, λ = left_canonical(A) + @test all(i -> space(i) == (χ * D * D' ← χ), AL) + @test all(i -> space(i) == (χ ← χ), L) + @test all(AL -> (AL' * AL ≈ isomorphism(χ, χ)), AL) + @test all(map((A, AL, L, λ) -> λ * AL * L ≈ _to_front(L * _to_tail(A)), A, AL, L, λ)) + + R, AR, λ = right_canonical(A) + @test all(i -> space(i) == (χ * D * D' ← χ), AR) + @test all(i -> space(i) == (χ ← χ), R) + @test all(AR -> (_to_tail(AR) * _to_tail(AR)' ≈ isomorphism(χ, χ)), AR) + @test all(map((A, R, AR, λ) -> _to_front(λ * R * _to_tail(AR)) ≈ A * R, A, R, AR, λ)) +end + +@testset "initialize FL FR for unitcell $Ni x $Nj" for Ni in 1:2, Nj in 1:2, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + FL = initial_FL(AL, ipeps) + @test all(i -> space(i) == (χ * D' * D ← χ), FL) + + FR = initial_FR(AR, ipeps) + @test all(i -> space(i) == (χ * D * D' ← χ), FR) +end + +@testset "leftenv and rightenv for unitcell $Ni x $Nj" for Ni in 1:3, Nj in 1:3, (d, D, χ) in zip(ds, Ds, χs), ifobs in [true, false] + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, FL = leftenv(AL, adjoint.(AL), ipeps; ifobs) + λR, FR = rightenv(AR, adjoint.(AR), ipeps; ifobs) + + @test all(i -> space(i) == (χ * D' * D ← χ), FL) + @test all(i -> space(i) == (χ * D * D' ← χ), FR) + + for i in 1:Ni + ir = ifobs ? Ni + 1 - i : mod1(i + 1, Ni) + @test λL[i] * FL[i,:] ≈ FLmap(FL[i,:], AL[i,:], adjoint.(AL)[ir,:], ipeps[i,:], adjoint.(ipeps[i,:])) rtol = 1e-12 + @test λR[i] * FR[i,:] ≈ FRmap(FR[i,:], AR[i,:], adjoint.(AR)[ir,:], ipeps[i,:], adjoint.(ipeps[i,:])) rtol = 1e-12 + end +end + +@testset "ACenv and Cenv for unitcell $Ni x $Nj" for Ni in 1:3, Nj in 1:3, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, FL = leftenv(AL, adjoint.(AL), ipeps) + λR, FR = rightenv(AR, adjoint.(AR), ipeps) + + C = LRtoC(L, R) + AC = ALCtoAC(AL, C) + + λAC, AC = ACenv(AC, FL, FR, ipeps) + λC, C = Cenv( C, FL, FR) + @test all(i -> space(i) == (χ * D * D' ← χ), AC) + @test all(i -> space(i) == (χ ← χ), C) + + for j in 1:Nj + jr = mod1(j + 1, Nj) + @test λAC[j] * AC[:,j] ≈ ACmap(AC[:,j], FL[:,j], FR[:,j], ipeps[:,j], adjoint.(ipeps[:,j])) rtol = 1e-12 + @test λC[j] * C[:,j] ≈ Cmap( C[:,j], FL[:,jr], FR[:,j]) rtol = 1e-10 + end +end + +@testset "ACCtoALAR for unitcell $Ni x $Nj" for Ni in 1:3, Nj in 1:3, (d, D, χ) in zip(ds, Ds, χs) + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, FL = leftenv(AL, adjoint.(AL), ipeps) + λR, FR = rightenv(AR, adjoint.(AR), ipeps) + + C = LRtoC(L, R) + AC = ALCtoAC(AL, C) + + λAC, AC = ACenv(AC, FL, FR, ipeps) + λC, C = Cenv( C, FL, FR) + + AL, AR, errL, errR = ACCtoALAR(AC, C) + @test all(i -> space(i) == (χ * D * D' ← χ), AL) + @test all(i -> space(i) == (χ * D * D' ← χ), AR) + @test all(AL -> (AL' * AL ≈ isomorphism(χ, χ)), AL) + @test all(AR -> (_to_tail(AR) * _to_tail(AR)' ≈ isomorphism(χ, χ)), AR) + @test errL isa Real + @test errR isa Real +end + +@testset "leftCenv and rightCenv for unitcell $Ni x $Nj" for Ni in 1:3, Nj in 1:3, (d, D, χ) in zip(ds, Ds, χs), ifobs in [true, false] + Random.seed!(42) + ipeps = InfinitePEPS(d, D; unitcell=(Ni, Nj)) + + A = initial_A(ipeps, χ) + AL, L, λ = left_canonical(A) + R, AR, λ = right_canonical(A) + + λL, L = leftCenv(AL, adjoint.(AL); ifobs) + λR, R = rightCenv(AR, adjoint.(AR); ifobs) + @test all(i -> space(i) == (χ ← χ), R) + @test all(i -> space(i) == (χ ← χ), L) + + for i in 1:Ni + ir = ifobs ? mod1(Ni + 2 - i, Ni) : i + @test λL[i] * L[i,:] ≈ Lmap(L[i,:], AL[i,:], adjoint.(AL)[ir,:]) rtol = 1e-12 + @test λR[i] * R[i,:] ≈ Rmap(R[i,:], AR[i,:], adjoint.(AR)[ir,:]) rtol = 1e-12 + end +end \ No newline at end of file