Refactor InfiniteEnvironments and GradientGrassmann

src/MPSKit.jl

@@ -106,10 +106,9 @@ include("transfermatrix/transfer.jl")
 
 include("environments/abstract_envs.jl")
 include("environments/finite_envs.jl")
-include("environments/infinitempo_envs.jl")
-include("environments/infinitempohamiltonian_envs.jl")
+include("environments/infinite_envs.jl")
+include("environments/multiline_envs.jl")
 include("environments/qp_envs.jl")
-include("environments/idmrg_envs.jl")
 include("environments/multiple_envs.jl")
 include("environments/lazylincocache.jl")
 
@@ -145,6 +144,7 @@ include("algorithms/excitation/dmrgexcitation.jl")
 include("algorithms/excitation/chepigaansatz.jl")
 include("algorithms/excitation/exci_transfer_system.jl")
 
+include("algorithms/statmech/leading_boundary.jl")
 include("algorithms/statmech/vumps.jl")
 include("algorithms/statmech/vomps.jl")
 include("algorithms/statmech/gradient_grassmann.jl")

src/algorithms/approximate/approximate.jl

@@ -25,3 +25,35 @@ of an MPS `ψ₀`.
 - `VOMPS`: Tangent space method for truncating uniform MPS.
 """
 approximate, approximate!
+
+# implementation in terms of Multiline
+function approximate(ψ::InfiniteMPS,
+                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS},
+                     algorithm,
+                     envs=environments(ψ, toapprox))
+    envs′ = Multiline([envs])
+    multi, envs = approximate(convert(MultilineMPS, ψ),
+                              (convert(MultilineMPO, toapprox[1]),
+                               convert(MultilineMPS, toapprox[2])), algorithm, envs′)
+    ψ = convert(InfiniteMPS, multi)
+    return ψ, envs
+end
+
+# dispatch to in-place method
+function approximate(ψ, toapprox, alg::Union{DMRG,DMRG2,IDMRG1,IDMRG2},
+                     envs=environments(ψ, toapprox))
+    return approximate!(copy(ψ), toapprox, alg, envs)
+end
+
+# disambiguate
+function approximate(ψ::InfiniteMPS,
+                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS},
+                     algorithm::Union{IDMRG1,IDMRG2},
+                     envs=environments(ψ, toapprox))
+    envs′ = Multiline([envs])
+    multi, envs = approximate(convert(MultilineMPS, ψ),
+                              (convert(MultilineMPO, toapprox[1]),
+                               convert(MultilineMPS, toapprox[2])), algorithm, envs′)
+    ψ = convert(InfiniteMPS, multi)
+    return ψ, envs
+end

src/algorithms/approximate/fvomps.jl

@@ -1,15 +1,5 @@
-# dispatch to in-place method
-function approximate(ψ, toapprox, alg::Union{DMRG,DMRG2}, envs...)
-    return approximate!(copy(ψ), toapprox, alg, envs...)
-end
-
-function approximate!(ψ::AbstractFiniteMPS, sq, alg, envs=environments(ψ, sq))
-    tor = approximate!(ψ, [sq], alg, [envs])
-    return (tor[1], tor[2][1], tor[3])
-end
-
-function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
-                      envs=[environments(ψ, sq) for sq in squash])
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2,
+                      envs=environments(ψ, Oϕ))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG2")
 
@@ -18,9 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
         for iter in 1:(alg.maxiter)
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = sum(zip(squash, envs)) do (sq, pr)
-                    return ac2_proj(pos, ψ, pr)
-                end
+                AC2′ = ac2_proj(pos, ψ, Oϕ, envs)
                 al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme)
 
                 AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
@@ -31,7 +19,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
             end
 
             # finalize
-            ψ, envs = alg.finalize(iter, ψ, squash, envs)::Tuple{typeof(ψ),typeof(envs)}
+            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ),typeof(envs)}
 
             if ϵ < alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
@@ -48,8 +36,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
     return ψ, envs, ϵ
 end
 
-function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG,
-                      envs=[environments(ψ, sq) for sq in squash])
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs=environments(ψ, Oϕ))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG")
 
@@ -58,18 +45,15 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG,
         for iter in 1:(alg.maxiter)
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = sum(zip(squash, envs)) do (sq, pr)
-                    return ac_proj(pos, ψ, pr)
-                end
-
+                AC′ = ac_proj(pos, ψ, Oϕ, envs)
                 AC = ψ.AC[pos]
                 ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
 
                 ψ.AC[pos] = AC′
             end
 
             # finalize
-            ψ, envs = alg.finalize(iter, ψ, squash, envs)::Tuple{typeof(ψ),typeof(envs)}
+            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ),typeof(envs)}
 
             if ϵ < alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)

src/algorithms/approximate/idmrg.jl

@@ -1,8 +1,5 @@
-function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
-                     alg::IDMRG1, oenvs=environments(ost, toapprox))
-    ψ = copy(ost)
-    mpo, above = toapprox
-    envs = IDMRGEnvironments(ost, oenvs)
+function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
+                      alg::IDMRG1, envs=environments(ψ, toapprox))
     log = IterLog("IDMRG")
     ϵ::Float64 = 2 * alg.tol
 
@@ -12,31 +9,27 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             C_current = ψ.C[:, 0]
 
             # left to right sweep
-            for col in 1:size(ψ, 2), row in 1:size(ψ, 1)
-                h = MPO_∂∂AC(mpo[row, col], leftenv(envs, row, col),
-                             rightenv(envs, row, col))
-                ψ.AC[row + 1, col] = h * above.AC[row, col]
-                normalize!(ψ.AC[row + 1, col])
-
-                ψ.AL[row + 1, col], ψ.C[row + 1, col] = leftorth(ψ.AC[row + 1, col])
-
-                tm = TransferMatrix(above.AL[row, col], mpo[row, col], ψ.AL[row + 1, col])
-                setleftenv!(envs, row, col + 1, normalize(leftenv(envs, row, col) * tm))
+            for col in 1:size(ψ, 2)
+                for row in 1:size(ψ, 1)
+                    ψ.AC[row + 1, col] = ac_proj(row, col, ψ, toapprox, envs)
+                    normalize!(ψ.AC[row + 1, col])
+                    ψ.AL[row + 1, col], ψ.C[row + 1, col] = leftorth!(ψ.AC[row + 1, col])
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
             end
 
             # right to left sweep
-            for col in size(ψ, 2):-1:1, row in 1:size(ψ, 1)
-                h = MPO_∂∂AC(mpo[row, col], leftenv(envs, row, col),
-                             rightenv(envs, row, col))
-                ψ.AC[row + 1, col] = h * above.AC[row, col]
-                normalize!(ψ.AC[row + 1, col])
-
-                ψ.C[row + 1, col - 1], temp = rightorth(_transpose_tail(ψ.AC[row + 1, col]))
-                ψ.AR[row + 1, col] = _transpose_front(temp)
-
-                tm = TransferMatrix(above.AR[row, col], mpo[row, col], ψ.AR[row + 1, col])
-                setrightenv!(envs, row, col - 1, normalize(tm * rightenv(envs, row, col)))
+            for col in size(ψ, 2):-1:1
+                for row in 1:size(ψ, 1)
+                    ψ.AC[row + 1, col] = ac_proj(row, col, ψ, toapprox, envs)
+                    normalize!(ψ.AC[row + 1, col])
+                    ψ.C[row + 1, col - 1], temp = rightorth!(_transpose_tail(ψ.AC[row + 1,
+                                                                                  col]))
+                    ψ.AR[row + 1, col] = _transpose_front(temp)
+                end
+                transfer_rightenv!(envs, ψ, toapprox, col - 1)
             end
+            normalize!(envs, ψ, toapprox)
 
             ϵ = norm(C_current - ψ.C[:, 0])
 
@@ -52,72 +45,62 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
         end
     end
 
+    # TODO: immediately compute in-place
     nst = MultilineMPS(map(x -> x, ψ.AR); tol=alg.tol_gauge)
-    nenvs = environments(nst, toapprox)
-    return nst, nenvs, ϵ
+    copy!(ψ, nst) # ensure output destination is unchanged
+
+    recalculate!(envs, ψ, toapprox)
+    return ψ, envs, ϵ
 end
 
-function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
-                     alg::IDMRG2, oenvs=environments(ost, toapprox))
-    length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    mpo, above = toapprox
-    ψ = copy(ost)
-    envs = IDMRGEnvironments(ost, oenvs)
+function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
+                      alg::IDMRG2, envs=environments(ψ, toapprox))
+    size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2"))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("IDMRG2")
+    O, ϕ = toapprox
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ)
         for iter in 1:(alg.maxiter)
             C_current = ψ.C[:, 0]
 
             # sweep from left to right
-            for col in 1:size(ψ, 2), row in 1:size(ψ, 1)
-                ac2 = above.AC[row, col] * _transpose_tail(above.AR[row, col + 1])
-                h = MPO_∂∂AC2(mpo[row, col], mpo[row, col + 1], leftenv(envs, row, col),
-                              rightenv(envs, row, col + 1))
-
-                al, c, ar, = tsvd!(h * ac2; trunc=alg.trscheme, alg=TensorKit.SVD())
-                normalize!(c)
-
-                ψ.AL[row + 1, col] = al
-                ψ.C[row + 1, col] = complex(c)
-                ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-
-                setleftenv!(envs, row, col + 1,
-                            normalize(leftenv(envs, row, col) *
-                                      TransferMatrix(above.AL[row, col], mpo[row, col],
-                                                     ψ.AL[row + 1, col])))
-                setrightenv!(envs, row, col,
-                             normalize(TransferMatrix(above.AR[row, col + 1],
-                                                      mpo[row, col + 1],
-                                                      ψ.AR[row + 1, col + 1]) *
-                                       rightenv(envs, row, col + 1)))
+            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())
+                    normalize!(c)
+
+                    ψ.AL[row + 1, col] = al
+                    ψ.C[row + 1, col] = complex(c)
+                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                    ψ.AC[row + 1, col + 1] = _transpose_front(c * ar)
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
+                transfer_rightenv!(envs, ψ, toapprox, col)
             end
+            normalize!(envs, ψ, toapprox)
 
             # sweep from right to left
-            for col in (size(ψ, 2) - 1):-1:0, row in 1:size(ψ, 1)
-                ac2 = above.AL[row, col] * _transpose_tail(above.AC[row, col + 1])
-                h = MPO_∂∂AC2(mpo[row, col], mpo[row, col + 1], leftenv(envs, row, col),
-                              rightenv(envs, row, col + 1))
-
-                al, c, ar, = tsvd!(h * ac2; trunc=alg.trscheme, alg=TensorKit.SVD())
-                normalize!(c)
-
-                ψ.AL[row + 1, col] = al
-                ψ.C[row + 1, col] = complex(c)
-                ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-
-                setleftenv!(envs, row, col + 1,
-                            normalize(leftenv(envs, row, col) *
-                                      TransferMatrix(above.AL[row, col], mpo[row, col],
-                                                     ψ.AL[row + 1, col])))
-                setrightenv!(envs, row, col,
-                             normalize(TransferMatrix(above.AR[row, col + 1],
-                                                      mpo[row, col + 1],
-                                                      ψ.AR[row + 1, col + 1]) *
-                                       rightenv(envs, row, col + 1)))
+            for col in (size(ψ, 2) - 1):-1:0
+                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())
+                    normalize!(c)
+
+                    ψ.AL[row + 1, col] = al
+                    ψ.C[row + 1, col] = complex(c)
+                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
+                transfer_rightenv!(envs, ψ, toapprox, col)
             end
+            normalize!(envs, ψ, toapprox)
 
             # update error
             ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
@@ -139,7 +122,10 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
         end
     end
 
+    # TODO: immediately compute in-place
     nst = MultilineMPS(map(x -> x, ψ.AR); tol=alg.tol_gauge)
-    nenvs = environments(nst, toapprox)
-    return nst, nenvs, ϵ
+    copy!(ψ, nst) # ensure output destination is unchanged
+    recalculate!(envs, ψ, toapprox)
+
+    return ψ, envs, ϵ
 end

src/algorithms/approximate/vomps.jl

@@ -1,14 +1,3 @@
-function approximate(ψ::InfiniteMPS,
-                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS}, algorithm,
-                     envs=environments(ψ, toapprox))
-    # PeriodicMPO's always act on MultilineMPS's. To avoid code duplication, define everything in terms of MultilineMPS's.
-    multi, envs = approximate(convert(MultilineMPS, ψ),
-                              (convert(MultilineMPO, toapprox[1]),
-                               convert(MultilineMPS, toapprox[2])), algorithm, envs)
-    ψ = convert(InfiniteMPS, multi)
-    return ψ, envs
-end
-
 Base.@deprecate(approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
                             alg::VUMPS, envs...; kwargs...),
                 approximate(ψ, toapprox,
@@ -18,10 +7,12 @@
 
 function approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
                      alg::VOMPS, envs=environments(ψ, toapprox))
-    ϵ::Float64 = calc_galerkin(ψ, envs)
+    ϵ::Float64 = calc_galerkin(ψ, toapprox..., envs)
     temp_ACs = similar.(ψ.AC)
     scheduler = Defaults.scheduler[]
     log = IterLog("VOMPS")
+    alg_environments = updatetol(alg.alg_environments, 0, ϵ)
+    recalculate!(envs, ψ, toapprox...; alg_environments.tol)
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ)
@@ -34,11 +25,11 @@
             ψ = MultilineMPS(temp_ACs, ψ.C[:, end]; alg_gauge.tol, alg_gauge.maxiter)
 
             alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-            recalculate!(envs, ψ; alg_environments.tol)
+            recalculate!(envs, ψ, toapprox...; alg_environments.tol)
 
             ψ, envs = alg.finalize(iter, ψ, toapprox, envs)::Tuple{typeof(ψ),typeof(envs)}
 
-            ϵ = calc_galerkin(ψ, envs)
+            ϵ = calc_galerkin(ψ, toapprox..., envs)
 
             if ϵ <= alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
@@ -55,18 +46,20 @@
     return ψ, envs, ϵ
 end
 
-function _vomps_localupdate(loc, ψ, (O, ψ₀), envs, factalg=QRpos())
+function _vomps_localupdate(loc, ψ, Oϕ, envs, factalg=QRpos())
     local tmp_AC, tmp_C
     if Defaults.scheduler[] isa SerialScheduler
-        tmp_AC = circshift([ac_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
-        tmp_C = circshift([c_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+        tmp_AC = circshift([ac_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)], 1)
+        tmp_C = circshift([c_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)], 1)
     else
         @sync begin
             Threads.@spawn begin
-                tmp_AC = circshift([ac_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+                tmp_AC = circshift([ac_proj(row, loc, ψ, Oϕ, envs)
+                                    for row in 1:size(ψ, 1)], 1)
             end
             Threads.@spawn begin
-                tmp_C = circshift([c_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+                tmp_C = circshift([c_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)],
+                                  1)
             end
         end
     end

src/algorithms/changebonds/changebonds.jl

@@ -9,6 +9,14 @@ See also: [`SvdCut`](@ref), [`RandExpand`](@ref), [`VUMPSSvdCut`](@ref), [`Optim
 function changebonds end
 function changebonds! end
 
+# write in terms of MultilineMPS
+function changebonds(ψ::InfiniteMPS, operator::InfiniteMPO, alg,
+                     envs=environments(ψ, operator))
+    ψ′, envs′ = changebonds(convert(MultilineMPS, ψ), convert(MultilineMPO, operator), alg,
+                            Multiline([envs]))
+    return convert(InfiniteMPS, ψ′), envs
+end
+
 _expand(ψ, AL′, AR′) = _expand!(copy(ψ), AL′, AR′)
 function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector)
     for i in 1:length(ψ)