Pass through more settings for environment and gauging sub-algorithms

docs/src/changelog.md

@@ -23,6 +23,17 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Changed
 
+- `environments` now follows a single positional contract for every state and operator kind:
+  `environments(below, operator, above, alg)`, where `alg` is the environment algorithm
+  (slot 4). The operator form requires an explicit `above`. Auxiliary inputs are keyword-only:
+  `leftstart`/`rightstart` for finite and window environments, and `lenvs`/`renvs` for window
+  and quasiparticle environments.
+  The two-argument form `environments(below, above)` (with two states) is reserved for the
+  operator-free overlap environments. There is no two-argument `environments(below, operator)`
+  shorthand: a two-argument call is always the overlap, since the second argument cannot be
+  disambiguated between a ket and an operator (undecidable for density matrices, where states
+  and operators share a representation). ([#436](https://github.com/QuantumKitHub/MPSKit.jl/pull/436))
+
 ### Deprecated
 
 ### Removed

src/MPSKit.jl

@@ -178,7 +178,6 @@ 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")
 include("algorithms/statmech/idmrg.jl")

src/algorithms/ED.jl

@@ -68,7 +68,7 @@ function exact_diagonalization(
     TB = tensormaptype(spacetype(H), 1, 1, T)
     Cs = Vector{Union{Missing, TB}}(missing, L + 1)
     state = FiniteMPS(ALs, ARs, ACs, Cs)
-    envs = environments(state, H)
+    envs = environments(state, H, state)
 
     # optimize the middle site
     # Because the MPS is full rank - this is equivalent to the full Hamiltonian

src/algorithms/approximate/approximate.jl

@@ -1,13 +1,17 @@
 @doc """
     approximate(ψ₀, (O, ψ), algorithm, [environments]; kwargs...) -> (ψ, environments)
     approximate!(ψ₀, (O, ψ), algorithm, [environments]; kwargs...) -> (ψ, environments)
+    approximate(ψ₀, ψ, algorithm, [environments]; kwargs...) -> (ψ, environments)
+    approximate!(ψ₀, ψ, algorithm, [environments]; kwargs...) -> (ψ, environments)
 
 Compute an approximation to the application of an operator `O` to the state `ψ` in the form
-of an MPS `ψ₀`.
+of an MPS `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, `ψ₀` is
+approximated directly to `ψ` (i.e. `O` is taken to be the identity).
 
 ## Arguments
 - `ψ₀::AbstractMPS`: initial guess of the approximated state
 - `(O::AbstractMPO, ψ::AbstractMPS)`: operator `O` and state `ψ` to be approximated
+- `ψ::AbstractMPS`: state to be approximated directly (without an operator)
 - `algorithm`: approximation algorithm. See below for a list of available algorithms.
 - `[environments]`: MPS environment manager
 
@@ -26,10 +30,15 @@ of an MPS `ψ₀`.
 """
 approximate, approximate!
 
+# the trailing `environments` arguments for an operator/ket bundle:
+# a tuple carries an explicit operator (3-argument form), a bare state means overlap (2-argument form).
+_environment_args(Oϕ::Tuple) = Oϕ
+_environment_args(ϕ) = (ϕ,)
+
 # implementation in terms of Multiline
 function approximate(
         ψ::InfiniteMPS, toapprox::Tuple{<:InfiniteMPO, <:InfiniteMPS}, algorithm,
-        envs = environments(ψ, toapprox)
+        envs = environments(ψ, toapprox...)
     )
     envs′ = Multiline([envs])
     multi, envs, δ = approximate(
@@ -43,16 +52,15 @@ end
 
 # dispatch to in-place method
 function approximate(
-        ψ, toapprox, alg::Union{DMRG, DMRG2, IDMRG, IDMRG2},
-        envs = environments(ψ, toapprox)
+        ψ, toapprox, alg::Union{DMRG, DMRG2, IDMRG, IDMRG2}, envs...
     )
-    return approximate!(copy(ψ), toapprox, alg, envs)
+    return approximate!(copy(ψ), toapprox, alg, envs...)
 end
 
 # disambiguate
 function approximate(
         ψ::InfiniteMPS, toapprox::Tuple{<:InfiniteMPO, <:InfiniteMPS},
-        algorithm::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox)
+        algorithm::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox...)
     )
     envs′ = Multiline([envs])
     multi, envs, δ = approximate(

src/algorithms/approximate/fvomps.jl

@@ -1,4 +1,4 @@
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, Oϕ))
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG2")
 
@@ -35,7 +35,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment
     return ψ, envs, ϵ
 end
 
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, Oϕ))
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG")
 

