diff --git a/Project.toml b/Project.toml index 1b4246ddb..25f73cd16 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MPSKit" uuid = "bb1c41ca-d63c-52ed-829e-0820dda26502" authors = "Lukas Devos, Maarten Van Damme and contributors" -version = "0.12.6" +version = "0.13.0-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" @@ -25,7 +25,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] Accessors = "0.1" Aqua = "0.8.9" -BlockTensorKit = "0.1.4" +BlockTensorKit = "0.1.6" Combinatorics = "1" Compat = "3.47, 4.10" DocStringExtensions = "0.9.3" diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 095da99b1..83647881f 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -60,6 +60,7 @@ using Compat: @compat using TensorKit using TensorKit: BraidingTensor using BlockTensorKit +using BlockTensorKit: TensorMapSumSpace using TensorOperations using KrylovKit using KrylovKit: KrylovAlgorithm @@ -75,7 +76,7 @@ using DocStringExtensions using LinearAlgebra: diag, Diagonal using LinearAlgebra: LinearAlgebra using Random -using Base: @kwdef +using Base: @kwdef, @propagate_inbounds using LoggingExtras using OhMyThreads @@ -112,6 +113,7 @@ include("states/ortho.jl") include("operators/abstractmpo.jl") include("operators/mpo.jl") +include("operators/jordanmpotensor.jl") include("operators/mpohamiltonian.jl") # the mpohamiltonian objects include("operators/multilinempo.jl") include("operators/projection.jl") @@ -131,7 +133,10 @@ include("environments/multiple_envs.jl") include("environments/lazylincocache.jl") include("algorithms/fixedpoint.jl") -include("algorithms/derivatives.jl") +include("algorithms/derivatives/derivatives.jl") +include("algorithms/derivatives/mpo_derivatives.jl") +include("algorithms/derivatives/hamiltonian_derivatives.jl") +include("algorithms/derivatives/projection_derivatives.jl") include("algorithms/expval.jl") include("algorithms/toolbox.jl") include("algorithms/grassmann.jl") diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 472215629..3d9ad6fb2 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -72,7 +72,7 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili for col in 1:size(ψ, 2) for row in 1:size(ψ, 1) AC2′ = ac2_proj(row, col, ψ, toapprox, envs) - al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AL[row + 1, col] = al @@ -90,10 +90,8 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili for row in 1:size(ψ, 1) # TODO: also write this as ac2_proj? AC2 = ϕ.AL[row, col] * _transpose_tail(ϕ.AC[row, col + 1]) - AC2′ = ∂AC2(AC2, O[row, col], O[row, col + 1], - leftenv(envs[row], col, ψ[row]), - rightenv(envs[row], col, ψ[row])) - al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + AC2′ = ∂∂AC2(row, col, ψ, O, envs) * AC2 + al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AL[row + 1, col] = al diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index 4a917afa8..1c273472c 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/src/algorithms/changebonds/optimalexpand.jl @@ -9,9 +9,12 @@ dominant contributions of a two-site updated MPS tensor, orthogonal to the origi $(TYPEDFIELDS) """ -@kwdef struct OptimalExpand <: Algorithm +@kwdef struct OptimalExpand{S} <: Algorithm + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + "algorithm used for truncating the expanded space" - trscheme::TruncationScheme = truncdim(1) + trscheme::TruncationScheme end function changebonds(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, alg::OptimalExpand, @@ -28,7 +31,7 @@ function changebonds(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, alg::OptimalExp VL = leftnull(ψ.AL[i]) VR = rightnull!(_transpose_tail(ψ.AR[i + 1])) intermediate = adjoint(VL) * AC2 * adjoint(VR) - U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=SVD()) + U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=alg.alg_svd) AL′[i] = VL * U AR′[i + 1] = V * VR @@ -54,7 +57,7 @@ function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs=environments( VL = leftnull(ψ.AL[i, j]) VR = rightnull!(_transpose_tail(ψ.AR[i, j + 1])) intermediate = adjoint(VL) * AC2 * adjoint(VR) - U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=SVD()) + U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=alg.alg_svd) AL′[i, j] = VL * U AR′[i, j + 1] = V * VR @@ -85,7 +88,7 @@ function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs=environ #Use this nullspaces and SVD decomposition to determine the optimal expansion space intermediate = adjoint(NL) * AC2 * adjoint(NR) - _, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=SVD()) + _, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=alg.alg_svd) ar_re = V * NR ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(V, 1))) diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index c2c5f7cd9..d63a9f5bc 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/src/algorithms/changebonds/randexpand.jl @@ -9,9 +9,12 @@ two-site MPS tensor, which is made orthogonal to the existing state. $(TYPEDFIELDS) """ -@kwdef struct RandExpand <: Algorithm +@kwdef struct RandExpand{S} <: Algorithm + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + "algorithm used for [truncation](@extref TensorKit.tsvd] the expanded space" - trscheme::TruncationScheme = truncdim(1) + trscheme::TruncationScheme end function changebonds(ψ::InfiniteMPS, alg::RandExpand) @@ -26,7 +29,7 @@ function changebonds(ψ::InfiniteMPS, alg::RandExpand) VL = leftnull(ψ.AL[i]) VR = rightnull!(_transpose_tail(ψ.AR[i + 1])) intermediate = adjoint(VL) * AC2 * adjoint(VR) - U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=SVD()) + U, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=alg.alg_svd) AL′[i] = VL * U AR′[i + 1] = V * VR @@ -50,7 +53,7 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::RandExpand) #Use this nullspaces and SVD decomposition to determine the optimal expansion space intermediate = adjoint(NL) * AC2 * adjoint(NR) - _, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=SVD()) + _, _, V, = tsvd!(intermediate; trunc=alg.trscheme, alg=alg.alg_svd) ar_re = V * NR ar_le = zerovector!(similar(ar_re, codomain(ψ.AC[i]) ← space(V, 1))) diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 330b7ef06..85497a2ea 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -7,9 +7,12 @@ An algorithm that uses truncated SVD to change the bond dimension of a ψ. $(TYPEDFIELDS) """ -@kwdef struct SvdCut <: Algorithm +@kwdef struct SvdCut{S} <: Algorithm + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + "algorithm used for [truncation][@extref TensorKit.tsvd] of the gauge tensors" - trscheme::TruncationScheme = notrunc() + trscheme::TruncationScheme end function changebonds(ψ::AbstractFiniteMPS, alg::SvdCut; kwargs...) @@ -17,7 +20,7 @@ function changebonds(ψ::AbstractFiniteMPS, alg::SvdCut; kwargs...) end function changebonds!(ψ::AbstractFiniteMPS, alg::SvdCut; normalize::Bool=true) for i in (length(ψ) - 1):-1:1 - U, S, V, = tsvd(ψ.C[i]; trunc=alg.trscheme, alg=TensorKit.SVD()) + U, S, V, = tsvd(ψ.C[i]; trunc=alg.trscheme, alg=alg.alg_svd) AL′ = ψ.AL[i] * U ψ.AC[i] = (AL′, complex(S)) AR′ = _transpose_front(V * _transpose_tail(ψ.AR[i + 1])) @@ -41,7 +44,7 @@ function changebonds!(mpo::FiniteMPO, alg::SvdCut) O_left = transpose(mpo[1], ((3, 1, 2), (4,))) local O_right for i in 2:length(mpo) - U, S, V, = tsvd!(O_left; trunc=alg.trscheme, alg=TensorKit.SVD()) + U, S, V, = tsvd!(O_left; trunc=alg.trscheme, alg=alg.alg_svd) @inbounds mpo[i - 1] = transpose(U, ((2, 3), (1, 4))) if i < length(mpo) @plansor O_left[-3 -1 -2; -4] := S[-1; 1] * V[1; 2] * mpo[i][2 -2; -3 -4] @@ -52,7 +55,7 @@ function changebonds!(mpo::FiniteMPO, alg::SvdCut) # right to left for i in (length(mpo) - 1):-1:1 - U, S, V, = tsvd!(O_right; trunc=alg.trscheme, alg=TensorKit.SVD()) + U, S, V, = tsvd!(O_right; trunc=alg.trscheme, alg=alg.alg_svd) @inbounds mpo[i + 1] = transpose(V, ((1, 4), (2, 3))) if i > 1 @plansor O_right[-1; -3 -4 -2] := mpo[i][-1 -2; -3 2] * U[2; 1] * S[1; -4] @@ -81,7 +84,7 @@ function changebonds(ψ::InfiniteMPS, alg::SvdCut) ncr = ψ.C[1] for i in 1:length(ψ) - U, ncr, = tsvd(ψ.C[i]; trunc=alg.trscheme, alg=TensorKit.SVD()) + U, ncr, = tsvd(ψ.C[i]; trunc=alg.trscheme, alg=alg.alg_svd) copied[i] = copied[i] * U copied[i + 1] = _transpose_front(U' * _transpose_tail(copied[i + 1])) end diff --git a/src/algorithms/changebonds/vumpssvd.jl b/src/algorithms/changebonds/vumpssvd.jl index 83f45837d..d9e2e4da7 100644 --- a/src/algorithms/changebonds/vumpssvd.jl +++ b/src/algorithms/changebonds/vumpssvd.jl @@ -8,12 +8,17 @@ An algorithm that uses a two-site update step to change the bond dimension of a $(TYPEDFIELDS) """ @kwdef struct VUMPSSvdCut <: Algorithm - "tolerance for gauging algorithm" - tol_gauge = Defaults.tolgauge - "tolerance for the eigenvalue solver" - tol_eigenval = Defaults.tol + "algorithm used for gauging the `InfiniteMPS`" + alg_gauge = Defaults.alg_gauge(; dynamic_tols=false) + + "algorithm used for the eigenvalue solvers" + alg_eigsolve = Defaults.alg_eigsolve(; dynamic_tols=false) + + "algorithm used for the singular value decomposition" + alg_svd = Defaults.alg_svd() + "algorithm used for [truncation][@extref TensorKit.tsvd] of the two-site update" - trscheme::TruncationScheme = notrunc() + trscheme::TruncationScheme end function changebonds_1(state::InfiniteMPS, H, alg::VUMPSSvdCut, @@ -31,11 +36,12 @@ function changebonds_1(state::InfiniteMPS, H, alg::VUMPSSvdCut, # collapse back to 1 site if D2 != D1 - (nstate, nenvs) = changebonds(nstate, nH, - SvdCut(; trscheme=truncspace(infimum(D1, D2))), nenvs) + cut_alg = SvdCut(; alg.alg_svd, trscheme=truncspace(infimum(D1, D2))) + nstate, nenvs = changebonds(nstate, nH, cut_alg, nenvs) end - collapsed = InfiniteMPS([nstate.AL[1]], nstate.C[1]; tol=alg.tol_gauge) + collapsed = InfiniteMPS([nstate.AL[1]], nstate.C[1]; alg.alg_gauge.tol, + alg.alg_gauge.maxiter) recalculate!(envs, collapsed, H, collapsed) return collapsed, envs @@ -47,17 +53,13 @@ function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs=environment @plansor AC2[-1 -2; -3 -4] := state.AC[loc][-1 -2; 1] * state.AR[loc + 1][1 -4; -3] h_ac2 = ∂∂AC2(loc, state, H, envs) - (vals, vecs, _) = eigsolve(h_ac2, AC2, 1, :SR; tol=alg.tol_eigenval, - ishermitian=false) - nAC2 = vecs[1] + _, nAC2 = fixedpoint(h_ac2, AC2, :SR, alg.alg_eigsolve) h_c = ∂∂C(loc + 1, state, H, envs) - (vals, vecs, _) = eigsolve(h_c, state.C[loc + 1], 1, :SR; tol=alg.tol_eigenval, - ishermitian=false) - nC2 = vecs[1] + _, nC2 = fixedpoint(h_c, state.C[loc + 1], :SR, alg.alg_eigsolve) #svd ac2, get new AL1 and S,V ---> AC - (AL1, S, V, eps) = tsvd(nAC2; trunc=alg.trscheme, alg=TensorKit.SVD()) + AL1, S, V, eps = tsvd!(nAC2; trunc=alg.trscheme, alg=alg.alg_svd) @plansor AC[-1 -2; -3] := S[-1; 1] * V[1; -3 -2] meps = max(eps, meps) @@ -72,16 +74,13 @@ function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs=environment copied = copy(state.AL) copied[loc] = AL1 copied[loc + 1] = AL2 - state = InfiniteMPS(copied; tol=alg.tol_gauge) + state = InfiniteMPS(copied; alg.alg_gauge.tol, alg.alg_gauge.maxiter) recalculate!(envs, state, H, state) end return state, envs end function changebonds(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs=environments(state, H)) - if (length(state) == 1) - return changebonds_1(state, H, alg, envs) - else - return changebonds_n(state, H, alg, envs) - end + return length(state) == 1 ? changebonds_1(state, H, alg, envs) : + changebonds_n(state, H, alg, envs) end diff --git a/src/algorithms/derivatives.jl b/src/algorithms/derivatives.jl deleted file mode 100644 index f75000885..000000000 --- a/src/algorithms/derivatives.jl +++ /dev/null @@ -1,264 +0,0 @@ -# Given a state and it's environments, we can act on it - -""" - Draft operators -""" -struct MPO_∂∂C{L,R} - leftenv::L - rightenv::R -end - -struct MPO_∂∂AC{O,L,R} - o::O - leftenv::L - rightenv::R -end - -struct MPO_∂∂AC2{O,L,R} - o1::O - o2::O - leftenv::L - rightenv::R -end - -const DerivativeOperator = Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2} - -Base.:*(h::Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}, v) = h(v); - -(h::MPO_∂∂C)(x) = ∂C(x, h.leftenv, h.rightenv); -(h::MPO_∂∂AC)(x) = ∂AC(x, h.o, h.leftenv, h.rightenv); -(h::MPO_∂∂AC2)(x) = ∂AC2(x, h.o1, h.o2, h.leftenv, h.rightenv); -(h::DerivativeOperator)(v, ::Number) = h(v) - -# draft operator constructors -function ∂∂C(pos::Int, mps, operator, envs) - return MPO_∂∂C(leftenv(envs, pos + 1, mps), rightenv(envs, pos, mps)) -end -function ∂∂C(row::Int, col::Int, mps, operator::MultilineMPO, envs) - return ∂∂C(col, mps[row], operator[row], envs[row]) -end - -function ∂∂AC(pos::Int, mps, operator, envs) - return MPO_∂∂AC(operator[pos], leftenv(envs, pos, mps), rightenv(envs, pos, mps)) -end -function ∂∂AC(row::Int, col::Int, mps, operator::MultilineMPO, envs) - return ∂∂AC(col, mps[row], operator[row], envs[row]) -end -function ∂∂AC(col::Int, mps, operator::MultilineMPO, envs) - return MPO_∂∂AC(operator[:, col], leftenv(envs, col, mps), rightenv(envs, col, mps)) -end - -function ∂∂AC2(pos::Int, mps, operator, envs) - return MPO_∂∂AC2(operator[pos], operator[pos + 1], leftenv(envs, pos, mps), - rightenv(envs, pos + 1, mps)) -end -function ∂∂AC2(col::Int, mps, operator::MultilineMPO, envs) - return MPO_∂∂AC2(operator[:, col], operator[:, col + 1], leftenv(envs, col, mps), - rightenv(envs, col + 1, mps)) -end -function ∂∂AC2(row::Int, col::Int, mps, operator::MultilineMPO, envs::MultilineEnvironments) - return ∂∂AC2(col, mps[row], operator[row], envs[row]) -end - -# allow calling them with CartesianIndices -∂∂C(pos::CartesianIndex, mps, operator, envs) = ∂∂C(Tuple(pos)..., mps, operator, envs) -∂∂AC(pos::CartesianIndex, mps, operator, envs) = ∂∂AC(Tuple(pos)..., mps, operator, envs) -∂∂AC2(pos::CartesianIndex, mps, operator, envs) = ∂∂AC2(Tuple(pos)..., mps, operator, envs) - -function ∂AC(x::MPSTensor{S}, operator::MPOTensor{S}, leftenv::MPSTensor{S}, - rightenv::MPSTensor{S})::typeof(x) where {S} - @plansor y[-1 -2; -3] := leftenv[-1 5; 4] * x[4 2; 1] * operator[5 -2; 2 3] * - rightenv[1 3; -3] - return y isa BlockTensorMap ? only(y) : y -end -function ∂AC(x::MPSTensor{S}, operator::Number, leftenv::MPSTensor{S}, - rightenv::MPSTensor{S})::typeof(x) where {S} - @plansor y[-1 -2; -3] := operator * (leftenv[-1 5; 4] * x[4 6; 1] * τ[6 5; 7 -2] * - rightenv[1 7; -3]) -end -function ∂AC(x::GenericMPSTensor{S,3}, operator::MPOTensor{S}, leftenv::MPSTensor{S}, - rightenv::MPSTensor{S})::typeof(x) where {S} - @plansor y[-1 -2 -3; -4] ≔ leftenv[-1 7; 6] * x[6 4 2; 1] * operator[7 -2; 4 5] * - τ[5 -3; 2 3] * rightenv[1 3; -4] -end - -# mpo multiline -function ∂AC(x::AbstractVector, opp, leftenv, rightenv) - return circshift(map(∂AC, x, opp, leftenv, rightenv), 1) -end - -function ∂AC(x::MPSTensor, ::Nothing, leftenv, rightenv) - return _transpose_front(leftenv * _transpose_tail(x * rightenv)) -end - -function ∂AC2(x::MPOTensor, operator1::MPOTensor, operator2::MPOTensor, leftenv, rightenv) - @plansor toret[-1 -2; -3 -4] := leftenv[-1 7; 6] * x[6 5; 1 3] * - operator1[7 -2; 5 4] * - operator2[4 -4; 3 2] * - rightenv[1 2; -3] - return toret isa BlockTensorMap ? only(toret) : toret -end -function ∂AC2(x::MPOTensor, ::Nothing, ::Nothing, leftenv, rightenv) - @plansor y[-1 -2; -3 -4] := x[1 -2; 2 -4] * leftenv[-1; 1] * rightenv[2; -3] -end -function ∂AC2(x::AbstractTensorMap{<:Any,<:Any,3,3}, operator1::MPOTensor, - operator2::MPOTensor, leftenv::MPSTensor, rightenv::MPSTensor)::typeof(x) - @plansor y[-1 -2 -3; -4 -5 -6] ≔ leftenv[-1 11; 10] * x[10 8 6; 1 2 4] * - rightenv[1 3; -4] * - operator1[11 -2; 8 9] * τ[9 -3; 6 7] * - operator2[7 -6; 4 5] * τ[5 -5; 2 3] -end - -function ∂AC2(x::AbstractVector, opp1, opp2, leftenv, rightenv) - return circshift(map(∂AC2, x, opp1, opp2, leftenv, rightenv), 1) -end - -function ∂C(x::MPSBondTensor, leftenv::MPSTensor, rightenv::MPSTensor) - @plansor y[-1; -2] := leftenv[-1 3; 1] * x[1; 2] * rightenv[2 3; -2] - return y isa BlockTensorMap ? only(y) : y -end - -function ∂C(x::MPSBondTensor, leftenv::MPSBondTensor, rightenv::MPSBondTensor) - @plansor toret[-1; -2] := leftenv[-1; 1] * x[1; 2] * rightenv[2; -2] -end - -function ∂C(x::AbstractVector, leftenv, rightenv) - return circshift(map(t -> ∂C(t...), zip(x, leftenv, rightenv)), 1) -end - -# downproject for approximate -function c_proj(pos::Int, ψ, (operator, ϕ)::Tuple, envs) - return ∂C(ϕ.C[pos], leftenv(envs, pos + 1, ψ), rightenv(envs, pos, ψ)) -end -function c_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) - return ∂C(ϕ.C[pos], leftenv(envs, pos + 1, ψ), rightenv(envs, pos, ψ)) -end -function c_proj(pos::Int, ψ, Oϕs::LazySum, envs) - return sum(zip(Oϕs.ops, envs.envs)) do x - return c_proj(pos, ψ, x...) - end -end -function c_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) - return c_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) -end - -function ac_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs) - return ∂AC(ϕ.AC[pos], O[pos], leftenv(envs, pos, ψ), rightenv(envs, pos, ψ)) -end -function ac_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) - return ∂AC(ϕ.AC[pos], nothing, leftenv(envs, pos, ψ), rightenv(envs, pos, ψ)) -end -function ac_proj(pos::Int, ψ, Oϕs::LazySum, envs) - return sum(zip(Oϕs.ops, envs.envs)) do x - return ac_proj(pos, ψ, x...) - end -end -function ac_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) - return ac_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) -end - -function ac2_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs) - AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1]) - return ∂AC2(AC2, O[pos], O[pos + 1], leftenv(envs, pos, ψ), rightenv(envs, pos + 1, ψ)) -end -function ac2_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) - AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1]) - return ∂AC2(AC2, nothing, nothing, leftenv(envs, pos, ψ), rightenv(envs, pos + 1, ψ)) -end -function ac2_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) - return ac2_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) -end - -function ∂∂C(pos::Int, mps, operator::LinearCombination, cache) - return LinearCombination(broadcast((h, e) -> ∂∂C(pos, mps, h, e), operator.opps, - cache.envs), - operator.coeffs) -end - -function ∂∂AC(pos::Int, mps, operator::LinearCombination, cache) - return LinearCombination(broadcast((h, e) -> ∂∂AC(pos, mps, h, e), operator.opps, - cache.envs), operator.coeffs) -end - -function ∂∂AC2(pos::Int, mps, operator::LinearCombination, cache) - return LinearCombination(broadcast((h, e) -> ∂∂AC2(pos, mps, h, e), operator.opps, - cache.envs), operator.coeffs) -end - -struct AC_EffProj{A,L} - a1::A - le::L - re::L -end -struct AC2_EffProj{A,L} - a1::A - a2::A - le::L - re::L -end -Base.:*(h::Union{AC_EffProj,AC2_EffProj}, v) = h(v); - -function (h::AC_EffProj)(x::MPSTensor) - @plansor v[-1; -2 -3 -4] := h.le[4; -1 -2 5] * h.a1[5 2; 1] * h.re[1; -3 -4 3] * - conj(x[4 2; 3]) - @plansor y[-1 -2; -3] := conj(v[1; 2 5 6]) * h.le[-1; 1 2 4] * h.a1[4 -2; 3] * - h.re[3; 5 6 -3] -end -function (h::AC2_EffProj)(x::MPOTensor) - @plansor v[-1; -2 -3 -4] := h.le[6; -1 -2 7] * h.a1[7 4; 5] * h.a2[5 2; 1] * - h.re[1; -3 -4 3] * conj(x[6 4; 3 2]) - @plansor y[-1 -2; -3 -4] := conj(v[2; 3 5 6]) * h.le[-1; 2 3 4] * h.a1[4 -2; 7] * - h.a2[7 -4; 1] * h.re[1; 5 6 -3] -end - -function ∂∂AC(pos::Int, state, operator::ProjectionOperator, env) - return AC_EffProj(operator.ket.AC[pos], leftenv(env, pos, state), - rightenv(env, pos, state)) -end -function ∂∂AC2(pos::Int, state, operator::ProjectionOperator, env) - return AC2_EffProj(operator.ket.AC[pos], operator.ket.AR[pos + 1], - leftenv(env, pos, state), - rightenv(env, pos + 1, state)) -end - -# time dependent derivate operators -(h::UntimedOperator{<:DerivativeOperator})(y, args...) = h.f * h.op(y) -(h::TimedOperator{<:DerivativeOperator})(y, t::Number) = h.f(t) * h.op(y) -function (x::LazySum{<:Union{MultipliedOperator{D},D} where {D<:DerivativeOperator}})(y, - t::Number) - return sum(O -> O(y, t), x) -end -function (x::LazySum{<:Union{MultipliedOperator{D},D} where {D<:DerivativeOperator}})(y) - return sum(O -> O(y), x) -end -function Base.:*(h::LazySum{<:Union{D,MultipliedOperator{D}} where {D<:DerivativeOperator}}, - v) - return h(v) -end - -function ∂∂C(pos::Int, mps, operator::MultipliedOperator, cache) - return MultipliedOperator(∂∂C(pos::Int, mps, operator.op, cache), operator.f) -end - -function ∂∂AC(pos::Int, mps, operator::MultipliedOperator, cache) - return MultipliedOperator(∂∂AC(pos::Int, mps, operator.op, cache), operator.f) -end - -function ∂∂AC2(pos::Int, mps, operator::MultipliedOperator, cache) - return MultipliedOperator(∂∂AC2(pos::Int, mps, operator.op, cache), operator.f) -end - -function ∂∂C(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) - suboperators = map((op, openv) -> ∂∂C(pos, mps, op, openv), operator.ops, cache.envs) - return LazySum{Union{MPO_∂∂C,MultipliedOperator{<:MPO_∂∂C}}}(suboperators) -end - -function ∂∂AC(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) - suboperators = map((op, openv) -> ∂∂AC(pos, mps, op, openv), operator.ops, cache.envs) - return LazySum{Union{MPO_∂∂AC,MultipliedOperator{<:MPO_∂∂AC}}}(suboperators) -end - -function ∂∂AC2(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) - suboperators = map((op, openv) -> ∂∂AC2(pos, mps, op, openv), operator.ops, cache.envs) - return LazySum{Union{MPO_∂∂AC2,MultipliedOperator{<:MPO_∂∂AC2}}}(suboperators) -end diff --git a/src/algorithms/derivatives/derivatives.jl b/src/algorithms/derivatives/derivatives.jl new file mode 100644 index 000000000..9dc7dfe9a --- /dev/null +++ b/src/algorithms/derivatives/derivatives.jl @@ -0,0 +1,248 @@ +# Given a state and it's environments, we can act on it +""" + DerivativeOperator + +Abstract supertype for derivative operators acting on MPS. These operators are used to represent +the effective local operators obtained from taking the partial derivative of an MPS-MPO-MPS sandwich. +""" +abstract type DerivativeOperator end + +Base.:*(h::DerivativeOperator, v) = h(v) +(h::DerivativeOperator)(v, ::Number) = h(v) + +@doc """ + ∂∂C(site::Int, mps, operator, envs) + ∂∂C(row::Int, col::Int, mps, operator, envs) + +Effective zero-site local operator acting at `site`. + +``` + ┌─ ─┐ + │ │ +┌┴┐ ┌┴┐ +│ ├───┤ │ +└┬┘ └┬┘ + │ │ + └─ ─┘ +``` + +See also [`∂C`](@ref). +""" ∂∂C + +∂∂C(pos::CartesianIndex, mps, operator, envs) = ∂∂C(Tuple(pos)..., mps, operator, envs) +function ∂∂C(row::Int, col::Int, mps::MultilineMPS, operator::MultilineMPO, envs) + return ∂∂C(col, mps[row], operator[row], envs[row]) +end +function ∂∂C(site::Int, mps, operator::MultilineMPO, envs) + return Multiline([∂∂C(row, site, mps, operator, envs) for row in 1:size(operator, 1)]) +end +function ∂∂C(pos::Int, mps, operator::MultipliedOperator, cache) + return MultipliedOperator(∂∂C(pos, mps, operator.op, cache), operator.f) +end +function ∂∂C(pos::Int, mps, operator::LinearCombination, cache) + return LinearCombination(broadcast((h, e) -> ∂∂C(pos, mps, h, e), operator.opps, + cache.envs), + operator.coeffs) +end +function ∂∂C(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) + suboperators = map((op, openv) -> ∂∂C(pos, mps, op, openv), operator.ops, cache.envs) + return LazySum{Union{MPO_∂∂C,MultipliedOperator{<:MPO_∂∂C}}}(suboperators) +end + +""" + ∂C(x, leftenv, rightenv) + +Application of the effective zero-site local operator on a local tensor `x`. + +``` + ┌─┐ + ┌─┤ ├─┐ + │ └─┘ │ +┌┴┐ ┌┴┐ +│ ├───┤ │ +└┬┘ └┬┘ + │ │ + └─ ─┘ +``` + +See also [`∂∂C`](@ref). +""" +∂C(x, leftenv, rightenv) = MPO_∂∂C(leftenv, rightenv)(x) + +@doc """ + ∂∂AC(site::Int, mps, operator, envs) + ∂∂AC(row::Int, col::Int, mps, operator, envs) + +Effective one-site local operator acting at `site`. + +``` + ┌─── ───┐ + │ │ │ +┌┴┐ ┌─┴─┐ ┌┴┐ +│ ├─┤ ├─┤ │ +└┬┘ └─┬─┘ └┬┘ + │ │ │ + └─── ───┘ +``` + +See also [`∂AC`](@ref). +""" ∂∂AC +∂∂AC(pos::CartesianIndex, mps, operator, envs) = ∂∂AC(Tuple(pos)..., mps, operator, envs) +function ∂∂AC(row::Int, col::Int, mps::MultilineMPS, operator::MultilineMPO, envs) + return ∂∂AC(col, mps[row], operator[row], envs[row]) +end +function ∂∂AC(site::Int, mps, operator::MultilineMPO, envs) + return Multiline([∂∂AC(row, site, mps, operator, envs) for row in 1:size(operator, 1)]) +end +function ∂∂AC(pos::Int, mps, operator::MultipliedOperator, cache) + return MultipliedOperator(∂∂AC(pos, mps, operator.op, cache), operator.f) +end +function ∂∂AC(pos::Int, mps, operator::LinearCombination, cache) + return LinearCombination(broadcast((h, e) -> ∂∂AC(pos, mps, h, e), operator.opps, + cache.envs), operator.coeffs) +end +function ∂∂AC(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) + suboperators = map((op, openv) -> ∂∂AC(pos, mps, op, openv), operator.ops, cache.envs) + elT = Union{D,MultipliedOperator{D}} where {D<:Union{MPO_∂∂AC,JordanMPO_∂∂AC}} + return LazySum{elT}(suboperators) +end + +""" + ∂AC(x, operator, leftenv, rightenv) + +Application of the effective one-site local operator on a local tensor `x`. + +``` + ┌───┐ + ┌──┤ ├──┐ + │ └─┬─┘ │ +┌┴┐ ┌─┴─┐ ┌┴┐ +│ ├─┤ ├─┤ │ +└┬┘ └─┬─┘ └┬┘ + │ │ │ + └── ──┘ +``` + +See also [`∂∂AC`](@ref). +""" +∂AC(x, operator, leftenv, rightenv) = MPO_∂∂AC(leftenv, operator, rightenv)(x) + +@doc """ + ∂∂AC2(site::Int, mps, operator, envs) + ∂∂AC2(row::Int, col::Int, mps, operator, envs) + +Effective two-site local operator acting at `site`. + +``` + ┌── ──┐ + │ │ │ │ +┌┴┐┌─┴─┐┌─┴─┐┌┴┐ +│ ├┤ ├┤ ├┤ │ +└┬┘└─┬─┘└─┬─┘└┬┘ + │ │ │ │ + └── ──┘ +``` + +See also [`∂AC2`](@ref). +""" ∂∂AC2 + +∂∂AC2(pos::CartesianIndex, mps, operator, envs) = ∂∂AC2(Tuple(pos)..., mps, operator, envs) +function ∂∂AC2(row::Int, col::Int, mps::MultilineMPS, operator::MultilineMPO, envs) + return ∂∂AC2(col, mps[row], operator[row], envs[row]) +end +function ∂∂AC2(site::Int, mps, operator::MultilineMPO, envs) + return Multiline([∂∂AC2(row, site, mps, operator, envs) for row in 1:size(operator, 1)]) +end +function ∂∂AC2(pos::Int, mps, operator::MultipliedOperator, cache) + return MultipliedOperator(∂∂AC2(pos, mps, operator.op, cache), operator.f) +end +function ∂∂AC2(pos::Int, mps, operator::LinearCombination, cache) + return LinearCombination(broadcast((h, e) -> ∂∂AC2(pos, mps, h, e), operator.opps, + cache.envs), operator.coeffs) +end +function ∂∂AC2(pos::Int, mps, operator::LazySum, cache::MultipleEnvironments) + suboperators = map((op, openv) -> ∂∂AC2(pos, mps, op, openv), operator.ops, cache.envs) + elT = Union{D,MultipliedOperator{D}} where {D<:Union{MPO_∂∂AC2,JordanMPO_∂∂AC2}} + return LazySum{elT}(suboperators) +end + +""" + ∂AC2(x, operator1, operator2, leftenv, rightenv) + +Application of the effective two-site local operator on a local tensor `x`. + +``` + ┌──────┐ + ┌──┤ ├──┐ + │ └┬────┬┘ │ +┌┴┐┌─┴─┐┌─┴─┐┌┴┐ +│ ├┤ ├┤ ├┤ │ +└┬┘└─┬─┘└─┬─┘└┬┘ + │ │ │ │ + └── ──┘ +``` + +See also [`∂∂AC2`](@ref). +""" +∂AC2(x, O₁, O₂, leftenv, rightenv) = MPO_∂∂AC2(leftenv, O₁, O₂, rightenv)(x) + +# Projection operators +# -------------------- +c_proj(pos::Int, ψ, (operator, ϕ)::Tuple, envs) = ∂∂C(pos, ψ, operator, envs) * ϕ.C[pos] +c_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) = ∂∂C(pos, ψ, nothing, envs) * ϕ.C[pos] +function c_proj(pos::Int, ψ, Oϕs::LazySum, envs) + return sum(zip(Oϕs.ops, envs.envs)) do x + return c_proj(pos, ψ, x...) + end +end +function c_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) + return c_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) +end + +ac_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs) = ∂∂AC(pos, ψ, O, envs) * ϕ.AC[pos] +ac_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) = ∂∂AC(pos, ψ, nothing, envs) * ϕ.AC[pos] + +function ac_proj(pos::Int, ψ, Oϕs::LazySum, envs) + return sum(zip(Oϕs.ops, envs.envs)) do x + return ac_proj(pos, ψ, x...) + end +end +function ac_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) + return ac_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) +end + +function ac2_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs) + AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1]) + return ∂∂AC2(pos, ψ, O, envs) * AC2 +end +function ac2_proj(pos::Int, ψ, ϕ::AbstractMPS, envs) + AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1]) + return ∂∂AC2(pos, ψ, nothing, envs) * AC2 +end +function ac2_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs) + return ac2_proj(col, ψ[row], (O[row], ϕ[row]), envs[row]) +end + +# Varia +# ----- + +# Multiline +function (H::Multiline{<:DerivativeOperator})(x::AbstractVector) + return [H[row - 1] * x[mod1(row - 1, end)] for row in 1:size(H, 1)] +end +Base.:*(H::Multiline{<:DerivativeOperator}, x::AbstractVector) = H(x) + +# time dependent derivative operators +(h::UntimedOperator{<:DerivativeOperator})(y, args...) = h.f * h.op(y) +(h::TimedOperator{<:DerivativeOperator})(y, t::Number) = h.f(t) * h.op(y) +function (x::LazySum{<:Union{MultipliedOperator{D},D} where {D<:DerivativeOperator}})(y, + t::Number) + return sum(O -> O(y, t), x) +end +function (x::LazySum{<:Union{MultipliedOperator{D},D} where {D<:DerivativeOperator}})(y) + return sum(O -> O(y), x) +end +function Base.:*(h::LazySum{<:Union{D,MultipliedOperator{D}} where {D<:DerivativeOperator}}, + v) + return h(v) +end diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl new file mode 100644 index 000000000..c7dbea837 --- /dev/null +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -0,0 +1,311 @@ +""" + JordanMPO_∂∂AC{O1,O2,O3} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_∂∂AC{O1,O2,O3} <: DerivativeOperator + D::Union{O1,Missing} # onsite + I::Union{O1,Missing} # not started + E::Union{O1,Missing} # finished + C::Union{O2,Missing} # starting + B::Union{O2,Missing} # ending + A::O3 # continuing + function JordanMPO_∂∂AC(onsite, not_started, finished, starting, ending, continuing) + tensor = coalesce(onsite, not_started, finished, starting, ending) + ismissing(tensor) && throw(ArgumentError("unable to determine type")) + S = spacetype(tensor) + M = storagetype(tensor) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + return new{O1,O2,typeof(continuing)}(onsite, not_started, finished, starting, + ending, continuing) + end + function JordanMPO_∂∂AC{O1,O2,O3}(onsite, not_started, finished, starting, ending, + continuing) where {O1,O2,O3} + return new{O1,O2,O3}(onsite, not_started, finished, starting, ending, continuing) + end +end + +""" + JordanMPO_∂∂AC2{O1,O2,O3,O4} + +Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. +In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. +""" +struct JordanMPO_∂∂AC2{O1,O2,O3,O4} <: DerivativeOperator + II::Union{O1,Missing} # not_started + IC::Union{O2,Missing} # starting right + ID::Union{O1,Missing} # onsite right + CB::Union{O2,Missing} # starting left - ending right + CA::Union{O3,Missing} # starting left - continuing right + AB::Union{O3,Missing} # continuing left - ending right + AA::O4 # continuing left - continuing right + BE::Union{O2,Missing} # ending left + DE::Union{O1,Missing} # onsite left + EE::Union{O1,Missing} # finished + function JordanMPO_∂∂AC2(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + tensor = coalesce(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + ismissing(tensor) && throw(ArgumentError("unable to determine type")) + S = spacetype(tensor) + M = storagetype(tensor) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = tensormaptype(S, 3, 3, M) + return new{O1,O2,O3,typeof(AA)}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end + function JordanMPO_∂∂AC2{O1,O2,O3,O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, + EE) where {O1,O2,O3,O4} + return new{O1,O2,O3,O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + end +end + +# Constructors +# ------------ +const _HAM_MPS_TYPES = Union{FiniteMPS{<:MPSTensor},WindowMPS{<:MPSTensor}, + InfiniteMPS{<:MPSTensor}} + +function ∂∂AC(site::Int, mps::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, + envs) + GL = leftenv(envs, site, mps) + GR = rightenv(envs, site, mps) + W = operator[site] + + # starting + if nonzero_length(W.C) > 0 + GR_2 = GR[2:(end - 1)] + @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] + starting = only(starting_) + else + starting = missing + end + + # ending + if nonzero_length(W.B) > 0 + GL_2 = GL[2:(end - 1)] + @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] + ending = only(ending_) + else + ending = missing + end + + # onsite + if nonzero_length(W.D) > 0 + if !ismissing(starting) + @plansor starting[-1 -2; -3 -4] += W.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + onsite = missing + elseif !ismissing(ending) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W.D[-2; -4] + onsite = missing + else + onsite = W.D + end + else + onsite = missing + end + + # not_started + if (!isfinite(operator) || site < length(operator)) && !ismissing(starting) + I = id(storagetype(GR[1]), physicalspace(W)) + @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + not_started = missing + else + not_started = removeunit(GR[1], 2) + end + + # finished + if (!isfinite(operator) || site > 1) && !ismissing(ending) + I = id(storagetype(GL[end]), physicalspace(W)) + @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + finished = missing + else + finished = removeunit(GL[end], 2) + end + + # continuing + A = W.A + continuing = (GL[2:(end - 1)], A, GR[2:(end - 1)]) + + return JordanMPO_∂∂AC(onsite, not_started, finished, starting, ending, continuing) +end + +function ∂∂AC2(site::Int, mps::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, + envs) + GL = leftenv(envs, site, mps) + GR = rightenv(envs, site + 1, mps) + W1 = operator[site] + W2 = operator[site + 1] + + # starting left - continuing right + if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 + @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * + GR[2:(end - 1)][-6 1; -3] + CA = only(CA_) + else + CA = missing + end + + # continuing left - ending right + if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 + @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * W1.A[2 -2; -5 1] * + W2.B[1 -3; -6] + AB = only(AB_) + else + AB = missing + end + + # middle + if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 + if !ismissing(CA) + @plansor CA[-1 -2 -3; -4 -5 -6] += W1.C[-1; -4 1] * W2.B[1 -2; -5] * + removeunit(GR[end], 2)[-6; -3] + CB = missing + elseif !ismissing(AB) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + W1.C[-2; -5 1] * + W2.B[1 -3; -6] + CB = missing + else + @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] + CB = only(CB_) + end + else + CB = missing + end + + # starting right + if nonzero_length(W2.C) > 0 + if !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * + GR[2:(end - 1)][-6 1; -3] + IC = missing + else + @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] + end + else + IC = missing + end + + # ending left + if nonzero_length(W1.B) > 0 + if !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += GL[2:(end - 1)][-1 1; -4] * + (W1.B[1 -2; -5] * I[-3; -6]) + BE = missing + else + @plansor BE[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4] + end + else + BE = missing + end + + # onsite left + if nonzero_length(W1.D) > 0 + if !ismissing(BE) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[1], 2)[-1; -3] * W1.D[-2; -4] + DE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[1], 2)[-1; -4] * + (W1.D[-2; -5] * I[-3; -6]) + DE = missing + # TODO: could also try in CA? + else + DE = only(W1.D) + end + else + DE = missing + end + + # onsite right + if nonzero_length(W2.D) > 0 + if !ismissing(IC) + @plansor IC[-1 -2; -3 -4] += W2.D[-1; -3] * removeunit(GR[end], 2)[-4; -2] + ID = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1)) + @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.D[-2; -5]) * + removeunit(GR[end], 2)[-6; -3] + ID = missing + else + ID = only(W2.D) + end + else + ID = missing + end + + # finished + if !ismissing(IC) + I = id(storagetype(GR[1]), physicalspace(W2)) + @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2] + II = missing + elseif !ismissing(CA) + I = id(storagetype(GR[1]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * removeunit(GR[1], 2)[-6; -3] + II = missing + else + II = transpose(removeunit(GR[1], 2)) + end + + # unstarted + if !ismissing(BE) + I = id(storagetype(GL[end]), physicalspace(W1)) + @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4] + EE = missing + elseif !ismissing(AB) + I = id(storagetype(GL[end]), physicalspace(W1) ⊗ physicalspace(W2)) + @plansor AB[-1 -2 -3; -4 -5 -6] += removeunit(GL[end], 2)[-1; -4] * I[-2 -3; -5 -6] + EE = missing + else + EE = removeunit(GL[end], 2) + end + + # continuing - continuing + # TODO: MPO_∂∂AC2 code reuse + optimization + AA = (GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)]) + + return JordanMPO_∂∂AC2(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) +end + +# Actions +# ------- +function (H::JordanMPO_∂∂AC)(x::MPSTensor) + y = zerovector(x) + + ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] + ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] + ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] + ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] + ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] + + GL, A, GR = H.A + if nonzero_length(A) > 0 + @plansor y[-1 -2; -3] += GL[-1 5; 4] * x[4 2; 1] * A[5 -2; 2 3] * GR[1 3; -3] + end + + return y +end + +function (H::JordanMPO_∂∂AC2)(x::MPOTensor) + y = zerovector(x) + ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] + ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] + ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] + ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] + ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] + ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] + ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] + ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] + ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] + + GL, A1, A2, GR = H.AA + if nonzero_length(A1) > 0 && nonzero_length(A2) > 0 + # TODO: there are too many entries here, this could be further optimized + @plansor y[-1 -2; -3 -4] += GL[-1 7; 6] * x[6 5; 1 3] * A1[7 -2; 5 4] * + A2[4 -4; 3 2] * GR[1 2; -3] + end + + return y +end diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl new file mode 100644 index 000000000..0f93b1d4c --- /dev/null +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -0,0 +1,102 @@ +""" + struct MPO_∂∂ACN{L,O<:Tuple,R} + +Effective local operator obtained from taking the partial derivative of an MPS-MPO-MPS sandwich. +""" +struct MPO_∂∂ACN{L,O<:Tuple,R} <: DerivativeOperator + leftenv::L + operators::O + rightenv::R +end + +Base.length(H::MPO_∂∂ACN) = length(H.operators) + +const MPO_∂∂C{L,R} = MPO_∂∂ACN{L,Tuple{},R} +MPO_∂∂C(GL, GR) = MPO_∂∂ACN(GL, (), GR) + +const MPO_∂∂AC{L,O,R} = MPO_∂∂ACN{L,Tuple{O},R} +MPO_∂∂AC(GL, O, GR) = MPO_∂∂ACN(GL, (O,), GR) + +const MPO_∂∂AC2{L,O₁,O₂,R} = MPO_∂∂ACN{L,Tuple{O₁,O₂},R} +MPO_∂∂AC2(GL, O1, O2, GR) = MPO_∂∂ACN(GL, (O1, O2), GR) + +# Constructors +# ------------ +function ∂∂C(pos::Int, mps, operator, envs) + return MPO_∂∂C(leftenv(envs, pos + 1, mps), rightenv(envs, pos, mps)) +end +function ∂∂AC(pos::Int, mps, operator, envs) + O = isnothing(operator) ? nothing : operator[pos] + return MPO_∂∂AC(leftenv(envs, pos, mps), O, rightenv(envs, pos, mps)) +end +function ∂∂AC2(pos::Int, mps, operator, envs) + O1, O2 = isnothing(operator) ? (nothing, nothing) : (operator[pos], operator[pos + 1]) + return MPO_∂∂AC2(leftenv(envs, pos, mps), O1, O2, rightenv(envs, pos + 1, mps)) +end + +# Properties +# ---------- +function TensorKit.domain(H::MPO_∂∂ACN) + V_l = right_virtualspace(H.leftenv) + V_r = left_virtualspace(H.rightenv) + V_o = prod(physicalspace, H.O; init=one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end +function TensorKit.codomain(H::MPO_∂∂ACN) + V_l = left_virtualspace(H.leftenv) + V_r = right_virtualspace(H.rightenv) + V_o = prod(physicalspace, H.O; init=one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end + +# Actions +# ------- +function (h::MPO_∂∂C{<:MPSBondTensor,<:MPSBondTensor})(x::MPSBondTensor) + @plansor y[-1; -2] ≔ h.leftenv[-1; 1] * x[1; 2] * h.rightenv[2; -2] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_∂∂C{<:MPSTensor,<:MPSTensor})(x::MPSBondTensor) + @plansor y[-1; -2] ≔ h.leftenv[-1 3; 1] * x[1; 2] * h.rightenv[2 3; -2] + return y isa AbstractBlockTensorMap ? only(y) : y +end + +function (h::MPO_∂∂AC{<:MPSTensor,<:MPOTensor,<:MPSTensor})(x::MPSTensor) + @plansor y[-1 -2; -3] ≔ h.leftenv[-1 5; 4] * x[4 2; 1] * + h.operators[1][5 -2; 2 3] * h.rightenv[1 3; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_∂∂AC{<:MPSTensor,<:Number,<:MPSTensor})(x::MPSTensor) + @plansor y[-1 -2; -3] ≔ (h.leftenv[-1 5; 4] * x[4 6; 1] * τ[6 5; 7 -2] * + h.rightenv[1 7; -3]) * only(h.operators) + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_∂∂AC{<:MPSBondTensor,Nothing,<:MPSBondTensor})(x::MPSTensor) + @plansor y[-1 -2; -3] ≔ h.leftenv[-1; 2] * x[2 -2; 1] * h.rightenv[1; -3] +end +function (h::MPO_∂∂AC{<:MPSTensor,<:MPOTensor,<:MPSTensor})(x::GenericMPSTensor{<:Any,3}) + @plansor y[-1 -2 -3; -4] ≔ h.leftenv[-1 7; 6] * x[6 4 2; 1] * + h.operators[1][7 -2; 4 5] * τ[5 -3; 2 3] * + h.rightenv[1 3; -4] + return y isa AbstractBlockTensorMap ? only(y) : y +end + +function (h::MPO_∂∂AC2{<:MPSBondTensor,Nothing,Nothing,<:MPSBondTensor})(x::MPOTensor) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2 -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_∂∂AC2{<:MPSTensor,<:MPOTensor,<:MPOTensor,<:MPSTensor})(x::MPOTensor) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 7; 6] * x[6 5; 1 3] * + h.operators[1][7 -2; 5 4] * h.operators[2][4 -4; 3 2] * + h.rightenv[1 2; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_∂∂AC2{<:MPSTensor,<:MPOTensor,<:MPOTensor,<:MPSTensor})(x::AbstractTensorMap{<:Any, + <:Any, + 3, + 3}) + @plansor y[-1 -2 -3; -4 -5 -6] ≔ h.leftenv[-1 11; 10] * x[10 8 6; 1 2 4] * + h.rightenv[1 3; -4] * + h.operators[1][11 -2; 8 9] * τ[9 -3; 6 7] * + h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] + return y isa AbstractBlockTensorMap ? only(y) : y +end diff --git a/src/algorithms/derivatives/projection_derivatives.jl b/src/algorithms/derivatives/projection_derivatives.jl new file mode 100644 index 000000000..48607460e --- /dev/null +++ b/src/algorithms/derivatives/projection_derivatives.jl @@ -0,0 +1,53 @@ +struct Projection_∂∂ACN{L,O<:Tuple,R} <: DerivativeOperator + leftenv::L + As::O + rightenv::R +end + +const Projection_∂∂AC{L,O,R} = Projection_∂∂ACN{L,Tuple{O},R} +Projection_∂∂AC(GL, A, GR) = Projection_∂∂ACN(GL, (A,), GR) + +const Projection_∂∂AC2{L,O₁,O₂,R} = Projection_∂∂ACN{L,Tuple{O₁,O₂},R} +Projection_∂∂AC2(GL, A1, A2, GR) = Projection_∂∂ACN(GL, (A1, A2), GR) + +struct AC_EffProj{A,L} <: DerivativeOperator + a1::A + le::L + re::L +end +struct AC2_EffProj{A,L} <: DerivativeOperator + a1::A + a2::A + le::L + re::L +end + +# Constructors +# ------------ +function ∂∂AC(pos::Int, state, operator::ProjectionOperator, env) + return Projection_∂∂AC(leftenv(env, pos, state), operator.ket.AC[pos], + rightenv(env, pos, state)) +end +function ∂∂AC2(pos::Int, state, operator::ProjectionOperator, env) + return Projection_∂∂AC2(leftenv(env, pos, state), operator.ket.AC[pos], + operator.ket.AR[pos + 1], + rightenv(env, pos + 1, state)) + return AC2_EffProj(operator.ket.AC[pos], operator.ket.AR[pos + 1], + leftenv(env, pos, state), + rightenv(env, pos + 1, state)) +end + +# Actions +# ------- +function (h::Projection_∂∂AC)(x::MPSTensor) + @plansor v[-1; -2 -3 -4] := h.leftenv[4; -1 -2 5] * h.As[1][5 2; 1] * + h.rightenv[1; -3 -4 3] * conj(x[4 2; 3]) + @plansor y[-1 -2; -3] := conj(v[1; 2 5 6]) * h.leftenv[-1; 1 2 4] * h.As[1][4 -2; 3] * + h.rightenv[3; 5 6 -3] +end +function (h::Projection_∂∂AC2)(x::MPOTensor) + @plansor v[-1; -2 -3 -4] := h.leftenv[6; -1 -2 7] * h.As[1][7 4; 5] * h.As[2][5 2; 1] * + h.rightenv[1; -3 -4 3] * conj(x[6 4; 3 2]) + @plansor y[-1 -2; -3 -4] := conj(v[2; 3 5 6]) * h.leftenv[-1; 2 3 4] * + h.As[1][4 -2; 7] * h.As[2][7 -4; 1] * h.rightenv[1; 5 6 -3] +end diff --git a/src/algorithms/fixedpoint.jl b/src/algorithms/fixedpoint.jl index 37023256f..81e7c69d9 100644 --- a/src/algorithms/fixedpoint.jl +++ b/src/algorithms/fixedpoint.jl @@ -11,7 +11,7 @@ function fixedpoint(A, x₀, which::Symbol, alg::Lanczos) vals, vecs, info = eigsolve(A, x₀, 1, which, alg) if info.converged == 0 - @warn "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])" + @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])" end return vals[1], vecs[1] @@ -21,10 +21,10 @@ function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi) TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg) if info.converged == 0 - @warn "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])" + @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])" end if size(TT, 2) > 1 && TT[2, 1] != 0 - @warn "non-unique fixedpoint detected" + @warnv 1 "non-unique fixedpoint detected" end return vals[1], vecs[1] diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 0ea686faa..9ab21e079 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -10,20 +10,24 @@ $(TYPEDFIELDS) struct DMRG{A,F} <: Algorithm "tolerance for convergence criterium" tol::Float64 + "maximal amount of iterations" maxiter::Int - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" - finalize::F + "setting for how much information is displayed" verbosity::Int + "algorithm used for the eigenvalue solvers" - eigalg::A + alg_eigsolve::A + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F end -function DMRG(; tol=Defaults.tol, maxiter=Defaults.maxiter, eigalg=(;), +function DMRG(; tol=Defaults.tol, maxiter=Defaults.maxiter, alg_eigsolve=(;), verbosity=Defaults.verbosity, finalize=Defaults._finalize) - eigalg′ = eigalg isa NamedTuple ? Defaults.alg_eigsolve(; eigalg...) : eigalg - return DMRG{typeof(eigalg′),typeof(finalize)}(tol, maxiter, finalize, verbosity, - eigalg′) + alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : + alg_eigsolve + return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize) end function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs=environments(ψ, H)) @@ -34,7 +38,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs=environment LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs)) for iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.eigalg, iter, ϵ) + alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) zerovector!(ϵs) for pos in [1:(length(ψ) - 1); length(ψ):-1:2] @@ -70,28 +74,35 @@ Two-site DMRG algorithm for finding the dominant eigenvector. $(TYPEDFIELDS) """ -struct DMRG2{A,F} <: Algorithm +struct DMRG2{A,S,F} <: Algorithm "tolerance for convergence criterium" tol::Float64 + "maximal amount of iterations" maxiter::Int - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" - finalize::F + "setting for how much information is displayed" verbosity::Int "algorithm used for the eigenvalue solvers" - eigalg::A + alg_eigsolve::A + + "algorithm used for the singular value decomposition" + alg_svd::S + "algorithm used for [truncation](@extref TensorKit.tsvd) of the two-site update" trscheme::TruncationScheme + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F end # TODO: find better default truncation -function DMRG2(; tol=Defaults.tol, maxiter=Defaults.maxiter, eigalg=(;), - trscheme=truncerr(1e-6), verbosity=Defaults.verbosity, +function DMRG2(; tol=Defaults.tol, maxiter=Defaults.maxiter, verbosity=Defaults.verbosity, + alg_eigsolve=(;), alg_svd=Defaults.alg_svd(), trscheme, finalize=Defaults._finalize) - eigalg′ = eigalg isa NamedTuple ? Defaults.alg_eigsolve(; eigalg...) : eigalg - return DMRG2{typeof(eigalg′),typeof(finalize)}(tol, maxiter, finalize, verbosity, - eigalg′, trscheme) + alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) : + alg_eigsolve + return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_svd, trscheme, finalize) end function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environments(ψ, H)) @@ -101,7 +112,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environmen LoggingExtras.withlevel(; alg.verbosity) do for iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.eigalg, iter, ϵ) + alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) zerovector!(ϵs) # left to right sweep @@ -110,7 +121,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environmen _, newA2center = fixedpoint(∂∂AC2(pos, ψ, H, envs), ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) @@ -126,7 +137,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environmen _, newA2center = fixedpoint(∂∂AC2(pos, ψ, H, envs), ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) @@ -153,6 +164,6 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environmen return ψ, envs, ϵ end -function find_groundstate(ψ, H, alg::Union{DMRG,DMRG2}, envs...) - return find_groundstate!(copy(ψ), H, alg, envs...) +function find_groundstate(ψ, H, alg::Union{DMRG,DMRG2}, envs...; kwargs...) + return find_groundstate!(copy(ψ), H, alg, envs...; kwargs...) end diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index ca8b0ea9b..551632720 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -85,7 +85,7 @@ Two-site infinite DMRG algorithm for finding the dominant eigenvector. $(TYPEDFIELDS) """ -@kwdef struct IDMRG2{A} <: Algorithm +@kwdef struct IDMRG2{A,S} <: Algorithm "tolerance for convergence criterium" tol::Float64 = Defaults.tol @@ -101,8 +101,11 @@ $(TYPEDFIELDS) "algorithm used for the eigenvalue solvers" alg_eigsolve::A = Defaults.alg_eigsolve() + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + "algorithm used for [truncation](@extref TensorKit.tsvd) of the two-site update" - trscheme::TruncationScheme = truncerr(1e-6) + trscheme::TruncationScheme end function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(ost, H)) @@ -125,7 +128,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os h_ac2 = ∂∂AC2(pos, ψ, H, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AL[pos] = al @@ -143,7 +146,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os h_ac2 = ∂∂AC2(0, ψ, H, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AC[end] = al * c @@ -165,7 +168,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os h_ac2 = ∂∂AC2(pos, ψ, H, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AL[pos] = al @@ -183,7 +186,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os inv(ψ.C[end])[2; 3] * ψ.AC[1][3 -4; -3] h_ac2 = ∂∂AC2(0, ψ, H, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=TensorKit.SVD()) + al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd) normalize!(c) ψ.AR[end] = _transpose_front(inv(ψ.C[end - 1]) * _transpose_tail(al * c)) diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl index 9e3b62034..029eada83 100644 --- a/src/algorithms/timestep/taylorcluster.jl +++ b/src/algorithms/timestep/taylorcluster.jl @@ -78,7 +78,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; # loopback step: Algorithm 1 # constructing the Nth order time evolution MPO - mpo = MPO(parent(H_n)) + mpo = MPO(map(SparseBlockTensorMap, parent(H_n))) for (i, slice) in enumerate(parent(mpo)) V_right = virtual_sz[i + 1] linds_right = linds[i + 1] diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 88fee9298..12c81e9fa 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -114,7 +114,7 @@ $(TYPEDFIELDS) * [Haegeman et al. Phys. Rev. Lett. 107 (2011)](@cite haegeman2011) """ -@kwdef struct TDVP2{A,F} <: Algorithm +@kwdef struct TDVP2{A,S,F} <: Algorithm "algorithm used in the exponential solvers" integrator::A = Defaults.alg_expsolve() @@ -124,8 +124,11 @@ $(TYPEDFIELDS) "maximal amount of iterations for gauging algorithm" gaugemaxiter::Int = Defaults.maxiter + "algorithm used for the singular value decomposition" + alg_svd::S = Defaults.alg_svd() + "algorithm used for truncation of the two-site update" - trscheme::TruncationScheme = truncerr(1e-3) + trscheme::TruncationScheme "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" finalize::F = Defaults._finalize @@ -140,7 +143,7 @@ function timestep!(ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, h_ac2 = ∂∂AC2(i, ψ, H, envs) nac2 = integrate(h_ac2, ac2, t, dt / 2, alg.integrator) - nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=TensorKit.SVD()) + nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=alg.alg_svd) ψ.AC[i] = (nal, complex(nc)) ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) @@ -156,7 +159,7 @@ function timestep!(ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, h_ac2 = ∂∂AC2(i - 1, ψ, H, envs) nac2 = integrate(h_ac2, ac2, t + dt / 2, dt / 2, alg.integrator) - nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=TensorKit.SVD()) + nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=alg.alg_svd) ψ.AC[i - 1] = (nal, complex(nc)) ψ.AC[i] = (complex(nc), _transpose_front(nar)) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 466338541..a004e433b 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -363,6 +363,8 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H)) F_right = fusers[site + 1][k, i, l] j′ = indmap[j, i, l] k′ = indmap[k, i, l] + ((j′ == 1 && k′ == 1) || + (j′ == size(output[site], 1) && k′ == size(output[site], 4))) && continue @plansor o[-1 -2; -3 -4] := h[1 2; -3 6] * F_left[-1; 1 3 5] * conj(F_right[-4; 6 7 8]) * @@ -377,6 +379,8 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H)) F_right = fusers[site + 1][i, l, k] j′ = indmap[i, l, j] k′ = indmap[i, l, k] + ((j′ == 1 && k′ == 1) || + (j′ == size(output[site], 1) && k′ == size(output[site], 4))) && continue @plansor o[-1 -2; -3 -4] := h[1 -2; 3 6] * F_left[-1; 4 2 1] * conj(F_right[-4; 8 7 6]) * @@ -393,6 +397,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H)) if j == 1 F_right = fusers[2][k, end, end] j′ = indmap[k, chi, chi] + j′ == 1 && continue @plansor o[-1 -2; -3 -4] := h[-1 -2; -3 2] * conj(F_right[-4; 2 3 3]) output[1][1, 1, 1, j′] = o end @@ -401,13 +406,13 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H)) if 1 < j < chi F_right = fusers[2][1, j, k] j′ = indmap[1, j, k] + j′ == 1 && continue @plansor o[-1 -2; -3 -4] := h[4 -2; 3 1] * conj(F_right[-4; 6 2 1]) * τ[5 4; 2 3] * τ[-3 -1; 6 5] output[1][1, 1, 1, j′] = o end end - output[1][1, 1, 1, 1] = BraidingTensor{scalartype(H)}(eachspace(output[1])[1]) # ender for (I, h) in nonzero_pairs(H[end]) @@ -415,6 +420,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H)) if k > 1 F_left = fusers[end][j, k, chi] k′ = indmap[j, k, chi] + k′ == size(output[end], 1) && continue @plansor o[-1 -2; -3 -4] := F_left[-1; 1 2 6] * h[1 3; -3 4] * τ[3 2; 4 5] * τ[5 6; -4 -2] output[end][k′, 1, 1, 1] = o diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 21adff0be..753b691c9 100644 --- a/src/operators/abstractmpo.jl +++ b/src/operators/abstractmpo.jl @@ -37,11 +37,6 @@ end # Utility functions # ----------------- -function jordanmpotensortype(::Type{S}, ::Type{T}) where {S<:VectorSpace,T<:Number} - TT = Base.promote_typejoin(tensormaptype(S, 2, 2, T), BraidingTensor{T,S}) - return SparseBlockTensorMap{TT} -end - remove_orphans!(mpo::AbstractMPO; tol=eps(real(scalartype(mpo)))^(3 / 4)) = mpo function remove_orphans!(mpo::SparseMPO; tol=eps(real(scalartype(mpo)))^(3 / 4)) droptol!.(mpo, tol) @@ -203,6 +198,11 @@ end function _fuse_mpo_mpo(O1::MPOTensor, O2::MPOTensor, Fₗ, Fᵣ) return if O1 isa BraidingTensor && O2 isa BraidingTensor + # shouldn't happen + T = promote_type(scalartype(O1), scalartype(O2)) + V = fuse(left_virtualspace(O2) ⊗ left_virtualspace(O1)) ⊗ physicalspace(O1) ← + physicalspace(O2) ⊗ fuse(right_virtualspace(O2) ⊗ right_virtualspace(O1)) + return BraidingTensor{T}(V) elseif O1 isa BraidingTensor @plansor O′[-1 -2; -3 -4] := Fₗ[-1; 1 2] * O2[1 3; -3 5] * diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl new file mode 100644 index 000000000..862cfc3a1 --- /dev/null +++ b/src/operators/jordanmpotensor.jl @@ -0,0 +1,343 @@ +""" + JordanMPOTensor{E,S,TA,TB,TC,TD} <: AbstractBlockTensorMap{E,S,2,2} + +A tensor map that represents the upper triangular block form of a matrix product operator (MPO). + +```math +\\begin{pmatrix} +1 & C & D \\\\ +0 & A & B \\\\ +0 & 0 & 1 +\\end{pmatrix} +``` +""" +struct JordanMPOTensor{E,S, + TA<:AbstractTensorMap{E,S,2,2}, + TB<:AbstractTensorMap{E,S,2,1}, + TC<:AbstractTensorMap{E,S,1,2}, + TD<:AbstractTensorMap{E,S,1,1}} <: AbstractBlockTensorMap{E,S,2,2} + V::TensorMapSumSpace{S,2,2} + A::SparseBlockTensorMap{TA,E,S,2,2,4} + B::SparseBlockTensorMap{TB,E,S,2,1,3} + C::SparseBlockTensorMap{TC,E,S,1,2,3} + D::SparseBlockTensorMap{TD,E,S,1,1,2} + # uninitialized constructor + function JordanMPOTensor{E,S,TA,TB,TC,TD}(::UndefInitializer, + V::TensorMapSumSpace{S,2,2}) where {E,S,TA,TB, + TC,TD} + allVs = eachspace(V) + + # Note that this is a bit of a hack using end to get the last index: + # it should be 1 or end depending on this being an "edge" tensor or a "bulk" tensor + VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)]) + A = SparseBlockTensorMap{TA}(undef, VA) + + VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4) + B = SparseBlockTensorMap{TB}(undef, VB) + + VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1) + C = SparseBlockTensorMap{TC}(undef, VC) + + VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1) + D = SparseBlockTensorMap{TD}(undef, VD) + + return new{E,S,TA,TB,TC,TD}(V, A, B, C, D) + end + + # constructor from data + function JordanMPOTensor{E,S,TA,TB,TC,TD}(V::TensorMapSumSpace, + A::SparseBlockTensorMap{TA,E,S,2,2}, + B::SparseBlockTensorMap{TB,E,S,2,1}, + C::SparseBlockTensorMap{TC,E,S,1,2}, + D::SparseBlockTensorMap{TD,E,S,1,1}) where {E, + S, + TA, + TB, + TC, + TD} + return new{E,S,TA,TB,TC,TD}(V, A, B, C, D) + end +end + +function JordanMPOTensor{E,S}(::UndefInitializer, V::TensorMapSumSpace{S}) where {E,S} + return jordanmpotensortype(S, E)(undef, V) +end +function JordanMPOTensor{E}(::UndefInitializer, V::TensorMapSumSpace{S}) where {E,S} + return JordanMPOTensor{E,S}(undef, V) +end + +function JordanMPOTensor(V::TensorMapSumSpace{S,2,2}, + A::SparseBlockTensorMap{TA,E,S,2,2}, + B::SparseBlockTensorMap{TB,E,S,2,1}, + C::SparseBlockTensorMap{TC,E,S,1,2}, + D::SparseBlockTensorMap{TD,E,S,1,1}) where {E,S,TA,TB,TC,TD} + allVs = eachspace(V) + VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)]) + VA == space(A) || throw(SpaceMismatch("A-block has incompatible spaces")) + + VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4) + VB == space(B) || throw(SpaceMismatch("B-block has incompatible spaces")) + + VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1) + VC == space(C) || throw(SpaceMismatch("C-block has incompatible spaces")) + + VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1) + VD == space(D) || throw(SpaceMismatch("D-block has incompatible spaces")) + + return JordanMPOTensor{E,S,TA,TB,TC,TD}(V, A, B, C, D) +end + +function JordanMPOTensor(W::SparseBlockTensorMap{TT,E,S,2,2}) where {TT,E,S} + @assert W[1, 1, 1, 1] isa BraidingTensor && W[end, 1, 1, end] isa BraidingTensor + # @assert all(I -> I[1] ≤ I[4], nonzero_keys(W)) + + A = W[2:(end - 1), 1, 1, 2:(end - 1)] + B = W[2:(end - 1), 1, 1, end] + C = W[1, 1, 1, 2:(end - 1)] + D = W[1, 1, 1, end:end] # ensure still blocktensor to allow for sparse + + return JordanMPOTensor(space(W), + A, + removeunit(B, 4), + removeunit(C, 1), + removeunit(removeunit(D, 4), 1)) +end + +function jordanmpotensortype(::Type{S}, ::Type{E}) where {S<:VectorSpace,E<:Number} + TA = Union{tensormaptype(S, 2, 2, E),BraidingTensor{E,S}} + TB = tensormaptype(S, 2, 1, E) + TC = tensormaptype(S, 1, 2, E) + TD = tensormaptype(S, 1, 1, E) + return JordanMPOTensor{E,S,TA,TB,TC,TD} +end +function jordanmpotensortype(::Type{O}) where {O<:MPOTensor} + return jordanmpotensortype(spacetype(O), scalartype(O)) +end + +function Base.similar(W::JordanMPOTensor, ::Type{T}) where {T<:Number} + return JordanMPOTensor{T}(undef, space(W)) +end + +# Properties +# ---------- +TensorKit.space(W::JordanMPOTensor) = W.V +Base.eltype(::Type{JordanMPOTensor{E,S,TA,TB,TC,TD}}) where {E,S,TA,TB,TC,TD} = TA + +function Base.haskey(W::JordanMPOTensor, I::CartesianIndex{4}) + Base.checkbounds(W, I.I...) + # only has braiding tensors if sizes are large enough + sz = size(W) + (sz[1] > 1 && I == CartesianIndex(1, 1, 1, 1) || + sz[4] > 1 && I == CartesianIndex(sz[1], 1, 1, sz[4])) && return true + + row, col = I.I[1], I.I[4] + + if row == 1 && col == sz[4] + return haskey(W.D, CartesianIndex(1, 1)) + elseif row == 1 + return haskey(W.C, CartesianIndex(1, 1, col - 1)) + elseif col == sz[4] + return haskey(W.B, CartesianIndex(row - 1, 1, 1)) + elseif 1 < row < sz[1] && 1 < col < sz[4] + return haskey(W.A, CartesianIndex(row - 1, 1, 1, col - 1)) + else + return false + end +end + +Base.parent(W::JordanMPOTensor) = parent(SparseBlockTensorMap(W)) + +BlockTensorKit.issparse(W::JordanMPOTensor) = true + +# Converters +# ---------- +function SparseBlockTensorMap(W::JordanMPOTensor) + τ = BraidingTensor{scalartype(W)}(eachspace(W)[1]) + W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W),spacetype(W),2,2}}(undef_blocks, + space(W)) + if size(W, 1) > 1 + W′[1, 1, 1, 1] = τ + end + if size(W, 4) > 1 + W′[end, 1, 1, end] = τ + end + + Ia = CartesianIndex(1, 0, 0, 1) + for (I, v) in nonzero_pairs(W.A) + W′[I + Ia] = v + end + Ib = CartesianIndex(1, 0, 0) + for (I, v) in nonzero_pairs(W.B) + W′[I + Ib, size(W′, 4)] = insertrightunit(v, 3) + end + Ic = CartesianIndex(0, 0, 1) + for (I, v) in nonzero_pairs(W.C) + W′[1, I + Ic] = insertleftunit(v, 1) + end + W′[1, 1, 1, end] = insertrightunit(insertleftunit(only(W.D), 1), 3) + + return W′ +end + +# Indexing +# -------- + +@inline Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4}) = W[I.I...] +@propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int,4}) + @assert I[2] == I[3] == 1 + i = I[1] + j = I[4] + if (size(W, 4) > 1 && i == 1 && j == 1) || + (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4)) + return BraidingTensor{scalartype(W)}(eachspace(W)[1]) + elseif i == 1 && j == size(W, 4) + return insertrightunit(insertleftunit(only(W.D), 1), 3) + elseif i == 1 + return insertleftunit(W.C[1, 1, j - 1], 1) + elseif j == size(W, 4) + return insertrightunit(W.B[i - 1, 1, 1], 3) + elseif 1 < i < size(W, 1) && 1 < j < size(W, 4) + return W.A[i - 1, 1, 1, j - 1] + else + return zeros(scalartype(W), eachspace(W)[i, 1, 1, j]) + end +end + +@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::CartesianIndex{4}) + return setindex!(W, v, I.I...) +end +@propagate_inbounds function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, + I::Vararg{Int,4}) + @assert I[2] == I[3] == 1 + i = I[1] + j = I[4] + if i == 1 && j == size(W, 4) + W.D[1] = removeunit(removeunit(v, 4), 1) + elseif i == 1 && 1 < j < size(W, 4) + W.C[1, 1, j - 1] = removeunit(v, 1) + elseif j == size(W, 4) && 1 < i < size(W, 1) + W.B[i - 1, 1, 1] = removeunit(v, 4) + elseif 1 < i < size(W, 1) && 1 < j < size(W, 4) + W.A[i - 1, 1, 1, j - 1] = v + elseif (size(W, 4) > 1 && i == 1 && j == 1) || + (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4)) + v isa BraidingTensor || throw(ArgumentError("Cannot set BraidingTensor")) + else + throw(ArgumentError("Cannot set index ($i, 1, 1, $j)")) + end + return W +end +@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::Int) + return setindex!(W, v, CartesianIndices(W)[I]) +end + +# Sparse functionality +# -------------------- +function BlockTensorKit.nonzero_keys(W::JordanMPOTensor) + nrows = size(W, 1) + ncols = size(W, 4) + p = CartesianIndex{4}[] + ncols > 1 && push!(p, CartesianIndex(1, 1, 1, 1)) + nrows > 1 && push!(p, CartesianIndex(nrows, 1, 1, ncols)) + + Ia = CartesianIndex(1, 0, 0, 1) + for I in nonzero_keys(W.A) + push!(p, I + Ia) + end + + Ib = CartesianIndex(1, 0, 0) + for I in nonzero_keys(W.B) + push!(p, CartesianIndex((I + Ib).I..., ncols)) + end + + Ic = CartesianIndex(0, 0, 1) + for I in nonzero_keys(W.C) + push!(p, CartesianIndex(1, (I + Ic).I...)) + end + + for I in nonzero_keys(W.D) + push!(p, CartesianIndex(1, 1, 1, ncols)) + end + + return p +end +function BlockTensorKit.nonzero_values(W::JordanMPOTensor) + return Iterators.map(I -> W[I], nonzero_keys(W)) +end +function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor) + return Iterators.map(I -> I => W[I], nonzero_keys(W)) +end +function BlockTensorKit.nonzero_length(W::JordanMPOTensor) + return nonzero_length(W.A) + nonzero_length(W.B) + nonzero_length(W.C) + + nonzero_length(W.D) + Int(size(W, 1) > 1) + Int(size(W, 4) > 1) +end + +# linalg +# ------ +# do we want this? +function Base.:+(W1::JordanMPOTensor, W2::JordanMPOTensor) + return SparseBlockTensorMap(W1) + SparseBlockTensorMap(W2) +end +function Base.:-(W1::JordanMPOTensor, W2::JordanMPOTensor) + return SparseBlockTensorMap(W1) - SparseBlockTensorMap(W2) +end + +function fuse_mul_mpo(O1::JordanMPOTensor, O2::JordanMPOTensor) + TT = promote_type((eltype(O1)), eltype((O2))) + V = fuse(left_virtualspace(O2) ⊗ left_virtualspace(O1)) ⊗ physicalspace(O1) ← + physicalspace(O2) ⊗ fuse(right_virtualspace(O2) ⊗ right_virtualspace(O1)) + O = jordanmpotensortype(TT)(undef, V) + cartesian_inds = reshape(CartesianIndices(O), + size(O2, 1), size(O1, 1), + size(O, 2), size(O, 3), + size(O2, 4), size(O1, 4)) + for (I, o2) in nonzero_pairs(O2), (J, o1) in nonzero_pairs(O1) + K = cartesian_inds[I[1], J[1], I[2], I[3], I[4], J[4]] + O[K] = fuse_mul_mpo(o1, o2) + end + return O +end + +function _conj_mpo(W::JordanMPOTensor) + V = left_virtualspace(W)' ⊗ physicalspace(W) ← physicalspace(W) ⊗ right_virtualspace(W)' + A = _conj_mpo(W.A) + @plansor B[-1 -2; -3] ≔ conj(W.B[-1 -3; -2]) + @plansor C[-1; -2 -3] ≔ conj(W.C[-2; -1 -3]) + D = copy(adjoint(W.D)) + return JordanMPOTensor(V, A, B, C, D) +end + +# Utility +# ------- +function Base.copy(W::JordanMPOTensor) + return JordanMPOTensor(W.V, copy(W.A), copy(W.B), copy(W.C), copy(W.D)) +end +function Base.copy!(Wdst::JordanMPOTensor, Wsrc::JordanMPOTensor) + space(Wdst) == space(Wsrc) || throw(SpaceMismatch()) + copy!(Wdst.A, Wsrc.A) + copy!(Wdst.B, Wsrc.B) + copy!(Wdst.C, Wsrc.C) + copy!(Wdst.D, Wsrc.D) + return Wdst +end + +# Avoid falling back to `norm(W1 - W2)` which has to convert to SparseBlockTensorMap +function Base.isapprox(W1::JordanMPOTensor, W2::JordanMPOTensor; kwargs...) + return isapprox(W1.A, W2.A; kwargs...) && + isapprox(W1.B, W2.B; kwargs...) && + isapprox(W1.C, W2.C; kwargs...) && + isapprox(W1.D, W2.D; kwargs...) +end + +function Base.summary(io::IO, W::JordanMPOTensor) + szstring = Base.dims2string(size(W)) + TT = eltype(W) + typeinfo = get(io, :typeinfo, Any) + if typeinfo <: typeof(W) || typeinfo <: TT + typestring = "" + else + typestring = "{$TT}" + end + V = space(W) + return print(io, "$szstring JordanMPOTensor$typestring($V)") +end diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index aef1c8d98..347b3bd48 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -28,7 +28,7 @@ H = MPOHamiltonian(lattice, (i, i+1) => O for i in 1:length(lattice)-1) See also [`instantiate_operator`](@ref), which is responsable for instantiating the local operators in a form that is compatible with this constructor. """ -struct MPOHamiltonian{TO,V<:AbstractVector{TO}} <: AbstractMPO{TO} +struct MPOHamiltonian{TO<:JordanMPOTensor,V<:AbstractVector{TO}} <: AbstractMPO{TO} W::V end @@ -58,14 +58,15 @@ end FiniteMPOHamiltonian(Ws::Vector{<:Matrix}) Create a `FiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]` represents -the the operator at site `i`, left level `j` and right level `k`. Here, the entries can be +the operator at site `i`, left level `j` and right level `k`. Here, the entries can be either `MPOTensor`, `Missing` or `Number`. """ function FiniteMPOHamiltonian(Ws::Vector{<:Matrix}) T = promote_type(_split_mpoham_types.(Ws)...) - return FiniteMPOHamiltonian{T}(Ws) + W = jordanmpotensortype(T) + return FiniteMPOHamiltonian{W}(Ws) end -function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor} +function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:JordanMPOTensor} T = scalartype(O) L = length(W_mats) # initialize sumspaces @@ -87,7 +88,7 @@ function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor} # ability to deduce spaces for (site, W_mat) in enumerate(W_mats) # physical space - operator_id = findfirst(x -> x isa O, W_mat) + operator_id = findfirst(x -> x isa MPOTensor, W_mat) @assert !isnothing(operator_id) "could not determine physical space at site $site" Pspaces[site] = physicalspace(W_mat[operator_id]) @@ -148,7 +149,8 @@ either `MPOTensor`, `Missing` or `Number`. """ function InfiniteMPOHamiltonian(Ws::Vector{<:Matrix}) T = promote_type(_split_mpoham_types.(Ws)...) - return InfiniteMPOHamiltonian{T}(Ws) + TW = jordanmpotensortype(T) + return InfiniteMPOHamiltonian{TW}(Ws) end function InfiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor} # InfiniteMPOHamiltonian only works for square matrices: @@ -167,7 +169,7 @@ function InfiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor # physical spaces Pspaces = map(W_mats) do W_mat - operator_id = findfirst(x -> x isa O, W_mat) + operator_id = findfirst(x -> x isa MPOTensor, W_mat) @assert !isnothing(operator_id) "could not determine physical space" return physicalspace(W_mat[operator_id]) end @@ -380,9 +382,9 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, E = scalartype(T) S = spacetype(T) - virtualspaces = Vector{SumSpace{S}}(undef, length(lattice) + 1) - virtualspaces[1] = SumSpace(oneunit(S)) - virtualspaces[end] = SumSpace(oneunit(S)) + virtualsumspaces = Vector{SumSpace{S}}(undef, length(lattice) + 1) + virtualsumspaces[1] = SumSpace(fill(oneunit(S), 1)) + virtualsumspaces[end] = SumSpace(fill(oneunit(S), 1)) for i in 1:(length(lattice) - 1) n_channels = maximum(last, nonzero_keys[i]; init=1) + 1 @@ -390,34 +392,28 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, if n_channels > 2 for ((key_L, key_R), O) in zip(nonzero_keys[i], nonzero_opps[i]) V[key_R == 0 ? end : key_R] = if O isa Number - virtualspaces[i][key_L] + virtualsumspaces[i][key_L] else right_virtualspace(O) end end end - virtualspaces[i + 1] = V + virtualsumspaces[i + 1] = V end - Otype = jordanmpotensortype(S, E) + # construct the tensor + TW = jordanmpotensortype(S, E) Os = map(1:length(lattice)) do site - # Initialize blocktensor - O = Otype(undef, virtualspaces[site] * lattice[site], - lattice[site] * virtualspaces[site + 1]) - fill!(O, zero(E)) - if site != length(lattice) - O[1, 1, 1, 1] = BraidingTensor{E}(eachspace(O)[1, 1, 1, 1]) - end - if site != 1 - O[end, end, end, end] = BraidingTensor{E}(eachspace(O)[end, end, end, end]) - end + V = virtualsumspaces[site] * lattice[site] ← + lattice[site] * virtualsumspaces[site + 1] + O = TW(undef, V) # Fill it - for ((key_L, key_R), o) in zip(nonzero_keys[site], nonzero_opps[site]) - key_R′ = key_R == 0 ? length(virtualspaces[site + 1]) : key_R - O[key_L, 1, 1, key_R′] += if o isa Number + for ((key_L, key_R′), o) in zip(nonzero_keys[site], nonzero_opps[site]) + key_R = key_R′ == 0 ? length(virtualsumspaces[site + 1]) : key_R′ + O[key_L, 1, 1, key_R] += if o isa Number iszero(o) && continue - τ = BraidingTensor{E}(eachspace(O)[key_L, 1, 1, key_R′]) + τ = BraidingTensor{E}(eachspace(O)[key_L, 1, 1, key_R]) isone(o) ? τ : τ * o else o @@ -522,13 +518,12 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, return SumSpace(collect(S, V)) end - # construct the tensors - Otype = jordanmpotensortype(S, E) + # construct the tensor + TW = jordanmpotensortype(S, E) Os = map(1:length(lattice)) do site - O = Otype(undef, virtualsumspaces[site - 1] * lattice[site], - lattice[site] * virtualsumspaces[site]) - O[1, 1, 1, 1] = BraidingTensor{E}(eachspace(O)[1, 1, 1, 1]) - O[end, end, end, end] = BraidingTensor{E}(eachspace(O)[end, end, end, end]) + V = virtualsumspaces[site - 1] * lattice[site] ← + lattice[site] * virtualsumspaces[site] + O = TW(undef, V) # Fill it for ((key_L, key_R′), o) in zip(nonzero_keys[site], nonzero_opps[site]) @@ -584,15 +579,19 @@ function Base.getproperty(H::MPOHamiltonian, sym::Symbol) end end -function isidentitylevel(H::InfiniteMPOHamiltonian, i::Int) - isemptylevel(H, i) && return false - return all(parent(H)) do h - return (h[i, 1, 1, i] isa BraidingTensor) +function isidentitylevel(H::InfiniteMPOHamiltonian{<:JordanMPOTensor}, i::Int) + if i == 1 || i == size(H[1], 1) + return true + else + return all(H.A) do A + return haskey(A, CartesianIndex(i - 1, 1, 1, i - 1)) && + A[i - 1, 1, 1, i - 1] isa BraidingTensor + end end end function isemptylevel(H::InfiniteMPOHamiltonian, i::Int) return any(parent(H)) do h - return !(CartesianIndex(i, 1, 1, i) in nonzero_keys(h)) + return !haskey(h, CartesianIndex(i, 1, 1, i)) end end @@ -603,6 +602,17 @@ function Base.convert(::Type{TensorMap}, H::FiniteMPOHamiltonian) return convert(TensorMap, _instantiate_finitempo(L, M, R)) end +function Base.convert(::Type{FiniteMPOHamiltonian{O1}}, + H::FiniteMPOHamiltonian{O2}) where {O1,O2} + O1 === O2 && return H + return FiniteMPOHamiltonian(convert.(O1, parent(H))) +end +function Base.convert(::Type{InfiniteMPOHamiltonian{O1}}, + H::InfiniteMPOHamiltonian{O2}) where {O1,O2} + O1 === O2 && return H + return InfiniteMPOHamiltonian(convert.(O1, parent(H))) +end + function add_physical_charge(H::MPOHamiltonian, charges::AbstractVector{<:Sector}) W = map(add_physical_charge, parent(H), charges) if isfinite(H) @@ -640,81 +650,48 @@ end # Linear Algebra # -------------- - -function Base.:+(H₁::MPOH, H₂::MPOH) where {MPOH<:MPOHamiltonian} - check_length(H₁, H₂) - @assert all(physicalspace.(parent(H₁)) .== physicalspace.(parent(H₂))) "physical spaces should match" - isinf = MPOH <: InfiniteMPOHamiltonian - +function Base.:+(H₁::FiniteMPOHamiltonian{O}, + H₂::FiniteMPOHamiltonian{O}) where {O<:JordanMPOTensor} + N = check_length(H₁, H₂) H = similar(parent(H₁)) - for i in 1:length(H) - # instantiate new blocktensor - Vₗ₁ = left_virtualspace(H₁, i) - Vₗ₂ = left_virtualspace(H₂, i) - @assert Vₗ₁[1] == Vₗ₂[1] && Vₗ₁[end] == Vₗ₂[end] "trivial spaces should match" - Vₗ = (!isinf && i == 1) ? Vₗ₁ : BlockTensorKit.oplus(Vₗ₁[1:(end - 1)], Vₗ₂[2:end]) - - Vᵣ₁ = right_virtualspace(H₁, i) - Vᵣ₂ = right_virtualspace(H₂, i) - @assert Vᵣ₁[1] == Vᵣ₂[1] && Vᵣ₁[end] == Vᵣ₂[end] "trivial spaces should match" - Vᵣ = (!isinf && i == length(H)) ? Vᵣ₁ : - BlockTensorKit.oplus(Vᵣ₁[1:(end - 1)], Vᵣ₂[2:end]) - - W = similar(eltype(H), Vₗ ⊗ physicalspace(H₁, i) ← physicalspace(H₁, i) ⊗ Vᵣ) - #= - (Direct) sum of two hamiltonians in Jordan form: - (1 C₁ D₁) (1 C₂ D₂) (1 C₁ C₂ D₁+D₂) - (0 A₁ B₁) + (0 A₂ B₂) = (0 A₁ 0 B₁ ) - (0 0 1 ) (0 0 1 ) (0 0 A₂ B₂ ) - (0 0 0 1 ) - =# - fill!(W, zero(scalartype(W))) - W[1, 1, 1, 1] = BraidingTensor{scalartype(W)}(eachspace(W)[1, 1, 1, 1]) - W[end, 1, 1, end] = BraidingTensor{scalartype(W)}(eachspace(W)[end, 1, 1, end]) - - if H₁ isa InfiniteMPOHamiltonian || i != length(H) - C₁ = H₁.C[i] - C₁_inds = CartesianIndices((1:1, 1:1, 1:1, 2:(size(H₁[i], 4) - 1))) - copyto!(W, C₁_inds, C₁, CartesianIndices(C₁)) - - C₂ = H₂.C[i] - C₂_inds = CartesianIndices((1:1, 1:1, 1:1, size(H₁[i], 4):(size(W, 4) - 1))) - copyto!(W, C₂_inds, C₂, CartesianIndices(C₂)) - end - - D₁ = H₁.D[i] - D₂ = H₂.D[i] - W[1, 1, 1, end] = D₁ + D₂ + Vtriv = oneunit(spacetype(H₁)) - if H₁ isa InfiniteMPOHamiltonian || i != 1 && i != length(H) - A₁ = H₁.A[i] - A₁_inds = CartesianIndices((2:(size(H₁[i], 1) - 1), 1:1, 1:1, - 2:(size(H₁[i], 4) - 1))) - copyto!(W, A₁_inds, A₁, CartesianIndices(A₁)) + for i in 1:N + A = cat(H₁[i].A, H₂[i].A; dims=(1, 4)) + B = cat(H₁[i].B, H₂[i].B; dims=1) + C = cat(H₁[i].C, H₂[i].C; dims=3) + D = H₁[i].D + H₂[i].D - A₂ = H₂.A[i] - A₂_inds = CartesianIndices((size(H₁[i], 1):(size(W, 1) - 1), 1:1, 1:1, - size(H₁[i], 4):(size(W, 4) - 1))) - copyto!(W, A₂_inds, A₂, CartesianIndices(A₂)) - end + Vleft = i == 1 ? left_virtualspace(H₁, 1) : + BlockTensorKit.oplus(Vtriv, left_virtualspace(A), Vtriv) + Vright = i == N ? right_virtualspace(H₁, N) : + BlockTensorKit.oplus(Vtriv, right_virtualspace(A), Vtriv) + V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright - if H₁ isa InfiniteMPOHamiltonian || i != 1 - B₁ = H₁.B[i] - B₁_inds = CartesianIndices((2:(size(H₁[i], 1) - 1), 1:1, 1:1, - size(W, 4):size(W, 4))) - copyto!(W, B₁_inds, B₁, CartesianIndices(B₁)) + H[i] = eltype(H)(V, A, B, C, D) + end + return FiniteMPOHamiltonian(H) +end +function Base.:+(H₁::InfiniteMPOHamiltonian{O}, + H₂::InfiniteMPOHamiltonian{O}) where {O<:JordanMPOTensor} + N = check_length(H₁, H₂) + H = similar(parent(H₁)) + Vtriv = oneunit(spacetype(H₁)) + for i in 1:N + A = cat(H₁[i].A, H₂[i].A; dims=(1, 4)) + B = cat(H₁[i].B, H₂[i].B; dims=1) + C = cat(H₁[i].C, H₂[i].C; dims=3) + D = H₁[i].D + H₂[i].D - B₂ = H₂.B[i] - B₂_inds = CartesianIndices((size(H₁[i], 1):(size(W, 1) - 1), 1:1, 1:1, - size(W, 4):size(W, 4))) - copyto!(W, B₂_inds, B₂, CartesianIndices(B₂)) - end + Vleft = BlockTensorKit.oplus(Vtriv, left_virtualspace(A), Vtriv) + Vright = BlockTensorKit.oplus(Vtriv, right_virtualspace(A), Vtriv) + V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright - H[i] = W + H[i] = eltype(H)(V, A, B, C, D) end - - return H₁ isa FiniteMPOHamiltonian ? FiniteMPOHamiltonian(H) : InfiniteMPOHamiltonian(H) + return InfiniteMPOHamiltonian(H) end + function Base.:+(H::FiniteMPOHamiltonian, λs::AbstractVector{<:Number}) check_length(H, λs) lattice = [physicalspace(H, i) for i in 1:length(H)] @@ -739,46 +716,24 @@ Base.:-(H::MPOHamiltonian, λs::AbstractVector{<:Number}) = H + (-λs) Base.:-(λs::AbstractVector{<:Number}, H::MPOHamiltonian) = λs + (-H) Base.:-(H1::MPOHamiltonian, H2::MPOHamiltonian) = H1 + (-H2) -function VectorInterface.scale!(H::InfiniteMPOHamiltonian, λ::Number) - foreach(parent(H)) do h - # multiply scalar with start of every interaction - # this avoids double counting - # 2:end to avoid multiplying the top left and bottom right corners - return scale!(h[1, 1, 1, 2:end], λ) - end - return H -end - -function VectorInterface.scale!(H::FiniteMPOHamiltonian, λ::Number) - foreach(enumerate(parent(H))) do (i, h) - if i != length(H) - scale!(h[1, 1, 1, 2:end], λ) # multiply top row (except BraidingTensor) - else - scale!(h[1, 1, 1, end], λ) # multiply right column (except BraidingTensor) - end +function VectorInterface.scale!(H::MPOHamiltonian{O}, + λ::Number) where {O<:JordanMPOTensor} + for i in 1:length(H) + scale!(H[i].C, λ) + scale!(H[i].D, λ) end return H end - -function VectorInterface.scale!(dst::MPOHamiltonian, src::MPOHamiltonian, - λ::Number) - N = check_length(dst, src) +function VectorInterface.scale!(Hdst::MPOHamiltonian{<:JordanMPOTensor}, + Hsrc::MPOHamiltonian{<:JordanMPOTensor}, λ::Number) + N = check_length(Hdst, Hsrc) for i in 1:N - space(dst[i]) == space(src[i]) || throw(SpaceMismatch()) - zerovector!(dst[i]) - for (I, v) in nonzero_pairs(src[i]) - # only scale "starting" terms - isstarting = I[1] == 1 && - ((isfinite(dst) && i == N && I[4] == size(src[i], 4)) || - ((!isfinite(dst) || i != N) && I[4] > 1)) - if v isa BraidingTensor && !isstarting - dst[i][I] = v - else - dst[i][I] = scale!(dst[i][I], v, isstarting ? λ : One()) - end - end + scale!(Hdst[i].C, Hsrc[i].C, λ) + scale!(Hdst[i].D, Hsrc[i].D, λ) + copy!(Hdst[i].A, Hsrc[i].A) + copy!(Hdst[i].B, Hsrc[i].B) end - return dst + return Hdst end function Base.:*(H1::MPOHamiltonian, H2::MPOHamiltonian) diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index ecd74e0fb..71b97a97d 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -8,6 +8,7 @@ module Defaults import KrylovKit: GMRES, Arnoldi, Lanczos using OhMyThreads using ..MPSKit: DynamicTol +using TensorKit: TensorKit const VERBOSE_NONE = 0 const VERBOSE_WARN = 1 @@ -53,6 +54,8 @@ function alg_eigsolve(; ishermitian=true, tol=tol, maxiter=maxiter, verbosity=0, return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end +alg_svd() = TensorKit.SDD() + # TODO: make verbosity and maxiter actually do something function alg_environments(; tol=tol, maxiter=maxiter, verbosity=0, dynamic_tols=dynamic_tols, tol_min=tol_min, tol_max=tol_max, @@ -60,6 +63,7 @@ function alg_environments(; tol=tol, maxiter=maxiter, verbosity=0, alg = (; tol, maxiter, verbosity) return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end + function alg_expsolve(; tol=tol, maxiter=maxiter, verbosity=0, ishermitian=true, krylovdim=krylovdim) return ishermitian ? Lanczos(; tol, maxiter, krylovdim, verbosity) : diff --git a/test/algorithms.jl b/test/algorithms.jl index f45319133..73d29ca71 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -322,7 +322,7 @@ end @testset "timestep" verbose = true begin dt = 0.1 - algs = [TDVP(), TDVP2()] + algs = [TDVP(), TDVP2(; trscheme=truncdim(10))] L = 10 H = force_planar(heisenberg_XXX(; spin=1 // 2, L)) @@ -387,7 +387,7 @@ end @testset "time_evolve" verbose = true begin t_span = 0:0.1:0.1 - algs = [TDVP(), TDVP2()] + algs = [TDVP(), TDVP2(; trscheme=truncdim(10))] L = 10 H = force_planar(heisenberg_XXX(; spin=1 // 2, L)) diff --git a/test/bench.jl b/test/bench.jl new file mode 100644 index 000000000..765c940bd --- /dev/null +++ b/test/bench.jl @@ -0,0 +1,45 @@ +using TestEnv +TestEnv.activate() + +using Profile + +include("setup.jl") +using .TestSetup +using TensorKit +using BlockTensorKit +using MPSKit +using MPSKit: ∂∂AC2 +using BenchmarkTools +using KrylovKit + +L = 50 +g = 1 +D = 128 +H = heisenberg_XXX(; L) +psi = FiniteMPS(randn, ComplexF64, physicalspace(H), fill(ComplexSpace(D), L - 1)) +envs = environments(psi, H); + +eigalg = Lanczos(; tol=1e-14, krylovdim=5, maxiter=1, verbosity=0, eager=false) +@benchmark find_groundstate(psi, H, + DMRG2(; eigalg, trscheme=truncdim(D), maxiter=3, verbosity=0); + svd_alg=TensorKit.SVD()) +@benchmark find_groundstate(psi, H, + DMRG2(; eigalg, trscheme=truncdim(D), maxiter=3, verbosity=0); + svd_alg=TensorKit.SDD()) +@profview find_groundstate(psi, H, + DMRG2(; eigalg, trscheme=truncdim(D), maxiter=3, verbosity=0)); + +pos = L ÷ 2 +Heff = ∂∂AC2(pos, psi, H, envs); +@plansor x0[-1 -2; -3 -4] := psi.AC[pos][-1 -2; 1] * psi.AR[pos + 1][1 -4; -3]; + +@benchmark $Heff * $x0 + +H2 = FiniteMPO(map(SparseBlockTensorMap, parent(H))) +envs2 = environments(psi, H2); +H2eff = ∂∂AC2(pos, psi, H2, envs2); +@plansor x0[-1 -2; -3 -4] := psi.AC[pos][-1 -2; 1] * psi.AR[pos + 1][1 -4; -3]; + +@benchmark $H2eff * $x0 + +VSCodeServer.view_profile() diff --git a/test/infinite.jl b/test/infinite.jl new file mode 100644 index 000000000..adc787cf3 --- /dev/null +++ b/test/infinite.jl @@ -0,0 +1,409 @@ +using TestEnv; +TestEnv.activate() +include("setup.jl") + +using .TestSetup + +using Test, TestExtras +using MPSKit +using MPSKit: fuse_mul_mpo +using TensorKit +using TensorKit: ℙ +using BlockTensorKit + +verbosity_full = 5 +verbosity_conv = 1 +tol = 1e-8 +g = 4.0 +D = 6 + +H_ref = transverse_field_ising(; g) +H2 = InfiniteMPOHamiltonian(map(SparseBlockTensorMap, parent(H_ref))) + +ψ = InfiniteMPS(ℂ^2, ℂ^D) +v₀ = variance(ψ, H_ref) + +for unit_cell_size in [1, 3] + ψ = unit_cell_size == 1 ? InfiniteMPS(ℂ^2, ℂ^D) : repeat(ψ, unit_cell_size) + H = repeat(H_ref, unit_cell_size) + H3 = repeat(H2, unit_cell_size) + + # test logging + ψ, envs, δ = find_groundstate(ψ, H, + VUMPS(; tol, verbosity=verbosity_full, maxiter=2)) + + ψ, envs, δ = find_groundstate(ψ, H, VUMPS(; tol, verbosity=verbosity_conv)) + v = variance(ψ, H, envs) + v2 = expectation_value(ψ, H3) + expectation_value(ψ, H, envs) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ + @test v < 1e-2 +end + +for unit_cell_size in [1, 3] + ψ = unit_cell_size == 1 ? InfiniteMPS(ℂ^2, ℂ^D) : repeat(ψ, unit_cell_size) + H = repeat(H_ref, unit_cell_size) + + # test logging + ψ, envs, δ = find_groundstate(ψ, H, + IDMRG(; tol, verbosity=verbosity_full, maxiter=2)) + + ψ, envs, δ = find_groundstate(ψ, H, IDMRG(; tol, verbosity=verbosity_conv)) + v = variance(ψ, H, envs) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ + @test v < 1e-2 +end + +begin + ψ = repeat(InfiniteMPS(2, D), 2) + H = repeat(H_ref, 2) + + trscheme = truncbelow(1e-8) + + # test logging + ψ, envs, δ = find_groundstate(ψ, H, + IDMRG2(; tol, verbosity=verbosity_full, maxiter=2, + trscheme)) + + ψ, envs, δ = find_groundstate(ψ, H, + IDMRG2(; tol, verbosity=verbosity_conv, trscheme)) + v = variance(ψ, H, envs) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ + @test v < 1e-2 +end + +for unit_cell_size in [1, 3] + ψ = unit_cell_size == 1 ? InfiniteMPS(ℂ^2, ℂ^D) : repeat(ψ, unit_cell_size) + H = repeat(H_ref, unit_cell_size) + + # test logging + ψ, envs, δ = find_groundstate(ψ, H, + GradientGrassmann(; tol, verbosity=verbosity_full, + maxiter=2)) + + ψ, envs, δ = find_groundstate(ψ, H, + GradientGrassmann(; tol, verbosity=verbosity_conv)) + v = variance(ψ, H, envs) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ + @test v < 1e-2 +end + +for unit_cell_size in [1, 3] + ψ = unit_cell_size == 1 ? InfiniteMPS(ℂ^2, ℂ^D) : repeat(ψ, unit_cell_size) + H = repeat(H_ref, unit_cell_size) + + alg = VUMPS(; tol=100 * tol, verbosity=verbosity_conv, maxiter=10) & + GradientGrassmann(; tol, verbosity=verbosity_conv, maxiter=50) + ψ, envs, δ = find_groundstate(ψ, H, alg) + + v = variance(ψ, H, envs) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ + @test v < 1e-2 +end + +######################## + +tol = 1e-8 +g = 4.0 +D = 6 +L = 10 + +H = transverse_field_ising(; g, L) + +H + H + +begin + ψ₀ = FiniteMPS(randn, ComplexF64, L, ℂ^2, ℂ^D) + v₀ = variance(ψ₀, H) + + # test logging + ψ, envs, δ = find_groundstate(ψ₀, H, + DMRG(; verbosity=verbosity_full, maxiter=2)) + + ψ, envs, δ = find_groundstate(ψ, H, + DMRG(; verbosity=verbosity_conv, maxiter=10), + envs) + v = variance(ψ, H) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ && v < 1e-2 +end + +begin + ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℂ^2, ℂ^D) + v₀ = variance(ψ₀, H) + trscheme = truncdim(floor(Int, D * 1.5)) + # test logging + ψ, envs, δ = find_groundstate(ψ₀, H, + DMRG2(; verbosity=verbosity_full, maxiter=2, + trscheme)) + + ψ, envs, δ = find_groundstate(ψ, H, + DMRG2(; verbosity=verbosity_conv, maxiter=10, + trscheme), envs) + v = variance(ψ, H) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + + @test v < v₀ && v < 1e-2 +end + +begin + ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℂ^2, ℂ^D) + v₀ = variance(ψ₀, H) + + # test logging + ψ, envs, δ = find_groundstate(ψ₀, H, + GradientGrassmann(; verbosity=verbosity_full, + maxiter=2)) + + ψ, envs, δ = find_groundstate(ψ, H, + GradientGrassmann(; verbosity=verbosity_conv, + maxiter=50), + envs) + v = variance(ψ, H) + + # test using low variance + @test sum(δ) ≈ 0 atol = 1e-3 + @test v < v₀ && v < 1e-2 +end + +L = 3 +T = ComplexF64 +V = ℂ^2 +lattice = fill(V, L) +O₁ = rand(T, V, V) +E = id(storagetype(O₁), domain(O₁)) +O₂ = rand(T, V^2 ← V^2) + +H1 = FiniteMPOHamiltonian(lattice, i => O₁ for i in 1:L) +H2 = FiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:(L - 1)) +H3 = FiniteMPOHamiltonian(lattice, 1 => O₁, (2, 3) => O₂, (1, 3) => O₂) + +# check if constructor works by converting back to tensormap +H1_tm = convert(TensorMap, H1) +operators = vcat(fill(E, L - 1), O₁) +@test H1_tm ≈ mapreduce(+, 1:L) do i + return reduce(⊗, circshift(operators, i)) +end +operators = vcat(fill(E, L - 2), O₂) +@test convert(TensorMap, H2) ≈ mapreduce(+, 1:(L - 1)) do i + return reduce(⊗, circshift(operators, i)) +end +@test convert(TensorMap, H3) ≈ + O₁ ⊗ E ⊗ E + E ⊗ O₂ + permute(O₂ ⊗ E, ((1, 3, 2), (4, 6, 5))) + +# check if adding terms on the same site works +single_terms = Iterators.flatten(Iterators.repeated((i => O₁ / 2 for i in 1:L), 2)) +H4 = FiniteMPOHamiltonian(lattice, single_terms) +@test H4 ≈ H1 atol = 1e-6 +double_terms = Iterators.flatten(Iterators.repeated(((i, i + 1) => O₂ / 2 + for i in 1:(L - 1)), 2)) +H5 = FiniteMPOHamiltonian(lattice, double_terms) +@test H5 ≈ H2 atol = 1e-6 + +# test linear algebra +@test H1 ≈ + FiniteMPOHamiltonian(lattice, 1 => O₁) + + FiniteMPOHamiltonian(lattice, 2 => O₁) + + FiniteMPOHamiltonian(lattice, 3 => O₁) +@test 0.8 * H1 + 0.2 * H1 ≈ H1 atol = 1e-6 +@test convert(TensorMap, H1 + H2) ≈ convert(TensorMap, H1) + convert(TensorMap, H2) atol = 1e-6 + +# test dot and application +state = rand(T, prod(lattice)) +mps = FiniteMPS(state) + +@test convert(TensorMap, H1 * mps) ≈ H1_tm * state +@test H1 * state ≈ H1_tm * state +@test dot(mps, H2, mps) ≈ dot(mps, H2 * mps) + +# test constructor from dictionary with mixed linear and Cartesian lattice indices as keys +grid = square = fill(V, 3, 3) + +local_operators = Dict((I,) => O₁ for I in eachindex(grid)) +I_vertical = CartesianIndex(1, 0) +vertical_operators = Dict((I, I + I_vertical) => O₂ + for I in eachindex(IndexCartesian(), square) + if I[1] < size(square, 1)) +operators = merge(local_operators, vertical_operators) +H4 = FiniteMPOHamiltonian(grid, operators) + +@test H4 ≈ + FiniteMPOHamiltonian(grid, local_operators) + + FiniteMPOHamiltonian(grid, vertical_operators) + +pspace = ℂ^2 +Dspace = ℂ^D +Os = map(1:3) do i + O = rand(ComplexF64, pspace^i, pspace^i) + return O += O' +end +fs = [t -> 1t, 1, 1] + +L = 5 +ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace) +lattice = fill(pspace, L) +Hs = map(enumerate(Os)) do (i, O) + return FiniteMPOHamiltonian(lattice, + ntuple(x -> x + j, i) => O for j in 0:(L - i)) +end +summedH = LazySum(Hs) + +envs = map(H -> environments(ψ, H), Hs) +summed_envs = environments(ψ, summedH) + +expval = sum(zip(Hs, envs)) do (H, env) + return expectation_value(ψ, H, env) +end +expval1 = expectation_value(ψ, sum(summedH)) +expval2 = expectation_value(ψ, summedH, summed_envs) +expval3 = expectation_value(ψ, summedH) +@test expval ≈ expval1 +@test expval ≈ expval2 +@test expval ≈ expval3 + +# test derivatives +summedhct = MPSKit.∂∂C(1, ψ, summedH, summed_envs) +sum1 = sum(zip(Hs, envs)) do (H, env) + return MPSKit.∂∂C(1, ψ, H, env)(ψ.C[1]) +end +@test summedhct(ψ.C[1], 0.0) ≈ sum1 + +summedhct = MPSKit.∂∂AC(1, ψ, summedH, summed_envs) +sum2 = sum(zip(Hs, envs)) do (H, env) + return MPSKit.∂∂AC(1, ψ, H, env)(ψ.AC[1]) +end +@test summedhct(ψ.AC[1], 0.0) ≈ sum2 + +v = _transpose_front(ψ.AC[1]) * _transpose_tail(ψ.AR[2]) +summedhct = MPSKit.∂∂AC2(1, ψ, summedH, summed_envs) +sum3 = sum(zip(Hs, envs)) do (H, env) + return MPSKit.∂∂AC2(1, ψ, H, env)(v) +end +@test summedhct(v, 0.0) ≈ sum3 + +Hts = [MultipliedOperator(Hs[1], fs[1]), MultipliedOperator(Hs[2], fs[2]), Hs[3]] +summedH = LazySum(Hts) +t = 1.1 +summedH_at = summedH(t) +@which summedH(t) + +envs = map(H -> environments(ψ, H), Hs) +summed_envs = environments(ψ, summedH) + +expval = sum(zip(fs, Hs, envs)) do (f, H, env) + return (f isa Function ? f(t) : f) * expectation_value(ψ, H, env) +end +expval1 = expectation_value(ψ, sum(summedH_at)) +expval2 = expectation_value(ψ, summedH_at, summed_envs) +expval3 = expectation_value(ψ, summedH_at) +@test expval ≈ expval1 +@test expval ≈ expval2 +@test expval ≈ expval3 + +# test derivatives +summedhct = MPSKit.∂∂C(1, ψ, summedH, summed_envs) +sum1 = sum(zip(fs, Hs, envs)) do (f, H, env) + if f isa Function + f = f(t) + end + return f * MPSKit.∂∂C(1, ψ, H, env)(ψ.C[1]) +end +@test summedhct(ψ.C[1], t) ≈ sum1 + +summedhct = MPSKit.∂∂AC(1, ψ, summedH, summed_envs) +sum2 = sum(zip(fs, Hs, envs)) do (f, H, env) + if f isa Function + f = f(t) + end + return f * MPSKit.∂∂AC(1, ψ, H, env)(ψ.AC[1]) +end +@test summedhct(ψ.AC[1], t) ≈ sum2 + +v = _transpose_front(ψ.AC[1]) * _transpose_tail(ψ.AR[2]) +summedhct = MPSKit.∂∂AC2(1, ψ, summedH, summed_envs) +sum3 = sum(zip(fs, Hs, envs)) do (f, H, env) + return (f isa Function ? f(t) : f) * MPSKit.∂∂AC2(1, ψ, H, env)(v) +end +@test summedhct(v, t) ≈ sum3 +# ------------------------------------------------------ +dt = 0.1 +algs = [TDVP(), TDVP2()] +L = 10 + +H = (TestSetup.heisenberg_XXX(; spin=1 // 2, L)) +ψ₀ = FiniteMPS(L, ℂ^2, ℂ^1) +E₀ = expectation_value(ψ₀, H) + +@testset "Finite $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, H, 0.0, dt, alg) + E = expectation_value(ψ, H, envs) + @test E₀ ≈ E atol = 1e-2 +end + +Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + +@testset "Finite LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, Hlazy, 0.0, dt, alg) + E = expectation_value(ψ, Hlazy, envs) + @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1e-2 +end + +Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45); + +@testset "Finite TimeDependent LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in + algs + ψ, envs = timestep(ψ₀, Ht(1.0), 0.0, dt, alg) + E = expectation_value(ψ, Ht(1.0), envs) + + ψt, envst = timestep(ψ₀, Ht, 1.0, dt, alg) + Et = expectation_value(ψt, Ht(1.0), envst) + @test E ≈ Et atol = 1e-8 +end + +H = repeat(TestSetup.heisenberg_XXX(; spin=1), 2) +ψ₀ = InfiniteMPS([ℂ^3, ℂ^3], [ℂ^50, ℂ^50]) +E₀ = expectation_value(ψ₀, H) + +@testset "Infinite TDVP" begin + ψ, envs = timestep(ψ₀, H, 0.0, dt, TDVP()) + E = expectation_value(ψ, H, envs) + @test E₀ ≈ E atol = 1e-2 +end + +Hlazy = LazySum([3 * deepcopy(H), 1.55 * deepcopy(H), -0.1 * deepcopy(H)]) + +@testset "Infinite LazySum TDVP" begin + ψ, envs = timestep(ψ₀, Hlazy, 0.0, dt, TDVP()) + E = expectation_value(ψ, Hlazy, envs) + @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1e-2 +end + +Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) + +@testset "Infinite TimeDependent LazySum" begin + ψ, envs = timestep(ψ₀, Ht(1.0), 0.0, dt, TDVP()) + E = expectation_value(ψ, Ht(1.0), envs) + + ψt, envst = timestep(ψ₀, Ht, 1.0, dt, TDVP()) + Et = expectation_value(ψt, Ht(1.0), envst) + @test E ≈ Et atol = 1e-8 +end diff --git a/test/setup.jl b/test/setup.jl index 804ca6e3a..9e0a570c1 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -4,6 +4,7 @@ module TestSetup # imports using MPSKit +using MPSKit: JordanMPOTensor using TensorKit using TensorKit: PlanarTrivial, ℙ, BraidingTensor using BlockTensorKit @@ -57,6 +58,18 @@ function force_planar(x::SparseBlockTensorMap) data = Dict(I => force_planar(v) for (I, v) in pairs(x.data)) return SparseBlockTensorMap{valtype(data)}(data, force_planar(space(x))) end +function force_planar(W::JordanMPOTensor) + V = force_planar(space(W)) + TW = MPSKit.jordanmpotensortype(eltype(V[1]), scalartype(W)) + dst = TW(undef, V) + + for t in (:A, :B, :C, :D) + for (I, v) in nonzero_pairs(getproperty(W, t)) + getproperty(dst, t)[I] = force_planar(v) + end + end + return dst +end force_planar(mpo::MPOHamiltonian) = MPOHamiltonian(map(force_planar, parent(mpo))) force_planar(mpo::MPO) = MPO(map(force_planar, parent(mpo))) diff --git a/test/states.jl b/test/states.jl index 4a7f1fee3..5c694da01 100644 --- a/test/states.jl +++ b/test/states.jl @@ -176,7 +176,7 @@ end @test real(e2) ≤ real(e1) - window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(), envs) + window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme=truncdim(20)), envs) window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs) e3 = expectation_value(window, (1, 2) => O)