From ea90a4e64efbbf429d422937eebde480919c9074 Mon Sep 17 00:00:00 2001 From: leburgel Date: Mon, 15 Jun 2026 12:27:15 +0200 Subject: [PATCH 01/12] Pass through more settings for environment and gauging sub-algorithms --- src/algorithms/approximate/vomps.jl | 12 +++++++----- src/algorithms/groundstate/vumps.jl | 11 ++++++----- src/algorithms/statmech/vomps.jl | 13 +++++++------ src/utility/defaults.jl | 14 +++++++------- 4 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 8893858c6..72da14c31 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -18,7 +18,7 @@ function approximate( 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) @@ -65,7 +65,8 @@ end function localupdate_step!( ::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 diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 4a67f5f8a..8ab6f033f 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -64,7 +64,7 @@ function dominant_eigsolve( mps = copy(mps) ϵ = calc_galerkin(mps, operator, mps, envs) alg_environments = updatetol(alg.alg_environments, iter, ϵ) - recalculate!(envs, mps, operator, mps; alg_environments.tol) + recalculate!(envs, mps, operator, mps; alg_environments...) state = VUMPSState(mps, operator, envs, iter, ϵ, which, timeroutput) it = IterativeSolver(alg, state) @@ -118,8 +118,9 @@ end function localupdate_step!( it::IterativeSolver{<:VUMPS}, state, scheduler = Defaults.scheduler[] ) + alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) - alg_orth = Defaults.alg_orth() + alg_orth = alg_gauge.alg_orth mps = state.mps src_Cs = mps isa Multiline ? eachcol(mps.C) : mps.C @@ -188,20 +189,20 @@ function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractVector) alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) mps = gaugefix!( state.mps, ACs, state.mps.C[end]; - alg_gauge.tol, alg_gauge.maxiter, order = :R, timeroutput = state.timeroutput, + order = :R, timeroutput = state.timeroutput, alg_gauge..., ) mul!.(mps.AC, mps.AL, mps.C) return mps end function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractMatrix) alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) - return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge.tol, alg_gauge.maxiter) + return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge.tol, alg_gauge.maxiter, alg_gauge.alg_orth) end function envs_step!(it::IterativeSolver{<:VUMPS}, state, mps) alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) return recalculate!( state.envs, mps, state.operator, mps; - alg_environments.tol, state.timeroutput, + state.timeroutput, alg_environments..., ) end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 866d2d75b..f9f2be454 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -57,7 +57,7 @@ function dominant_eigsolve( iter = 0 ϵ = calc_galerkin(mps, operator, mps, envs) alg_environments = updatetol(alg.alg_environments, iter, ϵ) - recalculate!(envs, mps, operator, mps; alg_environments.tol) + recalculate!(envs, mps, operator, mps; alg_environments...) state = VOMPSState(mps, operator, envs, iter, ϵ) it = IterativeSolver(alg, state) @@ -100,9 +100,10 @@ function Base.iterate(it::IterativeSolver{<:VOMPS}, state) end function localupdate_step!( - ::IterativeSolver{<:VOMPS}, state, scheduler = Defaults.scheduler[] + it::IterativeSolver{<:VOMPS}, state, scheduler = Defaults.scheduler[] ) - alg_orth = Defaults.alg_orth() + alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) + alg_orth = alg_gauge.alg_orth mps = state.mps ACs = similar(mps.AC) dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs @@ -137,14 +138,14 @@ end function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractVector) alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) - return InfiniteMPS(ACs, state.mps.C[end]; alg_gauge.tol, alg_gauge.maxiter) + return InfiniteMPS(ACs, state.mps.C[end]; alg_gauge...) end function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractMatrix) alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) - return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge.tol, alg_gauge.maxiter) + return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge...) end function envs_step!(it::IterativeSolver{<:VOMPS}, state, mps) alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) - return recalculate!(state.envs, mps, state.operator, mps; alg_environments.tol) + return recalculate!(state.envs, mps, state.operator, mps; alg_environments...) end diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index b49fa1719..9ff7bb9e3 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -38,12 +38,16 @@ const eigsolver = Arnoldi(; tol, maxiter, eager = true) # Default algorithms # ------------------ +alg_svd() = DefaultAlgorithm() +alg_orth() = Householder(; positive = true) + function alg_gauge(; tol = tolgauge, maxiter = maxiter, verbosity = VERBOSE_WARN, + alg_orth = alg_orth(), dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max, tol_factor = gauge_tolfactor ) - alg = (; tol, maxiter, verbosity) + alg = (; tol, maxiter, verbosity, alg_orth) return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end @@ -57,16 +61,12 @@ function alg_eigsolve(; return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end -alg_svd() = DefaultAlgorithm() -alg_orth() = Householder(; positive = true) - -# TODO: make verbosity and maxiter actually do something function alg_environments(; - tol = tol, maxiter = maxiter, verbosity = 0, + tol = tol, maxiter = maxiter, verbosity = 0, krylovdim = krylovdim, eager = true, dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max, tol_factor = envs_tolfactor ) - alg = (; tol, maxiter, verbosity) + alg = (; tol, maxiter, verbosity, eager, krylovdim) return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end From 9f33665e588a47651e7635f4b47f55f7db61ef9d Mon Sep 17 00:00:00 2001 From: Lander Burgelman <39218680+leburgel@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:39:28 +0200 Subject: [PATCH 02/12] Update src/algorithms/approximate/vomps.jl Co-authored-by: Lukas Devos --- src/algorithms/approximate/vomps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 72da14c31..807fd188f 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -63,7 +63,7 @@ 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_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ) alg_orth = alg_gauge.alg_orth From 1cd618cd6e54d810a6995e00d81f90f62b7a69b0 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 16 Jun 2026 12:50:32 +0200 Subject: [PATCH 03/12] Fix kwarg mismatch between eigsolve and linsolve algorithms --- src/environments/abstract_envs.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index 612a37207..37d43750e 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -87,7 +87,7 @@ end function environment_alg( below, ::InfiniteMPOHamiltonian, above; tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim, - verbosity = Defaults.VERBOSE_NONE + verbosity = Defaults.VERBOSE_NONE, kwargs... ) max_krylovdim = ceil(Int, dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1))) return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity) @@ -96,7 +96,7 @@ function environment_alg( ::Union{InfiniteQP, MultilineQP}, ::Union{InfiniteMPO, MultilineMPO}, ::Union{InfiniteQP, MultilineQP}; tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim, - verbosity = Defaults.VERBOSE_NONE + verbosity = Defaults.VERBOSE_NONE, kwargs... ) return GMRES(; tol, maxiter, krylovdim, verbosity) end From ef845fc51121c82809cd84a2c8ea62f6e2d4b167 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 16 Jun 2026 16:23:51 -0400 Subject: [PATCH 04/12] pass around environments more deeply --- src/MPSKit.jl | 1 - src/algorithms/groundstate/vumps.jl | 13 +++--- src/algorithms/statmech/gradient_grassmann.jl | 2 +- src/algorithms/statmech/leading_boundary.jl | 13 ++++-- src/algorithms/statmech/vomps.jl | 12 +++-- src/algorithms/statmech/vumps.jl | 15 ------- src/environments/infinite_envs.jl | 45 +++++++++++++++---- src/environments/multiline_envs.jl | 30 ++++++++++++- src/environments/multiple_envs.jl | 26 ++++++++--- src/utility/defaults.jl | 2 +- src/utility/dynamictols.jl | 5 +++ 11 files changed, 112 insertions(+), 52 deletions(-) delete mode 100644 src/algorithms/statmech/vumps.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index f4d2fc6f8..c39977731 100644 --- a/src/MPSKit.jl +++ b/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") diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 8ab6f033f..aa1c413cb 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -47,13 +47,13 @@ struct VUMPSState{S, O, E} end function find_groundstate( - mps::InfiniteMPS, operator, alg::VUMPS, envs = environments(mps, operator) + mps::InfiniteMPS, operator, alg::VUMPS, envs... ) - return dominant_eigsolve(operator, mps, alg, envs; which = :SR) + return dominant_eigsolve(operator, mps, alg, envs...; which = :SR) end function dominant_eigsolve( - operator, mps, alg::VUMPS, envs = environments(mps, operator); + operator, mps, alg::VUMPS, envs = environments(mps, operator, mps, alg.alg_environments); which ) log = IterLog("VUMPS") @@ -64,7 +64,7 @@ function dominant_eigsolve( mps = copy(mps) ϵ = calc_galerkin(mps, operator, mps, envs) alg_environments = updatetol(alg.alg_environments, iter, ϵ) - recalculate!(envs, mps, operator, mps; alg_environments...) + recalculate!(envs, mps, operator, mps, alg_environments; timeroutput) state = VUMPSState(mps, operator, envs, iter, ϵ, which, timeroutput) it = IterativeSolver(alg, state) @@ -201,8 +201,5 @@ end function envs_step!(it::IterativeSolver{<:VUMPS}, state, mps) alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) - return recalculate!( - state.envs, mps, state.operator, mps; - state.timeroutput, alg_environments..., - ) + return recalculate!(state.envs, mps, state.operator, mps, alg_environments; state.timeroutput) end diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl index c69b2bb6b..b00a6349e 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -2,7 +2,7 @@ function leading_boundary( state::MultilineMPS, operator::MultilineMPO, alg::GradientGrassmann, - envs::MultilineEnvironments = environments(state, H) + envs::MultilineEnvironments = environments(state, operator) ) fg(x) = GrassmannMPS.fg(x, operator, envs) x, _, _, _, normgradhistory = optimize( diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index 943e7ce03..d66088eea 100644 --- a/src/algorithms/statmech/leading_boundary.jl +++ b/src/algorithms/statmech/leading_boundary.jl @@ -25,9 +25,16 @@ optimization algorithm will be attempted based on the supplied keywords. # TODO: alg selector # implementation always in terms of Multiline objects -function leading_boundary( - state::InfiniteMPS, operator::InfiniteMPO, alg, envs = environments(state, operator) +function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg) + state_multi = convert(MultilineMPS, state) + operator_multi = convert(MultilineMPO, operator) + state_multi′, envs_multi′, err = leading_boundary( + state_multi, operator_multi, alg, ) + state′ = convert(InfiniteMPS, state_multi′) + return state′, only(envs_multi′), err +end +function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg, envs) state_multi = convert(MultilineMPS, state) operator_multi = convert(MultilineMPO, operator) envs_multi = Multiline([envs]) @@ -36,5 +43,5 @@ function leading_boundary( envs_multi ) state′ = convert(InfiniteMPS, state_multi′) - return state′, envs, err + return state′, only(envs_multi′), err end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index f9f2be454..9295f2b28 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -42,14 +42,12 @@ struct VOMPSState{S, O, E} ϵ::Float64 end -function leading_boundary( - ψ::MultilineMPS, O::MultilineMPO, alg::VOMPS, envs = environments(ψ, O) - ) - return dominant_eigsolve(O, ψ, alg, envs; which = :LM) +function leading_boundary(ψ::MultilineMPS, O::MultilineMPO, alg::Union{VOMPS, VUMPS}, envs...) + return dominant_eigsolve(O, ψ, alg, envs...; which = :LM) end function dominant_eigsolve( - operator, mps, alg::VOMPS, envs = environments(mps, operator); + operator, mps, alg::VOMPS, envs = environments(mps, operator, mps, alg.alg_environments); which ) @assert which === :LM "VOMPS only supports the LM eigenvalue problem" @@ -57,7 +55,7 @@ function dominant_eigsolve( iter = 0 ϵ = calc_galerkin(mps, operator, mps, envs) alg_environments = updatetol(alg.alg_environments, iter, ϵ) - recalculate!(envs, mps, operator, mps; alg_environments...) + recalculate!(envs, mps, operator, mps, alg_environments) state = VOMPSState(mps, operator, envs, iter, ϵ) it = IterativeSolver(alg, state) @@ -147,5 +145,5 @@ end function envs_step!(it::IterativeSolver{<:VOMPS}, state, mps) alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ) - return recalculate!(state.envs, mps, state.operator, mps; alg_environments...) + return recalculate!(state.envs, mps, state.operator, mps, alg_environments) end diff --git a/src/algorithms/statmech/vumps.jl b/src/algorithms/statmech/vumps.jl deleted file mode 100644 index 83abac739..000000000 --- a/src/algorithms/statmech/vumps.jl +++ /dev/null @@ -1,15 +0,0 @@ -""" - leading_boundary(ψ, opp, alg, envs=environments(ψ, opp)) - -Approximate the leading eigenvector for opp. -""" -function leading_boundary(ψ::InfiniteMPS, H, alg, envs = environments(ψ, H)) - st, pr, de = leading_boundary(convert(MultilineMPS, ψ), Multiline([H]), alg, envs) - return convert(InfiniteMPS, st), pr, de -end - -function leading_boundary( - mps::MultilineMPS, operator, alg::VUMPS, envs = environments(mps, operator) - ) - return dominant_eigsolve(operator, mps, alg, envs; which = :LM) -end diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl index 1029bad52..c215977fd 100644 --- a/src/environments/infinite_envs.jl +++ b/src/environments/infinite_envs.jl @@ -19,12 +19,19 @@ leftenv(envs::InfiniteEnvironments, site::Int, state) = envs.GLs[site] rightenv(envs::InfiniteEnvironments, site::Int, state) = envs.GRs[site] function environments( - below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, - above::InfiniteMPS = below; kwargs... + below::InfiniteMPS, operator, above = below; + timeroutput::TimerOutput = DISABLED_TIMER, kwargs... + ) + alg = environment_alg(below, operator, above; kwargs...) + return environments(below, operator, above, alg; timeroutput) +end +function environments( + below::InfiniteMPS, operator, above, alg; + timeroutput::TimerOutput = DISABLED_TIMER ) GLs, GRs = initialize_environments(below, operator, above) envs = InfiniteEnvironments(GLs, GRs) - return recalculate!(envs, below, operator, above; kwargs...) + return recalculate!(envs, below, operator, above, alg; timeroutput) end function issamespace( @@ -47,13 +54,36 @@ function issamespace( return true end +function recalculate!( + envs::InfiniteEnvironments, below, operator, above = below; + timeroutput::TimerOutput = DISABLED_TIMER, kwargs... + ) + alg = environment_alg(below, operator, above; kwargs...) + return recalculate!(envs, below, operator, above, alg; timeroutput) +end function recalculate!( envs::InfiniteEnvironments, below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, - above::InfiniteMPS = below; - timeroutput::TimerOutput = DISABLED_TIMER, + above::InfiniteMPS, alg::DefaultAlgorithm; kwargs... ) + alg = environment_alg(below, operator, above; alg.kwargs...) + return recalculate!(envs, below, operator, above, alg; kwargs...) +end +function recalculate!( + envs::InfiniteEnvironments, below::InfiniteMPS, + operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, + above::InfiniteMPS, alg::DynamicTol; + kwargs... + ) + return recalculate!(envs, below, operator, above, alg.alg; kwargs...) +end +function recalculate!( + envs::InfiniteEnvironments, below::InfiniteMPS, + operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, + above::InfiniteMPS, alg; + timeroutput::TimerOutput = DISABLED_TIMER, + ) if !issamespace(envs, below, operator, above) # TODO: in-place initialization? GLs, GRs = initialize_environments(below, operator, above) @@ -61,8 +91,6 @@ function recalculate!( copy!(envs.GRs, GRs) end - alg = environment_alg(below, operator, above; kwargs...) - tree_point = String[section.name for section in timeroutput.timer_stack] @sync begin @spawn begin @@ -117,8 +145,7 @@ function compute_rightenvs!( # push through unitcell for i in reverse(1:(length(operator) - 1)) envs.GRs[i] = TransferMatrix( - above.AR[i + 1], operator[i + 1], - below.AR[i + 1] + above.AR[i + 1], operator[i + 1], below.AR[i + 1] ) * envs.GRs[i + 1] end return λ, envs diff --git a/src/environments/multiline_envs.jl b/src/environments/multiline_envs.jl index 128f64b55..f636510c4 100644 --- a/src/environments/multiline_envs.jl +++ b/src/environments/multiline_envs.jl @@ -11,10 +11,21 @@ function environments( end return Multiline(PeriodicVector(envs)) end +function environments( + below::MultilineMPS, operator::MultilineMPO, above::MultilineMPS, alg; + kwargs... + ) + (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || + throw(ArgumentError("Incompatible sizes")) + envs = map(1:rows) do row + return environments(below[row + 1], operator[row], above[row], alg; kwargs...) + end + return Multiline(PeriodicVector(envs)) +end function recalculate!( envs::MultilineEnvironments, below::MultilineMPS, - operator::MultilineMPO, above::MultilineMPS = below; + operator::MultilineMPO, above::MultilineMPS; kwargs... ) (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || @@ -24,11 +35,28 @@ function recalculate!( end return envs end +function recalculate!( + envs::MultilineEnvironments, below::MultilineMPS, + operator::MultilineMPO, above::MultilineMPS, alg; + kwargs... + ) + (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || + throw(ArgumentError("Incompatible sizes")) + @threads for row in 1:rows + recalculate!(envs[row], below[row + 1], operator[row], above[row], alg; kwargs...) + end + return envs +end function recalculate!( envs::MultilineEnvironments, below, (operator, above)::Tuple; kwargs... ) return recalculate!(envs, below, operator, above; kwargs...) end +function recalculate!( + envs::MultilineEnvironments, below, (operator, above)::Tuple, alg; kwargs... + ) + return recalculate!(envs, below, operator, above, alg; kwargs...) +end function TensorKit.normalize!(envs::MultilineEnvironments, below, operator, above) for row in 1:size(below, 1) diff --git a/src/environments/multiple_envs.jl b/src/environments/multiple_envs.jl index cc2c6460e..8373e58f2 100644 --- a/src/environments/multiple_envs.jl +++ b/src/environments/multiple_envs.jl @@ -10,8 +10,18 @@ Base.iterate(x::MultipleEnvironments) = iterate(x.envs) Base.iterate(x::MultipleEnvironments, i) = iterate(x.envs, i) # we need constructor, agnostic of particular MPS -function environments(state, H::LazySum; kwargs...) - return MultipleEnvironments(map(Base.Fix1(environments, state), H.ops)) +function environments(below, H::LazySum, above = below; kwargs...) + return MultipleEnvironments(map(x -> environments(below, x, above; kwargs...), H.ops)) +end +function environments(below, H::LazySum, above, alg; kwargs...) + return MultipleEnvironments(map(x -> environments(below, x, above, alg; kwargs...), H.ops)) +end +# disambiguate: +function environments(below::InfiniteMPS, H::LazySum, above = below; kwargs...) + return MultipleEnvironments(map(x -> environments(below, x, above; kwargs...), H.ops)) +end +function environments(below::InfiniteMPS, H::LazySum, above, alg; kwargs...) + return MultipleEnvironments(map(x -> environments(below, x, above, alg; kwargs...), H.ops)) end function environments( @@ -26,10 +36,6 @@ function environments( ) end -# we need to define how to recalculate -""" - Recalculate in-place each sub-env in MultipleEnvironments -""" function recalculate!( envs::MultipleEnvironments, below, operator::LazySum, above = below; kwargs... ) @@ -38,6 +44,14 @@ function recalculate!( end return envs end +function recalculate!( + envs::MultipleEnvironments, below, operator::LazySum, above, alg; kwargs... + ) + for (subenvs, subO) in zip(envs.envs, operator) + recalculate!(subenvs, below, subO, above, alg; kwargs...) + end + return envs +end #maybe this can be used to provide compatibility with existing code? function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl index 9ff7bb9e3..05aa58c69 100644 --- a/src/utility/defaults.jl +++ b/src/utility/defaults.jl @@ -66,7 +66,7 @@ function alg_environments(; dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max, tol_factor = envs_tolfactor ) - alg = (; tol, maxiter, verbosity, eager, krylovdim) + alg = DefaultAlgorithm(; tol, maxiter, verbosity, eager, krylovdim) return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg end diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl index 75bea4d2f..7a67a9a01 100644 --- a/src/utility/dynamictols.jl +++ b/src/utility/dynamictols.jl @@ -3,6 +3,7 @@ module DynamicTols import ..MPSKit: Algorithm using Accessors using DocStringExtensions +using MatrixAlgebraKit: DefaultAlgorithm export updatetol, DynamicTol @@ -71,5 +72,9 @@ end function _updatetol(alg, tol::Real) return Accessors.@set alg.tol = tol end +function _updatetol(alg::DefaultAlgorithm, tol::Real) + kwargs = merge(alg.kwargs, (; tol = tol)) + return DefaultAlgorithm(; kwargs...) +end end From b201212263fcd13970ea8fabc256a86253a6746d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 19 Jun 2026 15:52:28 -0400 Subject: [PATCH 05/12] refactor environments functions for consistency --- .../excitation/quasiparticleexcitation.jl | 18 +++-- src/algorithms/propagator/corvector.jl | 2 +- src/environments/abstract_envs.jl | 22 +++++ src/environments/finite_envs.jl | 80 +++++++++---------- src/environments/infinite_envs.jl | 47 ++++++----- src/environments/multiple_envs.jl | 7 -- src/environments/qp_envs.jl | 40 +++++----- 7 files changed, 117 insertions(+), 99 deletions(-) diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl index 691b096ec..0563277be 100644 --- a/src/algorithms/excitation/quasiparticleexcitation.jl +++ b/src/algorithms/excitation/quasiparticleexcitation.jl @@ -47,7 +47,7 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renv E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs) H_eff = EffectiveExcitationHamiltonian(H, lenvs, renvs, E) Es, ϕs, convhist = eigsolve(ϕ₀, num, :SR, alg.alg) do ϕ - return H_eff(ϕ; alg.alg_environments...) + return H_eff(ϕ, alg.alg_environments) end convhist.converged < num && @warn "excitation failed to converge: normres = $(convhist.normres)" @@ -149,7 +149,7 @@ function excitations( E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs) H_eff = EffectiveExcitationHamiltonian(H, lenvs, renvs, E) Es, ϕs, convhist = eigsolve(ϕ₀, num, :SR, alg.alg) do ϕ - return H_eff(ϕ; alg.alg_environments...) + return H_eff(ϕ, alg.alg_environments) end convhist.converged < num && @@ -201,7 +201,7 @@ function excitations( H_eff = Multiline(H_effs) Es, ϕs, convhist = eigsolve(ϕ₀, num, :LM, alg.alg) do ϕ - return H_eff(ϕ; alg.alg_environments...) + return H_eff(ϕ, alg.alg_environments) end convhist.converged < num && @warn "excitation failed to converge: normres = $(convhist.normres)" @@ -217,7 +217,7 @@ function excitations( H_eff = EffectiveExcitationHamiltonian(H_eff, lenvs, renvs, E) Es, ϕs, convhist = eigsolve(ϕ₀, num, :LM, alg.alg) do ϕ - return H_eff(ϕ; alg.alg_environments...) + return H_eff(ϕ, alg.alg_environments) end convhist.converged < num && @warn "excitation failed to converge: normres = $(convhist.normres)" @@ -288,12 +288,14 @@ end # to allow Multiline checks Base.length(H::EffectiveExcitationHamiltonian) = length(H.operator) -function (H::EffectiveExcitationHamiltonian)(ϕ::QP; kwargs...) - qp_envs = environments(ϕ, H.operator, H.lenvs, H.renvs; kwargs...) +function (H::EffectiveExcitationHamiltonian)(ϕ::QP, alg_environments = DefaultAlgorithm()) + qp_envs = environments(ϕ, H.operator, ϕ, alg_environments; lenvs = H.lenvs, renvs = H.renvs) return effective_excitation_hamiltonian(H.operator, ϕ, qp_envs, H.energy) end -function (H::Multiline{<:EffectiveExcitationHamiltonian})(ϕ::MultilineQP; kwargs...) - return Multiline(map((x, y) -> x(y; kwargs...), parent(H), parent(ϕ))) +function (H::Multiline{<:EffectiveExcitationHamiltonian})( + ϕ::MultilineQP, alg_environments = DefaultAlgorithm() + ) + return Multiline(map((x, y) -> x(y, alg_environments), parent(H), parent(ϕ))) end function effective_excitation_hamiltonian(H, ϕ, envs = environments(ϕ, H)) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index f4936c002..ad32dfdce 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -187,7 +187,7 @@ function squaredenvs( # to construct the squared caches we will first initialize environments # then make all data invalid so it will be recalculated - envs² = environments(state, H², leftstart, rightstart) + envs² = environments(state, H²; leftstart, rightstart) for i in 1:L poison!(envs², i) end diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index 37d43750e..bc3450712 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -5,6 +5,13 @@ Abstract supertype for all environment types. """ abstract type AbstractMPSEnvironments end +@doc """ + environments(below, operator, [above], [alg]; kwargs...) + +Construct the environments for the combination of `operator` sandwiched between the states `below` (bra) and `above` (ket). +The keyword arguments or `alg` struct can be used to control the settings needed for computing the environments. +""" environments + # Locking # ------- Base.lock(f::Function, envs::AbstractMPSEnvironments) = lock(f, envs.lock) @@ -100,3 +107,18 @@ function environment_alg( ) return GMRES(; tol, maxiter, krylovdim, verbosity) end + +""" + resolve_environment_solver(alg, below, operator, above) + +Resolve an environment algorithm `alg` into a concrete iterative solver for the given problem. +A `KrylovKit` algorithm is returned as-is; the `DefaultAlgorithm`/`DynamicTol` configuration +wrappers are translated via [`environment_alg`](@ref). +""" +resolve_environment_solver(alg, below, operator, above) = alg +function resolve_environment_solver(alg::DefaultAlgorithm, below, operator, above) + return environment_alg(below, operator, above; alg.kwargs...) +end +function resolve_environment_solver(alg::DynamicTol, below, operator, above) + return resolve_environment_solver(alg.alg, below, operator, above) +end diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index 6372f65ab..d762ca2c7 100644 --- a/src/environments/finite_envs.jl +++ b/src/environments/finite_envs.jl @@ -16,71 +16,65 @@ struct FiniteEnvironments{A, B, C, D} <: AbstractMPSEnvironments GRs::Vector{D} end -function environments(below, (operator, above)::Tuple, args...; kwargs...) - return environments(below, operator, above, args...; kwargs...) -end -function environments(below, operator, leftstart, rightstart) - return environments(below, operator, nothing, leftstart, rightstart) -end -function environments(below, operator, above, leftstart, rightstart) - leftenvs = [i == 0 ? leftstart : similar(leftstart) for i in 0:length(below)] +function initialize_environments(below::FiniteMPS, operator, above, leftstart, rightstart) N = length(below) - rightenvs = [i == N ? rightstart : similar(rightstart) for i in 0:length(below)] + leftenvs = [i == 0 ? leftstart : similar(leftstart) for i in 0:N] + rightenvs = [i == N ? rightstart : similar(rightstart) for i in 0:N] t = similar(below.AL[1]) return FiniteEnvironments( - above, operator, fill(t, length(below)), - fill(t, length(below)), - leftenvs, - rightenvs + above, operator, fill(t, N), fill(t, N), leftenvs, rightenvs ) end function environments( - below::FiniteMPS{S}, O::Union{FiniteMPO, FiniteMPOHamiltonian}, above = nothing + below::FiniteMPS{S}, operator::Union{FiniteMPO, FiniteMPOHamiltonian}, above = nothing; + leftstart = nothing, rightstart = nothing ) where {S} - Vl_bot = left_virtualspace(below, 1) - Vl_mid = left_virtualspace(O, 1) - Vl_top = isnothing(above) ? left_virtualspace(below, 1) : left_virtualspace(above, 1) - leftstart = isomorphism(storagetype(S), Vl_bot ⊗ Vl_mid' ← Vl_top) - N = length(below) - Vr_bot = right_virtualspace(below, N) - Vr_mid = right_virtualspace(O, N) - Vr_top = isnothing(above) ? right_virtualspace(below, N) : right_virtualspace(above, N) - rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid ← Vr_bot) - - return environments(below, O, above, leftstart, rightstart) + if isnothing(leftstart) + Vl_bot = left_virtualspace(below, 1) + Vl_mid = left_virtualspace(operator, 1) + Vl_top = isnothing(above) ? left_virtualspace(below, 1) : left_virtualspace(above, 1) + leftstart = isomorphism(storagetype(S), Vl_bot ⊗ Vl_mid' ← Vl_top) + end + if isnothing(rightstart) + Vr_bot = right_virtualspace(below, N) + Vr_mid = right_virtualspace(operator, N) + Vr_top = isnothing(above) ? right_virtualspace(below, N) : right_virtualspace(above, N) + rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid ← Vr_bot) + end + return initialize_environments(below, operator, above, leftstart, rightstart) end function environments( - below::WindowMPS, O::Union{InfiniteMPOHamiltonian, InfiniteMPO}, above = nothing; - lenvs = environments(below.left_gs, O), - renvs = environments(below.right_gs, O) + below::WindowMPS, operator::Union{InfiniteMPOHamiltonian, InfiniteMPO}, above = nothing; + lenvs = environments(below.left_gs, operator), + renvs = environments(below.right_gs, operator), + leftstart = copy(lenvs.GLs[1]), + rightstart = copy(renvs.GRs[end]) ) - leftstart = copy(lenvs.GLs[1]) - rightstart = copy(renvs.GRs[end]) - - return environments(below, O, above, leftstart, rightstart) + return initialize_environments(below, operator, above, leftstart, rightstart) end -function environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} +environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} = + environments(below, nothing, above) +function environments(below::S, ::Nothing, above::S) where {S <: Union{FiniteMPS, WindowMPS}} S isa WindowMPS && (above.left_gs == below.left_gs || throw(ArgumentError("left gs differs"))) S isa WindowMPS && (above.right_gs == below.right_gs || throw(ArgumentError("right gs differs"))) - operator = fill(nothing, length(below)) - return environments(below, operator, above, l_LL(above), r_RR(above)) + ops = fill(nothing, length(below)) + return initialize_environments(below, ops, above, l_LL(above), r_RR(above)) end -function environments(state::Union{FiniteMPS, WindowMPS}, operator::ProjectionOperator) - @plansor leftstart[-1; -2 -3 -4] := l_LL(operator.ket)[-3; -4] * - l_LL(operator.ket)[-1; -2] - @plansor rightstart[-1; -2 -3 -4] := r_RR(operator.ket)[-1; -2] * - r_RR(operator.ket)[-3; -4] - return environments( - state, fill(nothing, length(state)), operator.ket, leftstart, - rightstart +function environments( + state::Union{FiniteMPS, WindowMPS}, operator::ProjectionOperator, above = nothing + ) + @plansor leftstart[-1; -2 -3 -4] := l_LL(operator.ket)[-3; -4] * l_LL(operator.ket)[-1; -2] + @plansor rightstart[-1; -2 -3 -4] := r_RR(operator.ket)[-1; -2] * r_RR(operator.ket)[-3; -4] + return initialize_environments( + state, fill(nothing, length(state)), operator.ket, leftstart, rightstart ) end diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl index c215977fd..d82d427de 100644 --- a/src/environments/infinite_envs.jl +++ b/src/environments/infinite_envs.jl @@ -19,14 +19,15 @@ leftenv(envs::InfiniteEnvironments, site::Int, state) = envs.GLs[site] rightenv(envs::InfiniteEnvironments, site::Int, state) = envs.GRs[site] function environments( - below::InfiniteMPS, operator, above = below; + below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above = below; timeroutput::TimerOutput = DISABLED_TIMER, kwargs... ) alg = environment_alg(below, operator, above; kwargs...) return environments(below, operator, above, alg; timeroutput) end function environments( - below::InfiniteMPS, operator, above, alg; + below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above, + alg; timeroutput::TimerOutput = DISABLED_TIMER ) GLs, GRs = initialize_environments(below, operator, above) @@ -35,8 +36,8 @@ function environments( end function issamespace( - envs::InfiniteEnvironments, below::InfiniteMPS, - operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above::InfiniteMPS + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above::InfiniteMPS ) L = check_length(below, operator, above) for i in 1:L @@ -55,7 +56,8 @@ function issamespace( end function recalculate!( - envs::InfiniteEnvironments, below, operator, above = below; + envs::InfiniteEnvironments, below, + operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above = below; timeroutput::TimerOutput = DISABLED_TIMER, kwargs... ) alg = environment_alg(below, operator, above; kwargs...) @@ -121,8 +123,9 @@ function initialize_environments( end function compute_leftenvs!( - envs::InfiniteEnvironments, below::InfiniteMPS, - operator::InfiniteMPO, above::InfiniteMPS, alg + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::InfiniteMPO, above::InfiniteMPS, + alg ) # compute eigenvector T = TransferMatrix(above.AL, operator, below.AL) @@ -136,8 +139,9 @@ function compute_leftenvs!( end function compute_rightenvs!( - envs::InfiniteEnvironments, below::InfiniteMPS, operator::InfiniteMPO, - above::InfiniteMPS, alg + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::InfiniteMPO, above::InfiniteMPS, + alg ) # compute eigenvector T = TransferMatrix(above.AR, operator, below.AR) @@ -157,8 +161,8 @@ end # this avoids catastrophic blow-up of norms, while keeping the total normalized # and does not lead to issues for negative overlaps and real entries. function TensorKit.normalize!( - envs::InfiniteEnvironments, below::InfiniteMPS, operator::InfiniteMPO, - above::InfiniteMPS + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::InfiniteMPO, above::InfiniteMPS ) for i in 1:length(operator) normalize!(envs.GRs[i]) @@ -171,8 +175,7 @@ end # InfiniteMPOHamiltonian environments # ----------------------------------- function initialize_environments( - below::InfiniteMPS, operator::InfiniteMPOHamiltonian, - above::InfiniteMPS = below + below::InfiniteMPS, operator::InfiniteMPOHamiltonian, above::InfiniteMPS = below ) L = check_length(above, operator, below) GLs = PeriodicVector([allocate_GL(below, operator, above, i) for i in 1:L]) @@ -202,8 +205,9 @@ function initialize_environments( end function compute_leftenvs!( - envs::InfiniteEnvironments, below::InfiniteMPS, - operator::InfiniteMPOHamiltonian, above::InfiniteMPS, alg + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::InfiniteMPOHamiltonian, above::InfiniteMPS, + alg ) L = check_length(below, above, operator) GLs = envs.GLs @@ -257,8 +261,8 @@ function compute_leftenvs!( end function left_cyclethrough!( - index::Int, GL, below::InfiniteMPS, H::InfiniteMPOHamiltonian, - above::InfiniteMPS = below + index::Int, GL, + below::InfiniteMPS, H::InfiniteMPOHamiltonian, above::InfiniteMPS = below ) # TODO: efficient transfer matrix slicing for large unitcells leftinds = 1:index @@ -271,8 +275,9 @@ function left_cyclethrough!( end function compute_rightenvs!( - envs::InfiniteEnvironments, below::InfiniteMPS, - operator::InfiniteMPOHamiltonian, above::InfiniteMPS, alg + envs::InfiniteEnvironments, + below::InfiniteMPS, operator::InfiniteMPOHamiltonian, above::InfiniteMPS, + alg ) L = check_length(above, operator, below) GRs = envs.GRs @@ -327,8 +332,8 @@ function compute_rightenvs!( end function right_cyclethrough!( - index::Int, GR, below::InfiniteMPS, operator::InfiniteMPOHamiltonian, - above::InfiniteMPS = below + index::Int, GR, + below::InfiniteMPS, operator::InfiniteMPOHamiltonian, above::InfiniteMPS = below ) # TODO: efficient transfer matrix slicing for large unitcells for site in reverse(eachindex(GR)) diff --git a/src/environments/multiple_envs.jl b/src/environments/multiple_envs.jl index 8373e58f2..d750c8107 100644 --- a/src/environments/multiple_envs.jl +++ b/src/environments/multiple_envs.jl @@ -16,13 +16,6 @@ end function environments(below, H::LazySum, above, alg; kwargs...) return MultipleEnvironments(map(x -> environments(below, x, above, alg; kwargs...), H.ops)) end -# disambiguate: -function environments(below::InfiniteMPS, H::LazySum, above = below; kwargs...) - return MultipleEnvironments(map(x -> environments(below, x, above; kwargs...), H.ops)) -end -function environments(below::InfiniteMPS, H::LazySum, above, alg; kwargs...) - return MultipleEnvironments(map(x -> environments(below, x, above, alg; kwargs...), H.ops)) -end function environments( st::WindowMPS, H::LazySum; diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index 4aacf0dae..e8f3ae993 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -25,29 +25,32 @@ function rightenv(envs::InfiniteQPEnvironments, site::Int, state) return rightenv(envs.rightenvs, site, state) end -# Explicitly define optional arguments as these depend on kwargs, -# which needs to come after these arguments. - -function environments(exci::Union{InfiniteQP, MultilineQP}, H; kwargs...) - lenvs = environments(exci.left_gs, H; kwargs...) - return environments(exci, H, lenvs; kwargs...) -end -function environments(exci::Union{InfiniteQP, MultilineQP}, H, lenvs; kwargs...) - renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H; kwargs...) - return environments(exci, H, lenvs, renvs; kwargs...) +function environments( + exci::Union{InfiniteQP, MultilineQP}, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian, MultilineMPO}, above = exci; + lenvs = environments(exci.left_gs, operator), renvs = istopological(exci) ? environments(exci.right_gs, operator) : lenvs, + kwargs... + ) + alg = environment_alg(exci, operator, above; kwargs...) + return environments(exci, operator, above, alg; lenvs, renvs) end -function environments(qp::MultilineQP, operator::MultilineMPO, lenvs, renvs; kwargs...) +function environments( + qp::MultilineQP, operator::MultilineMPO, above, alg; lenvs, renvs + ) (rows = size(qp, 1)) == size(operator, 1) || throw(ArgumentError("Incompatible sizes")) envs = map(1:rows) do row - return environments(qp[row], operator[row], lenvs[row], renvs[row]; kwargs...) + return environments( + qp[row], operator[row], qp[row], alg; lenvs = lenvs[row], renvs = renvs[row] + ) end return Multiline(PeriodicVector(envs)) end -function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; kwargs...) +function environments( + exci::InfiniteQP, H::InfiniteMPOHamiltonian, above, alg; lenvs, renvs + ) ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1)) .+ 1 - solver = environment_alg(exci, H, exci; kwargs...) + solver = resolve_environment_solver(alg, exci, H, exci) AL = exci.left_gs.AL AR = exci.right_gs.AR @@ -130,10 +133,9 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; end function environments( - exci::FiniteQP, H::FiniteMPOHamiltonian, + exci::FiniteQP, H::FiniteMPOHamiltonian, above = exci, alg = nothing; lenvs = environments(exci.left_gs, H), - renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H); - kwargs... + renvs = istopological(exci) ? environments(exci.right_gs, H) : lenvs ) AL = exci.left_gs.AL AR = exci.right_gs.AR @@ -160,10 +162,10 @@ function environments( return InfiniteQPEnvironments(lBs, rBs, lenvs, renvs) end -function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...) +function environments(exci::InfiniteQP, O::InfiniteMPO, above, alg; lenvs, renvs) istopological(exci) && @warn "there is a phase ambiguity in topologically nontrivial statmech excitations" - solver = environment_alg(exci, O, exci; kwargs...) + solver = resolve_environment_solver(alg, exci, O, exci) left_gs = exci.left_gs right_gs = exci.right_gs From 7326dcc78afe2481ffdd5f5837584bb0169219dd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 19 Jun 2026 15:53:11 -0400 Subject: [PATCH 06/12] update changelog --- docs/src/changelog.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 95fdaec65..a4340407e 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -23,10 +23,19 @@ 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 = below, alg)`, where `alg` is the environment + algorithm (slot 4). Auxiliary inputs are keyword-only: `leftstart`/`rightstart` for finite + and window environments, and `lenvs`/`renvs` for window and quasiparticle environments. + ### Deprecated ### Removed +- The positional boundary forms `environments(below, operator, leftstart, rightstart)` and + `environments(below, operator, above, leftstart, rightstart)` for finite environments. Pass + the boundary tensors as the `leftstart`/`rightstart` keyword arguments instead. + ### Fixed ### Performance From a9bfd5fff9e9116018af95dc4008a93f9394440d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 19 Jun 2026 17:08:46 -0400 Subject: [PATCH 07/12] more small fixes --- src/algorithms/approximate/approximate.jl | 9 ++++----- src/algorithms/fidelity_susceptibility.jl | 2 +- src/environments/finite_envs.jl | 2 +- src/environments/qp_envs.jl | 4 ++-- 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 3c564cde2..8da83afae 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -29,7 +29,7 @@ approximate, approximate! # 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 +43,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( diff --git a/src/algorithms/fidelity_susceptibility.jl b/src/algorithms/fidelity_susceptibility.jl index 22ef69912..06be64b9f 100644 --- a/src/algorithms/fidelity_susceptibility.jl +++ b/src/algorithms/fidelity_susceptibility.jl @@ -28,7 +28,7 @@ function fidelity_susceptibility( end vec, convhist = linsolve(Tos, Tos, GMRES(; maxiter = maxiter, tol = tol)) do x - return effective_excitation_hamiltonian(H₀, x, environments(x, H₀, henvs)) + return effective_excitation_hamiltonian(H₀, x, environments(x, H₀, x; lenvs = henvs)) end convhist.converged == 0 && @warn "failed to converge: normres = $(convhist.normres)" diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index d762ca2c7..70733e447 100644 --- a/src/environments/finite_envs.jl +++ b/src/environments/finite_envs.jl @@ -16,7 +16,7 @@ struct FiniteEnvironments{A, B, C, D} <: AbstractMPSEnvironments GRs::Vector{D} end -function initialize_environments(below::FiniteMPS, operator, above, leftstart, rightstart) +function initialize_environments(below::AbstractFiniteMPS, operator, above, leftstart, rightstart) N = length(below) leftenvs = [i == 0 ? leftstart : similar(leftstart) for i in 0:N] rightenvs = [i == N ? rightstart : similar(rightstart) for i in 0:N] diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index e8f3ae993..75638b332 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -35,7 +35,7 @@ function environments( end function environments( - qp::MultilineQP, operator::MultilineMPO, above, alg; lenvs, renvs + qp::MultilineQP, operator::MultilineMPO, above, alg; lenvs, renvs = lenvs ) (rows = size(qp, 1)) == size(operator, 1) || throw(ArgumentError("Incompatible sizes")) envs = map(1:rows) do row @@ -47,7 +47,7 @@ function environments( end function environments( - exci::InfiniteQP, H::InfiniteMPOHamiltonian, above, alg; lenvs, renvs + exci::InfiniteQP, H::InfiniteMPOHamiltonian, above, alg; lenvs, renvs = lenvs ) ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1)) .+ 1 solver = resolve_environment_solver(alg, exci, H, exci) From dc0e6f1135352850540e6ba96828c439285b99c3 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 20 Jun 2026 10:50:28 -0400 Subject: [PATCH 08/12] more fixes --- src/algorithms/approximate/fvomps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 0ac30a047..dcfffd27e 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/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(ψ, 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(ψ, Oϕ...)) ϵ::Float64 = 2 * alg.tol log = IterLog("DMRG") From d0de0d1a2d4aebe77fdb9ba326d959cf01280bf6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 24 Jun 2026 20:17:34 -0400 Subject: [PATCH 09/12] double down on syntax --- src/algorithms/ED.jl | 2 +- src/algorithms/approximate/approximate.jl | 11 ++++++- src/algorithms/approximate/fvomps.jl | 4 +-- src/algorithms/approximate/idmrg.jl | 4 +-- src/algorithms/approximate/vomps.jl | 2 +- src/algorithms/changebonds/optimalexpand.jl | 14 ++++---- src/algorithms/changebonds/randexpand.jl | 2 +- src/algorithms/changebonds/svdcut.jl | 2 +- src/algorithms/changebonds/vumpssvd.jl | 10 +++--- src/algorithms/excitation/chepigaansatz.jl | 4 +-- src/algorithms/excitation/dmrgexcitation.jl | 2 +- .../excitation/quasiparticleexcitation.jl | 32 +++++++++---------- src/algorithms/expval.jl | 10 +++--- src/algorithms/fidelity_susceptibility.jl | 6 ++-- src/algorithms/grassmann.jl | 12 +++---- src/algorithms/groundstate/dmrg.jl | 4 +-- .../groundstate/find_groundstate.jl | 2 +- .../groundstate/gradient_grassmann.jl | 2 +- src/algorithms/groundstate/idmrg.jl | 2 +- src/algorithms/propagator/corvector.jl | 8 ++--- src/algorithms/statmech/gradient_grassmann.jl | 2 +- src/algorithms/statmech/idmrg.jl | 4 +-- src/algorithms/timestep/tdvp.jl | 10 +++--- src/algorithms/timestep/time_evolve.jl | 2 +- src/algorithms/toolbox.jl | 10 +++--- src/algorithms/unionalg.jl | 4 +-- src/environments/abstract_envs.jl | 13 ++++++-- src/environments/finite_envs.jl | 8 ++--- src/environments/infinite_envs.jl | 2 +- src/environments/lazylincocache.jl | 4 +-- src/environments/multiline_envs.jl | 2 +- src/environments/multiple_envs.jl | 4 +-- src/environments/qp_envs.jl | 6 ++-- 33 files changed, 112 insertions(+), 94 deletions(-) diff --git a/src/algorithms/ED.jl b/src/algorithms/ED.jl index df145b064..32c0e800a 100644 --- a/src/algorithms/ED.jl +++ b/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 diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 8da83afae..f28e67e6b 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/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,6 +30,11 @@ 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, diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index dcfffd27e..5a9ec4d29 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/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") diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 77de9bc65..f437a55fa 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/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 diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 807fd188f..1f5d15e11 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -12,7 +12,7 @@ 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 diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl index f6daa2025..194808704 100644 --- a/src/algorithms/changebonds/optimalexpand.jl +++ b/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 diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl index 67a245a1c..346755e4c 100644 --- a/src/algorithms/changebonds/randexpand.jl +++ b/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 diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 59eabe1da..a94dfb3c6 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/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) diff --git a/src/algorithms/changebonds/vumpssvd.jl b/src/algorithms/changebonds/vumpssvd.jl index 3c3dff2d0..a8e531e1b 100644 --- a/src/algorithms/changebonds/vumpssvd.jl +++ b/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 diff --git a/src/algorithms/excitation/chepigaansatz.jl b/src/algorithms/excitation/chepigaansatz.jl index e212861e0..2358ee4ab 100644 --- a/src/algorithms/excitation/chepigaansatz.jl +++ b/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")) diff --git a/src/algorithms/excitation/dmrgexcitation.jl b/src/algorithms/excitation/dmrgexcitation.jl index 57fdd8440..e3a7c79ce 100644 --- a/src/algorithms/excitation/dmrgexcitation.jl +++ b/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) diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl index 0563277be..dab59fecd 100644 --- a/src/algorithms/excitation/quasiparticleexcitation.jl +++ b/src/algorithms/excitation/quasiparticleexcitation.jl @@ -56,12 +56,12 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renv end function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs; num = 1, kwargs...) # Infer `renvs` in function body as it depends on `solver`. - renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H, ϕ₀.right_gs; kwargs...) return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...) end function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP; num = 1, kwargs...) # Infer `lenvs` in function body as it depends on `solver`. - lenvs = environments(ϕ₀.left_gs, H; kwargs...) + lenvs = environments(ϕ₀.left_gs, H, ϕ₀.left_gs; kwargs...) return excitations(H, alg, ϕ₀, lenvs; num, kwargs...) end @@ -90,8 +90,8 @@ Create and optimize infinite quasiparticle states. """ function excitations( H, alg::QuasiparticleAnsatz, momentum::Number, lmps::InfiniteMPS, - lenvs = environments(lmps, H), rmps::InfiniteMPS = lmps, - renvs = lmps === rmps ? lenvs : environments(rmps, H); + lenvs = environments(lmps, H, lmps), rmps::InfiniteMPS = lmps, + renvs = lmps === rmps ? lenvs : environments(rmps, H, rmps); sector = leftunit(lmps), kwargs... ) ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector, momentum) @@ -99,8 +99,8 @@ function excitations( end function excitations( H, alg::QuasiparticleAnsatz, momenta, lmps, - lenvs = environments(lmps, H), rmps = lmps, - renvs = lmps === rmps ? lenvs : environments(rmps, H); + lenvs = environments(lmps, H, lmps), rmps = lmps, + renvs = lmps === rmps ? lenvs : environments(rmps, H, rmps); verbosity = Defaults.verbosity, num = 1, sector = leftunit(lmps), parallel = true, kwargs... ) @@ -142,8 +142,8 @@ end function excitations( H, alg::QuasiparticleAnsatz, ϕ₀::FiniteQP, - lenvs = environments(ϕ₀.left_gs, H), - renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H); + lenvs = environments(ϕ₀.left_gs, H, ϕ₀.left_gs), + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H, ϕ₀.right_gs); num = 1 ) E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs) @@ -178,8 +178,8 @@ Create and optimize finite quasiparticle states. """ function excitations( H, alg::QuasiparticleAnsatz, lmps::FiniteMPS, - lenvs = environments(lmps, H), rmps::FiniteMPS = lmps, - renvs = lmps === rmps ? lenvs : environments(rmps, H); + lenvs = environments(lmps, H, lmps), rmps::FiniteMPS = lmps, + renvs = lmps === rmps ? lenvs : environments(rmps, H, rmps); sector = leftunit(lmps), num = 1 ) ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector) @@ -230,7 +230,7 @@ function excitations( kwargs... ) # Infer `renvs` in function body as it depends on `solver`. - renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H, ϕ₀.right_gs; kwargs...) return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...) end function excitations( @@ -238,14 +238,14 @@ function excitations( num = 1, kwargs... ) # Infer `lenvs` in function body as it depends on `solver`. - lenvs = environments(ϕ₀.left_gs, H; kwargs...) + lenvs = environments(ϕ₀.left_gs, H, ϕ₀.left_gs; kwargs...) return excitations(H, alg, ϕ₀, lenvs; num, kwargs...) end function excitations( H::MPO, alg::QuasiparticleAnsatz, momentum::Real, lmps::InfiniteMPS, - lenvs = environments(lmps, H), rmps::InfiniteMPS = lmps, - renvs = lmps === rmps ? lenvs : environments(rmps, H); + lenvs = environments(lmps, H, lmps), rmps::InfiniteMPS = lmps, + renvs = lmps === rmps ? lenvs : environments(rmps, H, rmps); kwargs... ) multiline_H = convert(MultilineMPO, H) @@ -267,8 +267,8 @@ end function excitations( H::MultilineMPO, alg::QuasiparticleAnsatz, momentum::Real, lmps::MultilineMPS, - lenvs = environments(lmps, H), rmps = lmps, - renvs = lmps === rmps ? lenvs : environments(rmps, H); + lenvs = environments(lmps, H, lmps), rmps = lmps, + renvs = lmps === rmps ? lenvs : environments(rmps, H, rmps); sector = leftunit(lmps), kwargs... ) ϕ₀ = LeftGaugedQP(randn, lmps, rmps; sector, momentum) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index b5db28f12..c147ba391 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -123,14 +123,14 @@ end function expectation_value( ψ::FiniteMPS, H::FiniteMPOHamiltonian, - envs::AbstractMPSEnvironments = environments(ψ, H) + envs::AbstractMPSEnvironments = environments(ψ, H, ψ) ) return dot(ψ, H, ψ, envs) / dot(ψ, ψ) end function expectation_value( ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, - envs::AbstractMPSEnvironments = environments(ψ, H) + envs::AbstractMPSEnvironments = environments(ψ, H, ψ) ) T = TensorOperations.promote_contract(scalartype(ψ), scalartype(H)) s = zero(T) @@ -146,7 +146,7 @@ end #----------- function expectation_value( ψ::InfiniteMPS, mpo_pair::Tuple{InfiniteMPO, Pair{Int, <:MPOTensor}}, - envs::AbstractMPSEnvironments = environments(ψ, first(mpo_pair)) + envs::AbstractMPSEnvironments = environments(ψ, first(mpo_pair), ψ) ) O_mpo, (site, O) = mpo_pair return contract_mpo_expval(ψ.AC[site], envs.GLs[site], O, envs.GRs[site]) / @@ -166,7 +166,7 @@ function expectation_value(ψ::InfiniteMPS, mpo::InfiniteMPO, envs...) end function expectation_value( ψ::MultilineMPS, O::MultilineMPO{<:InfiniteMPO}, - envs::MultilineEnvironments = environments(ψ, O) + envs::MultilineEnvironments = environments(ψ, O, ψ) ) return prod(product(1:size(ψ, 1), 1:size(ψ, 2))) do (i, j) GL = envs[i].GLs[j] @@ -208,7 +208,7 @@ end # ------------------ function expectation_value( ψ::FiniteMPS, O::ProjectionOperator, - envs::FiniteEnvironments = environments(ψ, O) + envs::FiniteEnvironments = environments(ψ, O, ψ) ) ens = zeros(scalartype(ψ), length(ψ)) for i in 1:length(ψ) diff --git a/src/algorithms/fidelity_susceptibility.jl b/src/algorithms/fidelity_susceptibility.jl index 06be64b9f..63f028401 100644 --- a/src/algorithms/fidelity_susceptibility.jl +++ b/src/algorithms/fidelity_susceptibility.jl @@ -1,6 +1,6 @@ """ fidelity_susceptibility(state::Union{FiniteMPS,InfiniteMPS}, H₀::T, - Vs::AbstractVector{T}, [henvs=environments(state, H₀)]; + Vs::AbstractVector{T}, [henvs=environments(state, H₀, state)]; maxiter=Defaults.maxiter, tol=Defaults.tol) where {T<:MPOHamiltonian} @@ -14,11 +14,11 @@ corresponding to each of the perturbing Hamiltonians. """ function fidelity_susceptibility( state::Union{FiniteMPS, InfiniteMPS}, H₀::T, - Vs::AbstractVector{T}, henvs = environments(state, H₀); + Vs::AbstractVector{T}, henvs = environments(state, H₀, state); maxiter = Defaults.maxiter, tol = Defaults.tol ) where {T <: MPOHamiltonian} tangent_vecs = map(Vs) do V - venvs = environments(state, V) + venvs = environments(state, V, state) Tos = LeftGaugedQP(rand, state) for (i, ac) in enumerate(state.AC) diff --git a/src/algorithms/grassmann.jl b/src/algorithms/grassmann.jl index d83f14b52..8cd8b02e6 100644 --- a/src/algorithms/grassmann.jl +++ b/src/algorithms/grassmann.jl @@ -137,13 +137,13 @@ function transport!(h, state, g, α::Real, state′) end """ - fg(state, operator, envs=environments(state, operator)) + fg(state, operator, envs=environments(state, operator, state)) Compute the cost function and the tangent vector with respect to the `AL` parameters of the state. """ function fg( state::FiniteMPS, operator::Union{O, LazySum{O}}, - envs::AbstractMPSEnvironments = environments(state, operator); + envs::AbstractMPSEnvironments = environments(state, operator, state); timeroutput::TimerOutput = DISABLED_TIMER, ) where {O <: FiniteMPOHamiltonian} f = @timeit timeroutput "expval" expectation_value(state, operator, envs) @@ -157,7 +157,7 @@ function fg( end function fg( state::InfiniteMPS, operator::Union{O, LazySum{O}}, - envs::AbstractMPSEnvironments = environments(state, operator); + envs::AbstractMPSEnvironments = environments(state, operator, state); timeroutput::TimerOutput = DISABLED_TIMER, ) where {O <: InfiniteMPOHamiltonian} @timeit timeroutput "envs (parallel)" recalculate!(envs, state, operator, state; timeroutput) @@ -177,7 +177,7 @@ function fg( end function fg( state::InfiniteMPS, operator::Union{O, LazySum{O}}, - envs::AbstractMPSEnvironments = environments(state, operator); + envs::AbstractMPSEnvironments = environments(state, operator, state); timeroutput::TimerOutput = DISABLED_TIMER, ) where {O <: InfiniteMPO} @timeit timeroutput "envs (parallel)" recalculate!(envs, state, operator, state; timeroutput) @@ -197,7 +197,7 @@ function fg( end function fg( state::MultilineMPS, operator::MultilineMPO, - envs::MultilineEnvironments = environments(state, operator); + envs::MultilineEnvironments = environments(state, operator, state); timeroutput::TimerOutput = DISABLED_TIMER, ) @assert length(state) == 1 "not implemented" @@ -245,7 +245,7 @@ end # utility test function function optimtest( - ψ, O, envs = environments(ψ, O); + ψ, O, envs = environments(ψ, O, ψ); alpha = -0.1:0.001:0.1, retract = retract, inner = inner ) _fg(x) = fg(x, O, envs) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 055b90beb..ca8b09d18 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -32,7 +32,7 @@ function DMRG(; return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize) end -function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H)) +function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ)) ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) ϵ = maximum(ϵs) log = IterLog("DMRG") @@ -120,7 +120,7 @@ function DMRG2(; return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_svd, trscheme, finalize) end -function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environments(ψ, H)) +function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environments(ψ, H, ψ)) ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ)) ϵ = maximum(ϵs) log = IterLog("DMRG2") diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 375d25f39..28ebc182d 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -22,7 +22,7 @@ optimization algorithm will be attempted based on the supplied keywords. - `ϵ::Float64`: final convergence error upon terminating the algorithm """ function find_groundstate( - ψ::AbstractMPS, H, envs::AbstractMPSEnvironments = environments(ψ, H); + ψ::AbstractMPS, H, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, trscheme = nothing ) diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl index 714d6b432..049cfaae5 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -56,7 +56,7 @@ struct GradientGrassmann{O <: OptimKit.OptimizationAlgorithm, F} <: Algorithm end function find_groundstate( - ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H) + ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H, ψ) )::Tuple{S, P, Float64} where {S, P} !isa(ψ, FiniteMPS) || dim(ψ.C[end]) == 1 || @warn "This is not fully supported - split the mps up in a sum of mps's and optimize separately" diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 1d95873eb..ed83e664a 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -74,7 +74,7 @@ function IDMRGState{T}( return IDMRGState{S, O, E, T}(mps, operator, envs, iter, ϵ, T(energy), timeroutput) end -function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, operator)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} +function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, operator, mps)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= 2")) name = alg isa IDMRG ? "IDMRG" : "IDMRG2" log = IterLog(name) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index ad32dfdce..de04e8078 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -62,7 +62,7 @@ function propagator( A::AbstractFiniteMPS, z::Number, H::FiniteMPOHamiltonian, alg::DynamicalDMRG{NaiveInvert}; init = copy(A) ) - h_envs = environments(init, H) # environments for h + h_envs = environments(init, H, init) # environments for h mixedenvs = environments(init, A) # environments for ϵ = 2 * alg.tol @@ -125,7 +125,7 @@ function propagator( ω = real(z) η = imag(z) - envs1 = environments(init, H) # environments for h + envs1 = environments(init, H, init) # environments for h H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2 mixedenvs = environments(init, A) # environments for @@ -176,7 +176,7 @@ function propagator( end function squaredenvs( - state::AbstractFiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H) + state::AbstractFiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H, state) ) H² = conj(H) * H L = length(state) @@ -187,7 +187,7 @@ function squaredenvs( # to construct the squared caches we will first initialize environments # then make all data invalid so it will be recalculated - envs² = environments(state, H²; leftstart, rightstart) + envs² = environments(state, H², state; leftstart, rightstart) for i in 1:L poison!(envs², i) end diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl index b00a6349e..20bef7b38 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -2,7 +2,7 @@ function leading_boundary( state::MultilineMPS, operator::MultilineMPO, alg::GradientGrassmann, - envs::MultilineEnvironments = environments(state, operator) + envs::MultilineEnvironments = environments(state, operator, state) ) fg(x) = GrassmannMPS.fg(x, operator, envs) x, _, _, _, normgradhistory = optimize( diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 477583bb6..614e9f6bf 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -1,5 +1,5 @@ function leading_boundary( - ψ::MultilineMPS, operator, alg::IDMRG, envs = environments(ψ, operator) + ψ::MultilineMPS, operator, alg::IDMRG, envs = environments(ψ, operator, ψ) ) log = IterLog("IDMRG") ϵ::Float64 = 2 * alg.tol @@ -62,7 +62,7 @@ function leading_boundary( end function leading_boundary( - ψ::MultilineMPS, operator, alg::IDMRG2, envs = environments(ψ, operator) + ψ::MultilineMPS, operator, alg::IDMRG2, envs = environments(ψ, operator, ψ) ) size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2")) ϵ::Float64 = 2 * alg.tol diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 2d67f90a1..a6fc56c2f 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -27,7 +27,7 @@ end function timestep( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::AbstractMPSEnvironments = environments(ψ, H); + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); leftorthflag = true, imaginary_evolution::Bool = false ) # convert state to complex if necessary @@ -86,7 +86,7 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::AbstractMPSEnvironments = environments(ψ, H); + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); imaginary_evolution::Bool = false ) @@ -166,7 +166,7 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs::AbstractMPSEnvironments = environments(ψ, H); + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); imaginary_evolution::Bool = false ) @@ -221,12 +221,12 @@ function timestep( ψ′ = isreal ? complex(ψ) : copy(ψ) if length(envs) != 0 && isreal @warn "Currently cannot reuse real environments for complex evolution" - envs′ = environments(ψ′, H) + envs′ = environments(ψ′, H, ψ′) elseif length(envs) == 1 envs′ = only(envs) else @assert length(envs) == 0 "Invalid signature" - envs′ = environments(ψ′, H) + envs′ = environments(ψ′, H, ψ′) end return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, kwargs...) end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index a0918ad9f..51288abd7 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -26,7 +26,7 @@ function time_evolve end, function time_evolve! end for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!)) @eval function $time_evolve( ψ, H, t_span::AbstractVector{<:Number}, alg, - envs = environments(ψ, H); + envs = environments(ψ, H, ψ); verbosity::Int = 0, imaginary_evolution::Bool = false ) log = IterLog("TDVP") diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 3d9cad9dd..a319baede 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -185,14 +185,14 @@ function correlation_length(spectrum::AbstractVector{T}; kwargs...) where {T <: end """ - variance(state, hamiltonian, [envs=environments(state, hamiltonian)]) + variance(state, hamiltonian, [envs=environments(state, hamiltonian, state)]) Compute the variance of the energy of the state with respect to the Hamiltonian. """ function variance end function variance( - state::InfiniteMPS, H::InfiniteMPOHamiltonian, envs = environments(state, H) + state::InfiniteMPS, H::InfiniteMPOHamiltonian, envs = environments(state, H, state) ) e_local = map(1:length(state)) do i return contract_mpo_expval( @@ -206,7 +206,7 @@ function variance( return real(expectation_value(state, (H - H_renormalized)^2)) end -function variance(state::FiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H)) +function variance(state::FiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H, state)) H2 = H * H return real(expectation_value(state, H2) - expectation_value(state, H, envs)^2) end @@ -235,7 +235,7 @@ function variance(state::InfiniteQP, H::InfiniteMPOHamiltonian, envs = environme # TODO: this is probably broken E_ex = dot(state, effective_excitation_hamiltonian(H, state, envs)) - rescaled_envs = environments(gs, H_regularized) + rescaled_envs = environments(gs, H_regularized, gs) GL = leftenv(rescaled_envs, 1, gs) GR = rightenv(rescaled_envs, 0, gs) E_f = @plansor GL[5 3; 1] * gs.C[0][1; 4] * conj(gs.C[0][5; 2]) * GR[4 3; 2] @@ -247,7 +247,7 @@ function variance(state::InfiniteQP, H::InfiniteMPOHamiltonian, envs = environme ) end -function variance(ψ, H::LazySum, envs = environments(ψ, sum(H))) +function variance(ψ, H::LazySum, envs = environments(ψ, sum(H), ψ)) # TODO: avoid throwing error and just compute correct environments envs isa MultipleEnvironments && throw(ArgumentError("The environment cannot be Lazy i.e. use environments of sum(H)")) diff --git a/src/algorithms/unionalg.jl b/src/algorithms/unionalg.jl index 201f24d91..3e37f67f4 100644 --- a/src/algorithms/unionalg.jl +++ b/src/algorithms/unionalg.jl @@ -22,13 +22,13 @@ function changebonds(state, alg::UnionAlg) return state end -function changebonds(state, H, alg::UnionAlg, envs = environments(state, H)) +function changebonds(state, H, alg::UnionAlg, envs = environments(state, H, state)) state, envs = changebonds(state, H, alg.alg1, envs) state, envs = changebonds(state, H, alg.alg2, envs) return state, envs end -function find_groundstate(state, H, alg::UnionAlg, envs = environments(state, H)) +function find_groundstate(state, H, alg::UnionAlg, envs = environments(state, H, state)) state, envs = find_groundstate(state, H, alg.alg1, envs) return find_groundstate(state, H, alg.alg2, envs) end diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index bc3450712..b527412b6 100644 --- a/src/environments/abstract_envs.jl +++ b/src/environments/abstract_envs.jl @@ -6,10 +6,19 @@ Abstract supertype for all environment types. abstract type AbstractMPSEnvironments end @doc """ - environments(below, operator, [above], [alg]; kwargs...) + environments(below, operator, above, [alg]; kwargs...) + environments(below, above) -Construct the environments for the combination of `operator` sandwiched between the states `below` (bra) and `above` (ket). +Construct the environments for the `operator` sandwiched between the states `below` (bra) and `above` (ket). The keyword arguments or `alg` struct can be used to control the settings needed for computing the environments. + +The two-argument form `environments(below, above)` (with two states) constructs the overlap environments, i.e. the operator-free environments between the bra `below` and the ket `above`. + +!!! note + The operator form requires an explicit `above`; there is no two-argument + `environments(below, operator)` shorthand. Because a two-argument call is the overlap form, + the second argument cannot be disambiguated between a ket (overlap) and an operator — this is + genuinely undecidable for density matrices, where states and operators share a representation. """ environments # Locking diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index 70733e447..3dee085a9 100644 --- a/src/environments/finite_envs.jl +++ b/src/environments/finite_envs.jl @@ -28,7 +28,7 @@ function initialize_environments(below::AbstractFiniteMPS, operator, above, left end function environments( - below::FiniteMPS{S}, operator::Union{FiniteMPO, FiniteMPOHamiltonian}, above = nothing; + below::FiniteMPS{S}, operator::Union{FiniteMPO, FiniteMPOHamiltonian}, above; leftstart = nothing, rightstart = nothing ) where {S} N = length(below) @@ -47,9 +47,9 @@ function environments( return initialize_environments(below, operator, above, leftstart, rightstart) end function environments( - below::WindowMPS, operator::Union{InfiniteMPOHamiltonian, InfiniteMPO}, above = nothing; - lenvs = environments(below.left_gs, operator), - renvs = environments(below.right_gs, operator), + below::WindowMPS, operator::Union{InfiniteMPOHamiltonian, InfiniteMPO}, above; + lenvs = environments(below.left_gs, operator, below.left_gs), + renvs = environments(below.right_gs, operator, below.right_gs), leftstart = copy(lenvs.GLs[1]), rightstart = copy(renvs.GRs[end]) ) diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl index d82d427de..25285068e 100644 --- a/src/environments/infinite_envs.jl +++ b/src/environments/infinite_envs.jl @@ -19,7 +19,7 @@ leftenv(envs::InfiniteEnvironments, site::Int, state) = envs.GLs[site] rightenv(envs::InfiniteEnvironments, site::Int, state) = envs.GRs[site] function environments( - below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above = below; + below::InfiniteMPS, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian}, above; timeroutput::TimerOutput = DISABLED_TIMER, kwargs... ) alg = environment_alg(below, operator, above; kwargs...) diff --git a/src/environments/lazylincocache.jl b/src/environments/lazylincocache.jl index d4fafd913..8c9ada7f6 100644 --- a/src/environments/lazylincocache.jl +++ b/src/environments/lazylincocache.jl @@ -3,6 +3,6 @@ struct LazyLincoCache{A <: LinearCombination, C <: Tuple} <: AbstractMPSEnvironm envs::C end -function environments(state, H::LinearCombination) - return LazyLincoCache(H, broadcast(o -> environments(state, o), H.opps)) +function environments(below, H::LinearCombination, above = below) + return LazyLincoCache(H, broadcast(o -> environments(below, o, above), H.opps)) end diff --git a/src/environments/multiline_envs.jl b/src/environments/multiline_envs.jl index f636510c4..a8e1d73f4 100644 --- a/src/environments/multiline_envs.jl +++ b/src/environments/multiline_envs.jl @@ -1,7 +1,7 @@ const MultilineEnvironments{E <: AbstractMPSEnvironments} = Multiline{E} function environments( - below::MultilineMPS, operator::MultilineMPO, above::MultilineMPS = below; + below::MultilineMPS, operator::MultilineMPO, above::MultilineMPS; kwargs... ) (rows = size(above, 1)) == size(operator, 1) == size(below, 1) || diff --git a/src/environments/multiple_envs.jl b/src/environments/multiple_envs.jl index d750c8107..9935d52d5 100644 --- a/src/environments/multiple_envs.jl +++ b/src/environments/multiple_envs.jl @@ -18,12 +18,12 @@ function environments(below, H::LazySum, above, alg; kwargs...) end function environments( - st::WindowMPS, H::LazySum; + st::WindowMPS, H::LazySum, above = st; lenvs = environments(st.left_gs, H), renvs = environments(st.right_gs, H) ) return MultipleEnvironments( map( - (op, sublenv, subrenv) -> environments(st, op; lenvs = sublenv, renvs = subrenv), + (op, sublenv, subrenv) -> environments(st, op, above; lenvs = sublenv, renvs = subrenv), H.ops, lenvs, renvs ) ) diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index 75638b332..8ade2287c 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -27,7 +27,7 @@ end function environments( exci::Union{InfiniteQP, MultilineQP}, operator::Union{InfiniteMPO, InfiniteMPOHamiltonian, MultilineMPO}, above = exci; - lenvs = environments(exci.left_gs, operator), renvs = istopological(exci) ? environments(exci.right_gs, operator) : lenvs, + lenvs = environments(exci.left_gs, operator, exci.left_gs), renvs = istopological(exci) ? environments(exci.right_gs, operator, exci.right_gs) : lenvs, kwargs... ) alg = environment_alg(exci, operator, above; kwargs...) @@ -134,8 +134,8 @@ end function environments( exci::FiniteQP, H::FiniteMPOHamiltonian, above = exci, alg = nothing; - lenvs = environments(exci.left_gs, H), - renvs = istopological(exci) ? environments(exci.right_gs, H) : lenvs + lenvs = environments(exci.left_gs, H, exci.left_gs), + renvs = istopological(exci) ? environments(exci.right_gs, H, exci.right_gs) : lenvs ) AL = exci.left_gs.AL AR = exci.right_gs.AR From 7c80369ea7c6b3e3000168ecef70c1e5f456a50b Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 24 Jun 2026 20:19:01 -0400 Subject: [PATCH 10/12] update changelog --- docs/src/changelog.md | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index a4340407e..7e5dfbf04 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -24,18 +24,20 @@ 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 = below, alg)`, where `alg` is the environment - algorithm (slot 4). Auxiliary inputs are keyword-only: `leftstart`/`rightstart` for finite - and window environments, and `lenvs`/`renvs` for window and quasiparticle environments. + `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 -- The positional boundary forms `environments(below, operator, leftstart, rightstart)` and - `environments(below, operator, above, leftstart, rightstart)` for finite environments. Pass - the boundary tensors as the `leftstart`/`rightstart` keyword arguments instead. - ### Fixed ### Performance From 7991f1f0f19a642ee3283ef46b0904c8eafc5d8a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 24 Jun 2026 20:24:47 -0400 Subject: [PATCH 11/12] more fixes --- src/algorithms/approximate/vomps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 1f5d15e11..455265256 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -18,7 +18,7 @@ function approximate( iter = 0 ϵ = calc_galerkin(mps, toapprox..., envs) alg_environments = updatetol(alg.alg_environments, iter, ϵ) - recalculate!(envs, mps, toapprox...; alg_environments...) + recalculate!(envs, mps, toapprox..., alg_environments) state = VOMPSState(mps, toapprox, envs, iter, ϵ) it = IterativeSolver(alg, state) @@ -120,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...) + return recalculate!(state.envs, mps, state.operator..., alg_environments) end From a55cfc5ba8414cd985dc39300779aaed59fd454d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 25 Jun 2026 16:06:12 -0400 Subject: [PATCH 12/12] more fixes? --- src/environments/finite_envs.jl | 1 + test/operators/lazysum.jl | 16 ++++++++-------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index 3dee085a9..fee1555bb 100644 --- a/src/environments/finite_envs.jl +++ b/src/environments/finite_envs.jl @@ -17,6 +17,7 @@ struct FiniteEnvironments{A, B, C, D} <: AbstractMPSEnvironments end function initialize_environments(below::AbstractFiniteMPS, operator, above, leftstart, rightstart) + above = above === below ? nothing : above N = length(below) leftenvs = [i == 0 ? leftstart : similar(leftstart) for i in 0:N] rightenvs = [i == N ? rightstart : similar(rightstart) for i in 0:N] diff --git a/test/operators/lazysum.jl b/test/operators/lazysum.jl index 7da831e6d..1e0785210 100644 --- a/test/operators/lazysum.jl +++ b/test/operators/lazysum.jl @@ -96,8 +96,8 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 => 10, 3 => 5, 5 => 1)) end summedH = LazySum(Hs) - envs = map(H -> environments(ψ, H), Hs) - summed_envs = environments(ψ, summedH) + envs = map(H -> environments(ψ, H, ψ), Hs) + summed_envs = environments(ψ, summedH, ψ) expval = sum(zip(Hs, envs)) do (H, env) return expectation_value(ψ, H, env) @@ -134,8 +134,8 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 => 10, 3 => 5, 5 => 1)) t = 1.1 summedH_at = summedH(t) - envs = map(H -> environments(ψ, H), Hs) - summed_envs = environments(ψ, summedH) + 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) @@ -181,8 +181,8 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 => 10, 3 => 5, 5 => 1)) return repeat(H, 2) end summedH = LazySum(Hs) - envs = map(H -> environments(ψ, H), Hs) - summed_envs = environments(ψ, summedH) + envs = map(H -> environments(ψ, H, ψ), Hs) + summed_envs = environments(ψ, summedH, ψ) expval = sum(zip(Hs, envs)) do (H, Env) return expectation_value(ψ, H, Env) @@ -219,8 +219,8 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 => 10, 3 => 5, 5 => 1)) t = 1.1 summedH_at = summedH(t) - envs = map(H -> environments(ψ, H), Hs) - summed_envs = environments(ψ, summedH) + 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)