src/algorithms/approximate/idmrg.jl

@@ -1,6 +1,6 @@
 function approximate!(
         ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS}, alg::IDMRG,
-        envs = environments(ψ, toapprox)
+        envs = environments(ψ, toapprox...)
     )
     log = IterLog("IDMRG")
     ϵ::Float64 = 2 * alg.tol
@@ -62,7 +62,7 @@ end
 
 function approximate!(
         ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS},
-        alg::IDMRG2, envs = environments(ψ, toapprox)
+        alg::IDMRG2, envs = environments(ψ, toapprox...)
     )
     size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2"))
     ϵ::Float64 = 2 * alg.tol

src/algorithms/approximate/vomps.jl

@@ -12,13 +12,13 @@ Base.@deprecate(
 
 function approximate(
         mps::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS}, alg::VOMPS,
-        envs = environments(mps, toapprox)
+        envs = environments(mps, toapprox...)
     )
     log = IterLog("VOMPS")
     iter = 0
     ϵ = calc_galerkin(mps, toapprox..., envs)
     alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-    recalculate!(envs, mps, toapprox...; alg_environments.tol)
+    recalculate!(envs, mps, toapprox..., alg_environments)
 
     state = VOMPSState(mps, toapprox, envs, iter, ϵ)
     it = IterativeSolver(alg, state)
@@ -63,9 +63,10 @@ function Base.iterate(it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:T
 end
 
 function localupdate_step!(
-        ::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler
+        it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler
     )
-    alg_orth = Defaults.alg_orth()
+    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
     dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs
@@ -86,9 +87,10 @@ function localupdate_step!(
     return ACs
 end
 function localupdate_step!(
-        ::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler
+        it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler
     )
-    alg_orth = Defaults.alg_orth()
+    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
     dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs
@@ -118,5 +120,5 @@ end
 
 function envs_step!(it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, mps)
     alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ)
-    return recalculate!(state.envs, mps, state.operator...; alg_environments.tol)
+    return recalculate!(state.envs, mps, state.operator..., alg_environments)
 end

src/algorithms/changebonds/optimalexpand.jl

@@ -22,7 +22,7 @@ end
 
 # Simple wrapper to convert between diffrent type of InifniteMPS.
 function changebonds(
-        ψ::InfiniteMPS, operator::InfiniteMPO, alg::OptimalExpand, envs = environments(ψ, operator)
+        ψ::InfiniteMPS, operator::InfiniteMPO, alg::OptimalExpand, envs = environments(ψ, operator, ψ)
     )
     ψ′, envs′ = changebonds(
         convert(MultilineMPS, ψ), convert(MultilineMPO, operator), alg, Multiline([envs])
@@ -32,7 +32,7 @@ end
 
 function changebonds(
         ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, alg::OptimalExpand,
-        envs = environments(ψ, H)
+        envs = environments(ψ, H, ψ)
     )
     T = eltype(ψ.AL)
     AL′ = similar(ψ.AL)
@@ -52,11 +52,11 @@ function changebonds(
     end
 
     newψ = _expand(ψ, AL′, AR′)
-    envs = environments(newψ, H)
+    envs = environments(newψ, H, newψ)
     return newψ, envs
 end
 
-function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environments(ψ, H))
+function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
     TL = eltype(ψ.AL)
     AL′ = PeriodicMatrix{TL}(undef, size(ψ.AL))
     TR = tensormaptype(spacetype(TL), 1, numind(TL) - 1, storagetype(TL))
@@ -77,15 +77,15 @@ function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environment
     end
 
     newψ = _expand(ψ, AL′, AR′)
-    envs = environments(newψ, H) #  recalculate!(envs, newψ, H)
+    envs = environments(newψ, H, newψ) #  recalculate!(envs, newψ, H)
     return newψ, envs
 end
 
-function changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H))
+function changebonds(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
     return changebonds!(copy(ψ), H, alg, envs)
 end
 
-function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H))
+function changebonds!(ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs = environments(ψ, H, ψ))
     #inspired by the infinite mps algorithm, alternative is to use https://arxiv.org/pdf/1501.05504.pdf
 
     #the idea is that we always want to expand the state in such a way that there are zeros at site i

