Refactor derivatives and add JordanMPOTensor

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"

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")

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

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)))

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)))

src/algorithms/changebonds/svdcut.jl

@@ -7,17 +7,20 @@ 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...)
     return changebonds!(copy(ψ), alg; 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

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