src/algorithms/changebonds/randexpand.jl

@@ -56,7 +56,7 @@ changebonds(ψ::MultilineMPS, alg::RandExpand) = changebonds!(copy(ψ), alg)
 
 function changebonds(ψ, H, alg::RandExpand, envs)
     newψ = changebonds(ψ, alg)
-    return newψ, environments(newψ, H)
+    return newψ, environments(newψ, H, newψ)
 end
 
 

src/algorithms/changebonds/svdcut.jl

@@ -118,7 +118,7 @@ end
 
 function changebonds(ψ, H, alg::SvdCut, envs)
     newψ = changebonds(ψ, alg)
-    return newψ, environments(newψ, H)
+    return newψ, environments(newψ, H, newψ)
 end
 
 changebonds(mpo::FiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(mpo), alg)

src/algorithms/changebonds/vumpssvd.jl

@@ -25,7 +25,7 @@ $(TYPEDFIELDS)
 end
 
 function changebonds_1(
-        state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H)
+        state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H, state)
     ) # would be more efficient if we also repeated envs
     # the unitcell==1 case is unique, because there you have a sef-consistency condition
 
@@ -48,10 +48,10 @@ function changebonds_1(
         [nstate.AL[1]], nstate.C[1]; alg.alg_gauge.tol, alg.alg_gauge.maxiter
     )
 
-    return collapsed, environments(collapsed, H)
+    return collapsed, environments(collapsed, H, collapsed)
 end
 
-function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H))
+function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H, state))
     for loc in 1:length(state)
         @plansor AC2[-1 -2; -3 -4] := state.AC[loc][-1 -2; 1] * state.AR[loc + 1][1 -4; -3]
 
@@ -77,12 +77,12 @@ function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environme
         copied[loc] = AL1
         copied[loc + 1] = AL2
         state = InfiniteMPS(copied; alg.alg_gauge.tol, alg.alg_gauge.maxiter)
-        envs = environments(state, H)
+        envs = environments(state, H, state)
     end
     return state, envs
 end
 
-function changebonds(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H))
+function changebonds(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs = environments(state, H, state))
     return length(state) == 1 ? changebonds_1(state, H, alg, envs) :
         changebonds_n(state, H, alg, envs)
 end

src/algorithms/excitation/chepigaansatz.jl

@@ -34,7 +34,7 @@ function ChepigaAnsatz(; kwargs...)
 end
 
 function excitations(
-        H, alg::ChepigaAnsatz, ψ::FiniteMPS, envs = environments(ψ, H);
+        H, alg::ChepigaAnsatz, ψ::FiniteMPS, envs = environments(ψ, H, ψ);
         sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) || throw(ArgumentError("invalid position $pos"))
@@ -99,7 +99,7 @@ function ChepigaAnsatz2(; trscheme = notrunc(), kwargs...)
 end
 
 function excitations(
-        H, alg::ChepigaAnsatz2, ψ::FiniteMPS, envs = environments(ψ, H);
+        H, alg::ChepigaAnsatz2, ψ::FiniteMPS, envs = environments(ψ, H, ψ);
         sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) - 1 || throw(ArgumentError("invalid position $pos"))

src/algorithms/excitation/dmrgexcitation.jl

@@ -31,7 +31,7 @@ function excitations(
         tuple(H, ProjectionOperator.(states)...),
         tuple(1.0, broadcast(x -> alg.weight, states)...)
     )
-    envs = environments(init, super_op)
+    envs = environments(init, super_op, init)
     ne, _ = find_groundstate(init, super_op, alg.gsalg, envs)
 
     nstates = (states..., ne)