diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 338b43b44..129dc2e0c 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -106,10 +106,9 @@ include("transfermatrix/transfer.jl")
 
 include("environments/abstract_envs.jl")
 include("environments/finite_envs.jl")
-include("environments/infinitempo_envs.jl")
-include("environments/infinitempohamiltonian_envs.jl")
+include("environments/infinite_envs.jl")
+include("environments/multiline_envs.jl")
 include("environments/qp_envs.jl")
-include("environments/idmrg_envs.jl")
 include("environments/multiple_envs.jl")
 include("environments/lazylincocache.jl")
 
@@ -145,6 +144,7 @@ include("algorithms/excitation/dmrgexcitation.jl")
 include("algorithms/excitation/chepigaansatz.jl")
 include("algorithms/excitation/exci_transfer_system.jl")
 
+include("algorithms/statmech/leading_boundary.jl")
 include("algorithms/statmech/vumps.jl")
 include("algorithms/statmech/vomps.jl")
 include("algorithms/statmech/gradient_grassmann.jl")
diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl
index e107cb3ee..2ce627d07 100644
--- a/src/algorithms/approximate/approximate.jl
+++ b/src/algorithms/approximate/approximate.jl
@@ -25,3 +25,35 @@ of an MPS `ψ₀`.
 - `VOMPS`: Tangent space method for truncating uniform MPS.
 """
 approximate, approximate!
+
+# implementation in terms of Multiline
+function approximate(ψ::InfiniteMPS,
+                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS},
+                     algorithm,
+                     envs=environments(ψ, toapprox))
+    envs′ = Multiline([envs])
+    multi, envs = approximate(convert(MultilineMPS, ψ),
+                              (convert(MultilineMPO, toapprox[1]),
+                               convert(MultilineMPS, toapprox[2])), algorithm, envs′)
+    ψ = convert(InfiniteMPS, multi)
+    return ψ, envs
+end
+
+# dispatch to in-place method
+function approximate(ψ, toapprox, alg::Union{DMRG,DMRG2,IDMRG1,IDMRG2},
+                     envs=environments(ψ, toapprox))
+    return approximate!(copy(ψ), toapprox, alg, envs)
+end
+
+# disambiguate
+function approximate(ψ::InfiniteMPS,
+                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS},
+                     algorithm::Union{IDMRG1,IDMRG2},
+                     envs=environments(ψ, toapprox))
+    envs′ = Multiline([envs])
+    multi, envs = approximate(convert(MultilineMPS, ψ),
+                              (convert(MultilineMPO, toapprox[1]),
+                               convert(MultilineMPS, toapprox[2])), algorithm, envs′)
+    ψ = convert(InfiniteMPS, multi)
+    return ψ, envs
+end
diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index 7b85e9f54..7333721d4 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -1,15 +1,5 @@
-# dispatch to in-place method
-function approximate(ψ, toapprox, alg::Union{DMRG,DMRG2}, envs...)
-    return approximate!(copy(ψ), toapprox, alg, envs...)
-end
-
-function approximate!(ψ::AbstractFiniteMPS, sq, alg, envs=environments(ψ, sq))
-    tor = approximate!(ψ, [sq], alg, [envs])
-    return (tor[1], tor[2][1], tor[3])
-end
-
-function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
-                      envs=[environments(ψ, sq) for sq in squash])
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2,
+                      envs=environments(ψ, Oϕ))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG2")
 
@@ -18,9 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
         for iter in 1:(alg.maxiter)
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = sum(zip(squash, envs)) do (sq, pr)
-                    return ac2_proj(pos, ψ, pr)
-                end
+                AC2′ = ac2_proj(pos, ψ, Oϕ, envs)
                 al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme)
 
                 AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
@@ -31,7 +19,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
             end
 
             # finalize
-            ψ, envs = alg.finalize(iter, ψ, squash, envs)::Tuple{typeof(ψ),typeof(envs)}
+            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ),typeof(envs)}
 
             if ϵ < alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
@@ -48,8 +36,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG2,
     return ψ, envs, ϵ
 end
 
-function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG,
-                      envs=[environments(ψ, sq) for sq in squash])
+function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs=environments(ψ, Oϕ))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("DMRG")
 
@@ -58,10 +45,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG,
         for iter in 1:(alg.maxiter)
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = sum(zip(squash, envs)) do (sq, pr)
-                    return ac_proj(pos, ψ, pr)
-                end
-
+                AC′ = ac_proj(pos, ψ, Oϕ, envs)
                 AC = ψ.AC[pos]
                 ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
 
@@ -69,7 +53,7 @@ function approximate!(ψ::AbstractFiniteMPS, squash::Vector, alg::DMRG,
             end
 
             # finalize
-            ψ, envs = alg.finalize(iter, ψ, squash, envs)::Tuple{typeof(ψ),typeof(envs)}
+            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ),typeof(envs)}
 
             if ϵ < alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 6384be00d..04b5e3bba 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -1,8 +1,5 @@
-function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
-                     alg::IDMRG1, oenvs=environments(ost, toapprox))
-    ψ = copy(ost)
-    mpo, above = toapprox
-    envs = IDMRGEnvironments(ost, oenvs)
+function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
+                      alg::IDMRG1, envs=environments(ψ, toapprox))
     log = IterLog("IDMRG")
     ϵ::Float64 = 2 * alg.tol
 
@@ -12,31 +9,27 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             C_current = ψ.C[:, 0]
 
             # left to right sweep
-            for col in 1:size(ψ, 2), row in 1:size(ψ, 1)
-                h = MPO_∂∂AC(mpo[row, col], leftenv(envs, row, col),
-                             rightenv(envs, row, col))
-                ψ.AC[row + 1, col] = h * above.AC[row, col]
-                normalize!(ψ.AC[row + 1, col])
-
-                ψ.AL[row + 1, col], ψ.C[row + 1, col] = leftorth(ψ.AC[row + 1, col])
-
-                tm = TransferMatrix(above.AL[row, col], mpo[row, col], ψ.AL[row + 1, col])
-                setleftenv!(envs, row, col + 1, normalize(leftenv(envs, row, col) * tm))
+            for col in 1:size(ψ, 2)
+                for row in 1:size(ψ, 1)
+                    ψ.AC[row + 1, col] = ac_proj(row, col, ψ, toapprox, envs)
+                    normalize!(ψ.AC[row + 1, col])
+                    ψ.AL[row + 1, col], ψ.C[row + 1, col] = leftorth!(ψ.AC[row + 1, col])
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
             end
 
             # right to left sweep
-            for col in size(ψ, 2):-1:1, row in 1:size(ψ, 1)
-                h = MPO_∂∂AC(mpo[row, col], leftenv(envs, row, col),
-                             rightenv(envs, row, col))
-                ψ.AC[row + 1, col] = h * above.AC[row, col]
-                normalize!(ψ.AC[row + 1, col])
-
-                ψ.C[row + 1, col - 1], temp = rightorth(_transpose_tail(ψ.AC[row + 1, col]))
-                ψ.AR[row + 1, col] = _transpose_front(temp)
-
-                tm = TransferMatrix(above.AR[row, col], mpo[row, col], ψ.AR[row + 1, col])
-                setrightenv!(envs, row, col - 1, normalize(tm * rightenv(envs, row, col)))
+            for col in size(ψ, 2):-1:1
+                for row in 1:size(ψ, 1)
+                    ψ.AC[row + 1, col] = ac_proj(row, col, ψ, toapprox, envs)
+                    normalize!(ψ.AC[row + 1, col])
+                    ψ.C[row + 1, col - 1], temp = rightorth!(_transpose_tail(ψ.AC[row + 1,
+                                                                                  col]))
+                    ψ.AR[row + 1, col] = _transpose_front(temp)
+                end
+                transfer_rightenv!(envs, ψ, toapprox, col - 1)
             end
+            normalize!(envs, ψ, toapprox)
 
             ϵ = norm(C_current - ψ.C[:, 0])
 
@@ -52,19 +45,20 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
         end
     end
 
+    # TODO: immediately compute in-place
     nst = MultilineMPS(map(x -> x, ψ.AR); tol=alg.tol_gauge)
-    nenvs = environments(nst, toapprox)
-    return nst, nenvs, ϵ
+    copy!(ψ, nst) # ensure output destination is unchanged
+
+    recalculate!(envs, ψ, toapprox)
+    return ψ, envs, ϵ
 end
 
-function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
-                     alg::IDMRG2, oenvs=environments(ost, toapprox))
-    length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    mpo, above = toapprox
-    ψ = copy(ost)
-    envs = IDMRGEnvironments(ost, oenvs)
+function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
+                      alg::IDMRG2, envs=environments(ψ, toapprox))
+    size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2"))
     ϵ::Float64 = 2 * alg.tol
     log = IterLog("IDMRG2")
+    O, ϕ = toapprox
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ)
@@ -72,52 +66,41 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             C_current = ψ.C[:, 0]
 
             # sweep from left to right
-            for col in 1:size(ψ, 2), row in 1:size(ψ, 1)
-                ac2 = above.AC[row, col] * _transpose_tail(above.AR[row, col + 1])
-                h = MPO_∂∂AC2(mpo[row, col], mpo[row, col + 1], leftenv(envs, row, col),
-                              rightenv(envs, row, col + 1))
-
-                al, c, ar, = tsvd!(h * ac2; trunc=alg.trscheme, alg=TensorKit.SVD())
-                normalize!(c)
-
-                ψ.AL[row + 1, col] = al
-                ψ.C[row + 1, col] = complex(c)
-                ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-
-                setleftenv!(envs, row, col + 1,
-                            normalize(leftenv(envs, row, col) *
-                                      TransferMatrix(above.AL[row, col], mpo[row, col],
-                                                     ψ.AL[row + 1, col])))
-                setrightenv!(envs, row, col,
-                             normalize(TransferMatrix(above.AR[row, col + 1],
-                                                      mpo[row, col + 1],
-                                                      ψ.AR[row + 1, col + 1]) *
-                                       rightenv(envs, row, col + 1)))
+            for col in 1:size(ψ, 2)
+                for row in 1:size(ψ, 1)
+                    AC2′ = ac2_proj(row, col, ψ, toapprox, envs)
+                    al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=TensorKit.SVD())
+                    normalize!(c)
+
+                    ψ.AL[row + 1, col] = al
+                    ψ.C[row + 1, col] = complex(c)
+                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                    ψ.AC[row + 1, col + 1] = _transpose_front(c * ar)
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
+                transfer_rightenv!(envs, ψ, toapprox, col)
             end
+            normalize!(envs, ψ, toapprox)
 
             # sweep from right to left
-            for col in (size(ψ, 2) - 1):-1:0, row in 1:size(ψ, 1)
-                ac2 = above.AL[row, col] * _transpose_tail(above.AC[row, col + 1])
-                h = MPO_∂∂AC2(mpo[row, col], mpo[row, col + 1], leftenv(envs, row, col),
-                              rightenv(envs, row, col + 1))
-
-                al, c, ar, = tsvd!(h * ac2; trunc=alg.trscheme, alg=TensorKit.SVD())
-                normalize!(c)
-
-                ψ.AL[row + 1, col] = al
-                ψ.C[row + 1, col] = complex(c)
-                ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-
-                setleftenv!(envs, row, col + 1,
-                            normalize(leftenv(envs, row, col) *
-                                      TransferMatrix(above.AL[row, col], mpo[row, col],
-                                                     ψ.AL[row + 1, col])))
-                setrightenv!(envs, row, col,
-                             normalize(TransferMatrix(above.AR[row, col + 1],
-                                                      mpo[row, col + 1],
-                                                      ψ.AR[row + 1, col + 1]) *
-                                       rightenv(envs, row, col + 1)))
+            for col in (size(ψ, 2) - 1):-1:0
+                for row in 1:size(ψ, 1)
+                    # TODO: also write this as ac2_proj?
+                    AC2 = ϕ.AL[row, col] * _transpose_tail(ϕ.AC[row, col + 1])
+                    AC2′ = ∂AC2(AC2, O[row, col], O[row, col + 1],
+                                leftenv(envs[row], col, ψ[row]),
+                                rightenv(envs[row], col, ψ[row]))
+                    al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=TensorKit.SVD())
+                    normalize!(c)
+
+                    ψ.AL[row + 1, col] = al
+                    ψ.C[row + 1, col] = complex(c)
+                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                end
+                transfer_leftenv!(envs, ψ, toapprox, col + 1)
+                transfer_rightenv!(envs, ψ, toapprox, col)
             end
+            normalize!(envs, ψ, toapprox)
 
             # update error
             ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
@@ -139,7 +122,10 @@ function approximate(ost::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
         end
     end
 
+    # TODO: immediately compute in-place
     nst = MultilineMPS(map(x -> x, ψ.AR); tol=alg.tol_gauge)
-    nenvs = environments(nst, toapprox)
-    return nst, nenvs, ϵ
+    copy!(ψ, nst) # ensure output destination is unchanged
+    recalculate!(envs, ψ, toapprox)
+
+    return ψ, envs, ϵ
 end
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index d8d6f5269..f1b2fd448 100644
--- a/src/algorithms/approximate/vomps.jl
+++ b/src/algorithms/approximate/vomps.jl
@@ -1,14 +1,3 @@
-function approximate(ψ::InfiniteMPS,
-                     toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS}, algorithm,
-                     envs=environments(ψ, toapprox))
-    # PeriodicMPO's always act on MultilineMPS's. To avoid code duplication, define everything in terms of MultilineMPS's.
-    multi, envs = approximate(convert(MultilineMPS, ψ),
-                              (convert(MultilineMPO, toapprox[1]),
-                               convert(MultilineMPS, toapprox[2])), algorithm, envs)
-    ψ = convert(InfiniteMPS, multi)
-    return ψ, envs
-end
-
 Base.@deprecate(approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
                             alg::VUMPS, envs...; kwargs...),
                 approximate(ψ, toapprox,
@@ -18,10 +7,12 @@ Base.@deprecate(approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:M
 
 function approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:MultilineMPS},
                      alg::VOMPS, envs=environments(ψ, toapprox))
-    ϵ::Float64 = calc_galerkin(ψ, envs)
+    ϵ::Float64 = calc_galerkin(ψ, toapprox..., envs)
     temp_ACs = similar.(ψ.AC)
     scheduler = Defaults.scheduler[]
     log = IterLog("VOMPS")
+    alg_environments = updatetol(alg.alg_environments, 0, ϵ)
+    recalculate!(envs, ψ, toapprox...; alg_environments.tol)
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ)
@@ -34,11 +25,11 @@ function approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multilin
             ψ = MultilineMPS(temp_ACs, ψ.C[:, end]; alg_gauge.tol, alg_gauge.maxiter)
 
             alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-            recalculate!(envs, ψ; alg_environments.tol)
+            recalculate!(envs, ψ, toapprox...; alg_environments.tol)
 
             ψ, envs = alg.finalize(iter, ψ, toapprox, envs)::Tuple{typeof(ψ),typeof(envs)}
 
-            ϵ = calc_galerkin(ψ, envs)
+            ϵ = calc_galerkin(ψ, toapprox..., envs)
 
             if ϵ <= alg.tol
                 @infov 2 logfinish!(log, iter, ϵ)
@@ -55,18 +46,20 @@ function approximate(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multilin
     return ψ, envs, ϵ
 end
 
-function _vomps_localupdate(loc, ψ, (O, ψ₀), envs, factalg=QRpos())
+function _vomps_localupdate(loc, ψ, Oϕ, envs, factalg=QRpos())
     local tmp_AC, tmp_C
     if Defaults.scheduler[] isa SerialScheduler
-        tmp_AC = circshift([ac_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
-        tmp_C = circshift([c_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+        tmp_AC = circshift([ac_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)], 1)
+        tmp_C = circshift([c_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)], 1)
     else
         @sync begin
             Threads.@spawn begin
-                tmp_AC = circshift([ac_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+                tmp_AC = circshift([ac_proj(row, loc, ψ, Oϕ, envs)
+                                    for row in 1:size(ψ, 1)], 1)
             end
             Threads.@spawn begin
-                tmp_C = circshift([c_proj(row, loc, ψ, envs) for row in 1:size(ψ, 1)], 1)
+                tmp_C = circshift([c_proj(row, loc, ψ, Oϕ, envs) for row in 1:size(ψ, 1)],
+                                  1)
             end
         end
     end
diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl
index 009351328..f2cab1b87 100644
--- a/src/algorithms/changebonds/changebonds.jl
+++ b/src/algorithms/changebonds/changebonds.jl
@@ -9,6 +9,14 @@ See also: [`SvdCut`](@ref), [`RandExpand`](@ref), [`VUMPSSvdCut`](@ref), [`Optim
 function changebonds end
 function changebonds! end
 
+# write in terms of MultilineMPS
+function changebonds(ψ::InfiniteMPS, operator::InfiniteMPO, alg,
+                     envs=environments(ψ, operator))
+    ψ′, envs′ = changebonds(convert(MultilineMPS, ψ), convert(MultilineMPO, operator), alg,
+                            Multiline([envs]))
+    return convert(InfiniteMPS, ψ′), envs
+end
+
 _expand(ψ, AL′, AR′) = _expand!(copy(ψ), AL′, AR′)
 function _expand!(ψ::InfiniteMPS, AL′::PeriodicVector, AR′::PeriodicVector)
     for i in 1:length(ψ)
diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 492e30f4c..4b542eb69 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -13,7 +13,8 @@ original ψ.
     trscheme::TruncationScheme = truncdim(1)
 end
 
-function changebonds(ψ::InfiniteMPS, H, alg::OptimalExpand, envs=environments(ψ, H))
+function changebonds(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian, alg::OptimalExpand,
+                     envs=environments(ψ, H))
     T = eltype(ψ.AL)
     AL′ = similar(ψ.AL)
     AR′ = similar(ψ.AR, tensormaptype(spacetype(T), 1, numind(T) - 1, storagetype(T)))
@@ -36,12 +37,6 @@ function changebonds(ψ::InfiniteMPS, H, alg::OptimalExpand, envs=environments(
     return newψ, envs
 end
 
-function changebonds(ψ::InfiniteMPS, H::DenseMPO, alg::OptimalExpand,
-                     envs=environments(ψ, H))
-    (nmψ, envs) = changebonds(convert(MultilineMPS, ψ), convert(MultilineMPO, H), alg, envs)
-    return (convert(InfiniteMPS, nmψ), envs)
-end
-
 function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs=environments(ψ, H))
     TL = eltype(ψ.AL)
     AL′ = PeriodicMatrix{TL}(undef, size(ψ.AL))
diff --git a/src/algorithms/changebonds/vumpssvd.jl b/src/algorithms/changebonds/vumpssvd.jl
index 3380ac603..e79dfaf25 100644
--- a/src/algorithms/changebonds/vumpssvd.jl
+++ b/src/algorithms/changebonds/vumpssvd.jl
@@ -26,8 +26,8 @@ function changebonds_1(state::InfiniteMPS, H, alg::VUMPSSvdCut,
 
     nstate, nenvs = changebonds(nstate, nH, alg)
 
-    D1 = space(nstate.AL[1], 1)
-    D2 = space(nstate.AL[2], 1)
+    D1 = left_virtualspace(nstate, 1)
+    D2 = left_virtualspace(nstate, 2)
 
     # collapse back to 1 site
     if D2 != D1
@@ -36,6 +36,7 @@ function changebonds_1(state::InfiniteMPS, H, alg::VUMPSSvdCut,
     end
 
     collapsed = InfiniteMPS([nstate.AL[1]], nstate.C[1]; tol=alg.tol_gauge)
+    recalculate!(envs, collapsed, H, collapsed)
 
     return collapsed, envs
 end
@@ -72,8 +73,8 @@ function changebonds_n(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs=environment
         copied[loc] = AL1
         copied[loc + 1] = AL2
         state = InfiniteMPS(copied; tol=alg.tol_gauge)
+        recalculate!(envs, state, H, state)
     end
-
     return state, envs
 end
 
diff --git a/src/algorithms/derivatives.jl b/src/algorithms/derivatives.jl
index ad562a799..d9351b14e 100644
--- a/src/algorithms/derivatives.jl
+++ b/src/algorithms/derivatives.jl
@@ -31,42 +31,36 @@ Base.:*(h::Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}, v) = h(v);
 (h::DerivativeOperator)(v, ::Number) = h(v)
 
 # draft operator constructors
-function ∂∂C(pos::Int, mps, operator::AbstractMPO, cache)
-    return MPO_∂∂C(leftenv(cache, pos + 1, mps), rightenv(cache, pos, mps))
-end
-function ∂∂C(col::Int, mps, operator::MultilineMPO, envs)
-    return MPO_∂∂C(leftenv(envs, col + 1, mps), rightenv(envs, col, mps))
+function ∂∂C(pos::Int, mps, operator, envs)
+    return MPO_∂∂C(leftenv(envs, pos + 1, mps), rightenv(envs, pos, mps))
 end
 function ∂∂C(row::Int, col::Int, mps, operator::MultilineMPO, envs)
-    return MPO_∂∂C(leftenv(envs, row, col + 1, mps), rightenv(envs, row, col, mps))
+    return ∂∂C(col, mps[row], operator[row], envs[row])
 end
 
-function ∂∂AC(pos::Int, mps, operator::AbstractMPO, cache)
-    return MPO_∂∂AC(operator[pos], leftenv(cache, pos, mps), rightenv(cache, pos, mps))
+function ∂∂AC(pos::Int, mps, operator, envs)
+    return MPO_∂∂AC(operator[pos], leftenv(envs, pos, mps), rightenv(envs, pos, mps))
 end
 function ∂∂AC(row::Int, col::Int, mps, operator::MultilineMPO, envs)
-    return MPO_∂∂AC(operator[row, col], leftenv(envs, row, col, mps),
-                    rightenv(envs, row, col, mps))
+    return ∂∂AC(col, mps[row], operator[row], envs[row])
 end
 function ∂∂AC(col::Int, mps, operator::MultilineMPO, envs)
-    return MPO_∂∂AC(envs.operator[:, col], leftenv(envs, col, mps),
-                    rightenv(envs, col, mps))
-end;
-
-function ∂∂AC2(pos::Int, mps, operator::AbstractMPO, cache)
-    return MPO_∂∂AC2(operator[pos], operator[pos + 1], leftenv(cache, pos, mps),
-                     rightenv(cache, pos + 1, mps))
-end;
+    return MPO_∂∂AC(operator[:, col], leftenv(envs, col, mps), rightenv(envs, col, mps))
+end
+
+function ∂∂AC2(pos::Int, mps, operator, envs)
+    return MPO_∂∂AC2(operator[pos], operator[pos + 1], leftenv(envs, pos, mps),
+                     rightenv(envs, pos + 1, mps))
+end
 function ∂∂AC2(col::Int, mps, operator::MultilineMPO, envs)
     return MPO_∂∂AC2(operator[:, col], operator[:, col + 1], leftenv(envs, col, mps),
                      rightenv(envs, col + 1, mps))
 end
-function ∂∂AC2(row::Int, col::Int, mps, operator::MultilineMPO, envs)
-    return MPO_∂∂AC2(operator[row, col], operator[row, col + 1],
-                     leftenv(envs, row, col, mps), rightenv(envs, row, col + 1, mps))
+function ∂∂AC2(row::Int, col::Int, mps, operator::MultilineMPO, envs::MultilineEnvironments)
+    return ∂∂AC2(col, mps[row], operator[row], envs[row])
 end
 
-#allow calling them with CartesianIndices
+# allow calling them with CartesianIndices
 ∂∂C(pos::CartesianIndex, mps, operator, envs) = ∂∂C(Tuple(pos)..., mps, operator, envs)
 ∂∂AC(pos::CartesianIndex, mps, operator, envs) = ∂∂AC(Tuple(pos)..., mps, operator, envs)
 ∂∂AC2(pos::CartesianIndex, mps, operator, envs) = ∂∂AC2(Tuple(pos)..., mps, operator, envs)
@@ -93,8 +87,10 @@ function ∂AC(x::MPSTensor, ::Nothing, leftenv, rightenv)
 end
 
 function ∂AC2(x::MPOTensor, operator1::MPOTensor, operator2::MPOTensor, leftenv, rightenv)
-    @plansor toret[-1 -2; -3 -4] := leftenv[-1 7; 6] * x[6 5; 1 3] * operator1[7 -2; 5 4] *
-                                    operator2[4 -4; 3 2] * rightenv[1 2; -3]
+    @plansor toret[-1 -2; -3 -4] := leftenv[-1 7; 6] * x[6 5; 1 3] *
+                                    operator1[7 -2; 5 4] *
+                                    operator2[4 -4; 3 2] *
+                                    rightenv[1 2; -3]
     return toret isa BlockTensorMap ? only(toret) : toret
 end
 function ∂AC2(x::MPOTensor, ::Nothing, ::Nothing, leftenv, rightenv)
@@ -118,41 +114,47 @@ function ∂C(x::Vector, leftenv, rightenv)
     return circshift(map(t -> ∂C(t...), zip(x, leftenv, rightenv)), 1)
 end
 
-#downproject for approximate
-function c_proj(pos, below, envs::FiniteEnvironments)
-    return ∂C(envs.above.C[pos], leftenv(envs, pos + 1, below), rightenv(envs, pos, below))
+# downproject for approximate
+function c_proj(pos::Int, ψ, (operator, ϕ)::Tuple, envs)
+    return ∂C(ϕ.C[pos], leftenv(envs, pos + 1, ψ), rightenv(envs, pos, ψ))
 end
-
-function c_proj(row, col, below, envs::InfiniteMPOEnvironments)
-    return ∂C(envs.above.C[row, col], leftenv(envs, row, col + 1, below),
-              rightenv(envs, row, col, below))
+function c_proj(pos::Int, ψ, ϕ::AbstractMPS, envs)
+    return ∂C(ϕ.C[pos], leftenv(envs, pos + 1, ψ), rightenv(envs, pos, ψ))
+end
+function c_proj(pos::Int, ψ, Oϕs::LazySum, envs)
+    return sum(zip(Oϕs.ops, envs.envs)) do x
+        return c_proj(pos, ψ, x...)
+    end
+end
+function c_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs)
+    return c_proj(col, ψ[row], (O[row], ϕ[row]), envs[row])
 end
 
-# TODO: rewrite this to not use operator from cache?
-function ac_proj(pos, below, envs)
-    le = leftenv(envs, pos, below)
-    re = rightenv(envs, pos, below)
-
-    return ∂AC(envs.above.AC[pos], envs.operator[pos], le, re)
+function ac_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs)
+    return ∂AC(ϕ.AC[pos], O[pos], leftenv(envs, pos, ψ), rightenv(envs, pos, ψ))
+end
+function ac_proj(pos::Int, ψ, ϕ::AbstractMPS, envs)
+    return ∂AC(ϕ.AC[pos], nothing, leftenv(envs, pos, ψ), rightenv(envs, pos, ψ))
 end
-function ac_proj(row, col, below, envs::InfiniteMPOEnvironments)
-    return ∂AC(envs.above.AC[row, col], envs.operator[row, col],
-               leftenv(envs, row, col, below),
-               rightenv(envs, row, col, below))
+function ac_proj(pos::Int, ψ, Oϕs::LazySum, envs)
+    return sum(zip(Oϕs.ops, envs.envs)) do x
+        return ac_proj(pos, ψ, x...)
+    end
+end
+function ac_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs)
+    return ac_proj(col, ψ[row], (O[row], ϕ[row]), envs[row])
 end
-function ac2_proj(pos, below, envs)
-    le = leftenv(envs, pos, below)
-    re = rightenv(envs, pos + 1, below)
 
-    return ∂AC2(envs.above.AC[pos] * _transpose_tail(envs.above.AR[pos + 1]),
-                envs.operator[pos],
-                envs.operator[pos + 1], le, re)
+function ac2_proj(pos::Int, ψ, (O, ϕ)::Tuple, envs)
+    AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1])
+    return ∂AC2(AC2, O[pos], O[pos + 1], leftenv(envs, pos, ψ), rightenv(envs, pos + 1, ψ))
+end
+function ac2_proj(pos::Int, ψ, ϕ::AbstractMPS, envs)
+    AC2 = ϕ.AC[pos] * _transpose_tail(ϕ.AR[pos + 1])
+    return ∂AC2(AC2, nothing, nothing, leftenv(envs, pos, ψ), rightenv(envs, pos + 1, ψ))
 end
-function ac2_proj(row, col, below, envs::InfiniteMPOEnvironments)
-    @plansor ac2[-1 -2; -3 -4] := envs.above.AC[row, col][-1 -2; 1] *
-                                  envs.above.AR[row, col + 1][1 -4; -3]
-    return ∂AC2(ac2, leftenv(envs, row, col + 1, below),
-                rightenv(envs, row, col + 1, below))
+function ac2_proj(row::Int, col::Int, ψ::MultilineMPS, (O, ϕ)::Tuple, envs)
+    return ac2_proj(col, ψ[row], (O[row], ϕ[row]), envs[row])
 end
 
 function ∂∂C(pos::Int, mps, operator::LinearCombination, cache)
diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl
index e93ddf961..2a5c26925 100644
--- a/src/algorithms/excitation/quasiparticleexcitation.jl
+++ b/src/algorithms/excitation/quasiparticleexcitation.jl
@@ -33,8 +33,8 @@ end
 ################################################################################
 
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renvs;
-                     num=1, solver=Defaults.linearsolver)
-    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; solver)
+                     num=1, kwargs...)
+    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; kwargs...)
     E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs)
     H_eff = @closure(ϕ -> effective_excitation_hamiltonian(H, ϕ, qp_envs(ϕ), E))
 
@@ -45,16 +45,16 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renv
     return Es, ϕs
 end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs;
-                     num=1, solver=Defaults.linearsolver)
+                     num=1, kwargs...)
     # Infer `renvs` in function body as it depends on `solver`.
-    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; solver=solver)
-    return excitations(H, alg, ϕ₀, lenvs, renvs; num, solver)
+    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
+    return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...)
 end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP;
-                     num=1, solver=Defaults.linearsolver)
+                     num=1, kwargs...)
     # Infer `lenvs` in function body as it depends on `solver`.
-    lenvs = environments(ϕ₀.left_gs, H; solver=solver)
-    return excitations(H, alg, ϕ₀, lenvs; num, solver)
+    lenvs = environments(ϕ₀.left_gs, H; kwargs...)
+    return excitations(H, alg, ϕ₀, lenvs; num, kwargs...)
 end
 
 """
@@ -83,19 +83,20 @@ Create and optimise 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);
-                     sector=one(sectortype(lmps)), num=1, solver=Defaults.linearsolver)
+                     sector=one(sectortype(lmps)), num=1, kwargs...)
     ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector, momentum)
-    return excitations(H, alg, ϕ₀, lenvs, renvs; num, solver)
+    return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...)
 end
 function excitations(H, alg::QuasiparticleAnsatz, momenta, lmps,
                      lenvs=environments(lmps, H), rmps=lmps,
                      renvs=lmps === rmps ? lenvs : environments(rmps, H);
-                     verbosity=Defaults.verbosity, num=1, solver=Defaults.linearsolver,
-                     sector=one(sectortype(lmps)), parallel=true)
+                     verbosity=Defaults.verbosity, num=1,
+                     sector=one(sectortype(lmps)), parallel=true, kwargs...)
     if parallel
         tasks = map(momenta) do momentum
             Threads.@spawn begin
-                E, ϕ = excitations(H, alg, momentum, lmps, lenvs, rmps, renvs; num, solver,
+                E, ϕ = excitations(H, alg, momentum, lmps, lenvs, rmps, renvs; num,
+                                   kwargs...,
                                    sector)
                 verbosity ≥ VERBOSE_CONV &&
                     @info "Found excitations for momentum = $(momentum)"
@@ -106,7 +107,7 @@ function excitations(H, alg::QuasiparticleAnsatz, momenta, lmps,
         fetched = fetch.(tasks)
     else
         fetched = map(momenta) do momentum
-            E, ϕ = excitations(H, alg, momentum, lmps, lenvs, rmps, renvs; num, solver,
+            E, ϕ = excitations(H, alg, momentum, lmps, lenvs, rmps, renvs; num, kwargs...,
                                sector)
             verbosity ≥ VERBOSE_CONV && @info "Found excitations for momentum = $(momentum)"
             return E, ϕ
@@ -167,9 +168,9 @@ end
 #                           Statmech Excitations                               #
 ################################################################################
 
-function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::Multiline{<:InfiniteQP},
-                     lenvs, renvs; num=1, solver=Defaults.linearsolver)
-    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; solver)
+function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::MultilineQP,
+                     lenvs, renvs; num=1, kwargs...)
+    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; kwargs...)
     function H_eff(ϕ′)
         ϕ = Multiline(ϕ′)
         return effective_excitation_hamiltonian(H, ϕ, qp_envs(ϕ)).data.data
@@ -183,8 +184,8 @@ function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::Multiline
 end
 
 function excitations(H::InfiniteMPO, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renvs;
-                     num=1, solver=Defaults.linearsolver)
-    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; solver)
+                     num=1, kwargs...)
+    qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; kwargs...)
     H_eff(ϕ) = effective_excitation_hamiltonian(H, ϕ, qp_envs(ϕ))
 
     Es, ϕs, convhist = eigsolve(H_eff, ϕ₀, num, :LM, alg.alg)
@@ -194,14 +195,14 @@ function excitations(H::InfiniteMPO, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP
     return Es, ϕs
 end
 
-function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::Multiline{<:InfiniteQP},
-                     lenvs; num=1, solver=Defaults.linearsolver)
+function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::MultilineQP,
+                     lenvs; num=1, kwargs...)
     # Infer `renvs` in function body as it depends on `solver`.
-    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; solver)
-    return excitations(H, alg, ϕ₀, lenvs, renvs; num, solver)
+    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
+    return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...)
 end
-function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::Multiline{<:InfiniteQP};
-                     num=1, solver=Defaults.linearsolver)
+function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::MultilineQP;
+                     num=1, kwargs...)
     # Infer `lenvs` in function body as it depends on `solver`.
     lenvs = environments(ϕ₀.left_gs, H; solver)
     return excitations(H, alg, ϕ₀, lenvs; num, solver)
@@ -211,15 +212,17 @@ function excitations(H::DenseMPO, alg::QuasiparticleAnsatz, momentum::Real,
                      lmps::InfiniteMPS,
                      lenvs=environments(lmps, H), rmps::InfiniteMPS=lmps,
                      renvs=lmps === rmps ? lenvs : environments(rmps, H);
-                     sector=one(sectortype(lmps)), num=1, solver=Defaults.linearsolver)
+                     sector=one(sectortype(lmps)), num=1, kwargs...)
     multiline_lmps = convert(MultilineMPS, lmps)
+    lenvs′ = Multiline([lenvs])
     if lmps === rmps
-        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs,
+        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs′,
                     multiline_lmps,
-                    lenvs; sector, num, solver)
+                    lenvs′; sector, num, kwargs...)
     else
-        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs,
-                    convert(MultilineMPS, rmps), renvs; sector, num, solver)
+        renvs′ = Multiline([renvs])
+        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs′,
+                    convert(MultilineMPS, rmps), renvs′; sector, num, kwargs...)
     end
 end
 
@@ -227,12 +230,12 @@ function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, momentum::Real,
                      lmps::MultilineMPS,
                      lenvs=environments(lmps, H), rmps=lmps,
                      renvs=lmps === rmps ? lenvs : environments(rmps, H);
-                     sector=one(sectortype(lmps)), num=1, solver=Defaults.linearsolver)
+                     sector=one(sectortype(lmps)), num=1, kwargs...)
     ϕ₀ = Multiline(map(1:size(lmps, 1)) do row
                        return LeftGaugedQP(rand, lmps[row], rmps[row]; sector, momentum)
                    end)
 
-    return excitations(H, alg, ϕ₀, lenvs, renvs; num, solver)
+    return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...)
 end
 
 ################################################################################
@@ -253,41 +256,10 @@ function effective_excitation_hamiltonian(H::Union{InfiniteMPOHamiltonian,
     return ϕ′
 end
 
-function effective_excitation_hamiltonian(H::MultilineMPO, ϕ::Multiline{<:InfiniteQP},
+function effective_excitation_hamiltonian(H::MultilineMPO, ϕ::MultilineQP,
                                           envs=environments(ϕ, H))
-    ϕ′ = Multiline(similar.(ϕ.data))
-    left_gs = ϕ.left_gs
-    right_gs = ϕ.right_gs
-
-    for row in 1:size(H, 1)
-        Bs = [ϕ[row][i] for i in 1:size(H, 2)]
-        for col in 1:size(H, 2)
-            en = @plansor conj(left_gs.AC[row, col][2 6; 4]) *
-                          leftenv(envs.lenvs, row, col, left_gs)[2 5; 3] *
-                          left_gs.AC[row + 1, col][3 7; 1] *
-                          H[row, col][5 6; 7 8] *
-                          rightenv(envs.lenvs, row, col, left_gs)[1 8; 4]
-
-            @plansor T[-1 -2; -3 -4] := leftenv(envs.lenvs, row, col, left_gs)[-1 5; 4] *
-                                        Bs[col][4 2; -3 1] *
-                                        H[row, col][5 -2; 2 3] *
-                                        rightenv(envs.renvs, row, col, right_gs)[1 3; -4]
-
-            @plansor T[-1 -2; -3 -4] += envs.lBs[row, col - 1][-1 4; -3 5] *
-                                        right_gs.AR[row, col][5 2; 1] *
-                                        H[row, col][4 -2; 2 3] *
-                                        rightenv(envs.renvs, row, col, right_gs)[1 3; -4]
-
-            @plansor T[-1 -2; -3 -4] += leftenv(envs.lenvs, row, col, left_gs)[-1 2; 1] *
-                                        left_gs.AL[row, col][1 3; 4] *
-                                        H[row, col][2 -2; 3 5] *
-                                        envs.rBs[row, col + 1][4 5; -3 -4]
-
-            ϕ′[row + 1][col] = T / en
-        end
-    end
-
-    return ϕ′
+    return Multiline(map(effective_excitation_hamiltonian,
+                         parent(H), parent(ϕ), parent(envs)))
 end
 
 function effective_excitation_hamiltonian(H::InfiniteMPO, ϕ::InfiniteQP,
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index 1a711db21..60ce1044a 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -103,9 +103,8 @@ end
 function expectation_value(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
                            envs::AbstractMPSEnvironments=environments(ψ, H))
     return sum(1:length(ψ)) do i
-        util = fill_data!(similar(ψ.AL[1], right_virtualspace(H, i)[end]), one)
-        @plansor GR[-1 -2; -3] := r_LL(ψ, i)[-1; -3] * util[-2]
-        return contract_mpo_expval(ψ.AL[i], leftenv(envs, i, ψ), H[i][:, 1, 1, end], GR)
+        return contract_mpo_expval(ψ.AC[i], envs.GLs[i], H[i][:, 1, 1, end],
+                                   envs.GRs[i][end])
     end
 end
 
@@ -120,14 +119,12 @@ end
 function expectation_value(ψ::InfiniteMPS, mpo::InfiniteMPO, envs...)
     return expectation_value(convert(MultilineMPS, ψ), convert(MultilineMPO, mpo), envs...)
 end
-function expectation_value(ψ::MultilineMPS, O::MultilineMPO{<:Union{DenseMPO,SparseMPO}},
-                           envs::InfiniteMPOEnvironments=environments(ψ, O))
+function expectation_value(ψ::MultilineMPS, O::MultilineMPO{<:InfiniteMPO},
+                           envs::MultilineEnvironments=environments(ψ, O))
     return prod(product(1:size(ψ, 1), 1:size(ψ, 2))) do (i, j)
-        GL = leftenv(envs, i, j, ψ)
-        GR = rightenv(envs, i, j, ψ)
-        @plansor GL[1 2; 3] * O[i, j][2 4; 6 5] *
-                 ψ.AC[i, j][3 6; 7] * GR[7 5; 8] *
-                 conj(ψ.AC[i + 1, j][1 4; 8])
+        GL = envs[i].GLs[j]
+        GR = envs[i].GRs[j]
+        return contract_mpo_expval(ψ.AC[i, j], GL, O[i, j], GR, ψ.AC[i + 1, j])
     end
 end
 function expectation_value(ψ::MultilineMPS, mpo::MultilineMPO, envs...)
diff --git a/src/algorithms/grassmann.jl b/src/algorithms/grassmann.jl
index b1efe5a42..8ff4a3ea6 100644
--- a/src/algorithms/grassmann.jl
+++ b/src/algorithms/grassmann.jl
@@ -11,260 +11,224 @@ The module exports nothing, and all references to it should be qualified, e.g.
 module GrassmannMPS
 
 using ..MPSKit
+using ..MPSKit: AbstractMPSEnvironments, InfiniteEnvironments, MultilineEnvironments, ∂∂AC
 using TensorKit
 using OhMyThreads
 import TensorKitManifolds.Grassmann
-
-function TensorKit.rmul!(a::Grassmann.GrassmannTangent, b::AbstractTensorMap)
-    #rmul!(a.Z,b);
-    Base.setfield!(a, :Z, a.Z * b)
-    Base.setfield!(a, :U, nothing)
-    Base.setfield!(a, :S, nothing)
-    Base.setfield!(a, :V, nothing)
-    return a
-end
-function Base.:/(a::Grassmann.GrassmannTangent, b::AbstractTensorMap)
-    return Grassmann.GrassmannTangent(a.W, a.Z / b)
+using TensorKitManifolds.Grassmann: GrassmannTangent, checkbase
+using VectorInterface: VectorInterface, One
+using OptimKit: OptimKit
+
+# utility
+# -------
+function rmul(Δ::GrassmannTangent, C::AbstractTensorMap)
+    return GrassmannTangent(Δ.W, Δ.Z * C)
+end
+
+# TODO: implement VectorInterface support for TensorKitManifolds
+function add!(x::GrassmannTangent, y::GrassmannTangent, α::Number=One(), β::Number=One())
+    checkbase(x, y)
+    VectorInterface.add!(x.Z, y.Z, α, β)
+    Base.setfield!(x, :U, nothing)
+    Base.setfield!(x, :S, nothing)
+    Base.setfield!(x, :V, nothing)
+    return x
+end
+function add!(x::AbstractArray, y::AbstractArray, α::Number=One(), β::Number=One())
+    return (add!.(x, y, α, β); x)
+end
+
+function scale!(x::GrassmannTangent, α::Number)
+    VectorInterface.scale!(x.Z, α)
+    if !isnothing(Base.getfield(x, :S))
+        sα = sign(α)
+        if sα != 1
+            VectorInterface.scale!(x.U, sα)
+        end
+        VectorInterface.scale!(x.S, abs(α))
+    end
+    return x
 end
+scale!(x::AbstractArray, α::Number) = (scale!.(x, α); x)
 
-# preconditioned gradient
-struct PrecGrad{A,B}
-    Pg::A
-    g::A
-    rho::B
-end
+# manifold methods
+# ----------------
+"""
+    inner(state, g1, g2)
 
-function PrecGrad(v::Grassmann.GrassmannTangent)
-    return PrecGrad(v, v, isometry(storagetype(v.Z), domain(v.Z), domain(v.Z)))
-end
-PrecGrad(v::Grassmann.GrassmannTangent, rho) = PrecGrad(v / rho, v, rho)
+The inner product between tangent vectors on the left-canonical MPS manifold.
+This is given by the (Euclidean) inner product between the tangent vectors at each site.
+Because the cost function is assumed holomorphic, we take the twice the real part of the result.
 
-Grassmann.base(g::PrecGrad) = Grassmann.base(g.Pg)
+For the effective inner product lifted from the Hilbert space metric,
+see also [`precondition`](@ref).
+"""
+inner(state, g1, g2) = 2 * real(sum(x -> Grassmann.inner(x...), zip(state.AL, g1, g2)))
 
-function inner(g1::PrecGrad, g2::PrecGrad, rho=one(g1.rho))
-    Grassmann.base(g1) == Grassmann.base(g2) || throw(ArgumentError("incompatible base"))
-    if g1.rho == rho
-        dot(g1.g.Z, g2.Pg.Z)
-    elseif g2.rho == rho
-        dot(g1.Pg.Z, g2.g.Z)
-    else
-        dot(g1.Pg.Z, g2.Pg.Z * rho)
-    end
-end
+"""
+    precondition(state, g)
 
-Base.:*(g::PrecGrad, alpha::Number) = PrecGrad(g.Pg * alpha, g.g * alpha, g.rho)
-function Base.:+(a::PrecGrad, b::PrecGrad)
-    if a.rho == b.rho
-        PrecGrad(a.Pg + b.Pg, a.g + b.g, a.rho)
-    else
-        PrecGrad(a.Pg + b.Pg)
+In order to obtain an effective inner product between tangent vectors that is lifted from
+the inner product of the Hilbert space of the state, we additionally include the inverse
+metric: ``g ↦ g / ρ`` where `ρ` is the (regularised) right fixed point of the MPS transfer
+matrix.
+"""
+function precondition(state, g)
+    g′ = similar(g)
+    rtolmin = eps(real(scalartype(state)))^(3 / 4)
+    tforeach(eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        rtol = max(rtolmin, norm(g[i]))
+        ρ = rho_inv_regularized(state.C[i]; rtol)
+        g′[i] = rmul(g[i], ρ)
+        return nothing
     end
+    return g′
 end
 
-struct ManifoldPoint{T,E,G,C}
-    state::T # the state at that point
-    envs::E # the environments
-    g::G # the MPS gradient, which is not equivalent with the grassmann gradient!
-    Rhoreg::C # the regularized density matrices
-end
+"""
+    retract(state, g, α) -> state′, ξ
 
-function ManifoldPoint(state::FiniteMPS, envs)
-    al_d = similar(state.AL)
-    O = envs.operator
-    for i in 1:length(state)
-        al_d[i] = MPSKit.∂∂AC(i, state, O, envs) * state.AC[i]
+Retract a state a distance `α` along a direction `g`, obtaining a new state and the local tangent vector. 
+"""
+function retract(state::FiniteMPS, g, α::Real)
+    state′ = copy(state)
+    h = map(eachindex(state)) do i
+        AL′, ξ = Grassmann.retract(state.AL[i], g[i], α)
+        state′.AC[i] = (AL′, state.C[i])
+        return ξ
     end
-    g = Grassmann.project.(al_d, state.AL)
-
-    Rhoreg = Vector{eltype(state.C)}(undef, length(state))
-    δmin = sqrt(eps(real(scalartype(state))))
-    tmap!(Rhoreg, 1:length(state); scheduler=MPSKit.Defaults.scheduler[]) do i
-        return regularize(state.C[i], max(norm(g[i]) / 10, δmin))
+    normalize!(state′)
+    return state′, h
+end
+function retract(state::InfiniteMPS, g, α::Real)
+    AL′ = similar(state.AL)
+    g′ = similar(g)
+    tforeach(eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        AL′[i], g′[i] = Grassmann.retract(state.AL[i], g[i], α)
+        return nothing
     end
-
-    return ManifoldPoint(state, envs, g, Rhoreg)
-end
-function ManifoldPoint(state::InfiniteMPS, envs)
-    δmin = sqrt(eps(real(scalartype(state))))
-    Tg = Core.Compiler.return_type(Grassmann.project,
-                                   Tuple{eltype(state.AL),eltype(state.AL)})
-    g = similar(state.AL, Tg)
-    ρ = similar(state.C)
-
-    MPSKit.check_recalculate!(envs, state)
-    tforeach(1:length(state); scheduler=MPSKit.Defaults.scheduler[]) do i
-        AC′ = MPSKit.∂∂AC(i, state, envs.operator, envs) * state.AC[i]
-        g[i] = Grassmann.project(AC′, state.AL[i])
-        ρ[i] = regularize(state.C[i], max(norm(g[i]) / 10, δmin))
+    state′ = InfiniteMPS(AL′, state.C[end])
+    return state′, g′
+end
+function retract(state::MultilineMPS, g, α::Real)
+    AL′ = similar(state.AL)
+    g′ = similar(g)
+    tforeach(eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        AL′[i], g′[i] = Grassmann.retract(state.AL[i], g[i], α)
         return nothing
     end
-    return ManifoldPoint(state, envs, g, ρ)
+    state′ = MultilineMPS(AL′, state.C[:, end])
+    return state′, g′
 end
 
-function ManifoldPoint(state::MultilineMPS, envs)
-    # FIXME: add support for unitcells
-    @assert length(state.AL) == 1 "GradientGrassmann only supports MultilineMPS without unitcells for now"
-
-    # TODO: this really should not use the operator from the environment
-    f = expectation_value(state, envs.operator, envs)
-    imag(f) > MPSKit.Defaults.tol && @warn "MPO might not be Hermitian $f"
-    real(f) > 0 || @warn "MPO might not be positive definite $f"
-
-    grad = map(CartesianIndices(state.AC)) do I
-        AC′ = MPSKit.∂∂AC(I, state, envs.operator, envs) * state.AC[I]
-        # the following formula is wrong when unitcells are involved
-        # actual costfunction should be F = -log(prod(f)) => ∂F = -2 * g / |f|
-        return rmul!(Grassmann.project(AC′, state.AL[I]), -2 / f)
-    end
+"""
+    transport!(h, state, g, α, state′) -> h
 
-    δmin = sqrt(eps(real(scalartype(state))))
-    ρ_regularized = map(state.C, grad) do ρ, g
-        return regularize(ρ, max(norm(g) / 10, δmin))
+In-place transport of a tangent vector `g` at a point `state`, to a new point `state′`.
+"""
+function transport!(h, state, g, α::Real, state′)
+    tforeach(eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        h[i] = Grassmann.transport!(h[i], state.AL[i], g[i], α, state′.AL[i])
+        return nothing
     end
-
-    return ManifoldPoint(state, envs, grad, ρ_regularized)
+    return h
 end
 
 """
-Compute the expectation value, and its gradient with respect to the tensors in the unit
-cell as tangent vectors on Grassmann manifolds.
-"""
-function fg(x::ManifoldPoint{T}) where {T<:Union{InfiniteMPS,FiniteMPS}}
-    # the gradient I want to return is the preconditioned gradient!
-    Tg = Core.Compiler.return_type(PrecGrad, Tuple{eltype(x.g),eltype(x.Rhoreg)})
-    g_prec = similar(x.g, Tg)
-    tmap!(g_prec, eachindex(x.g); scheduler=MPSKit.Defaults.scheduler[]) do i
-        return PrecGrad(rmul!(copy(x.g[i]), x.state.C[i]'), x.Rhoreg[i])
-    end
+    fg(state, operator, envs=environments(state, operator))
 
-    # TODO: the operator really should not be part of the environments, and this should
-    # be passed as an explicit argument
-    f = expectation_value(x.state, x.envs.operator, x.envs)
+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)) where {O<:FiniteMPOHamiltonian}
+    f = expectation_value(state, operator, envs)
     isapprox(imag(f), 0; atol=eps(abs(f))^(3 / 4)) || @warn "MPO might not be Hermitian: $f"
-
-    return real(f), g_prec
-end
-function fg(x::ManifoldPoint{<:MultilineMPS})
-    @assert length(x.state) == 1 "GradientGrassmann only supports MultilineMPS without unitcells for now"
-    # the gradient I want to return is the preconditioned gradient!
-    g_prec = map(enumerate(x.g)) do (i, cg)
-        return PrecGrad(rmul!(copy(cg), x.state.C[i]'), x.Rhoreg[i])
+    gs = map(1:length(state)) do i
+        AC′ = ∂∂AC(i, state, operator, envs) * state.AC[i]
+        g = Grassmann.project(AC′, state.AL[i])
+        return rmul(g, state.C[i]')
     end
-
-    # TODO: the operator really should not be part of the environments, and this should
-    # be passed as an explicit argument
-    f = expectation_value(x.state, x.envs.operator, x.envs)
+    return real(f), gs
+end
+function fg(state::InfiniteMPS, operator::Union{O,LazySum{O}},
+            envs::AbstractMPSEnvironments=environments(state, operator)) where {O<:InfiniteMPOHamiltonian}
+    recalculate!(envs, state, operator, state)
+    f = expectation_value(state, operator, envs)
     isapprox(imag(f), 0; atol=eps(abs(f))^(3 / 4)) || @warn "MPO might not be Hermitian: $f"
-    real(f) > 0 || @warn "MPO might not be positive definite: $f"
 
-    return -log(real(f)), g_prec[:]
+    A = Core.Compiler.return_type(Grassmann.project, Tuple{eltype(state),eltype(state)})
+    gs = Vector{A}(undef, length(state))
+    tmap!(gs, 1:length(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        AC′ = ∂∂AC(i, state, operator, envs) * state.AC[i]
+        g = Grassmann.project(AC′, state.AL[i])
+        return rmul(g, state.C[i]')
+    end
+    return real(f), gs
 end
+function fg(state::InfiniteMPS, operator::Union{O,LazySum{O}},
+            envs::AbstractMPSEnvironments=environments(state, operator)) where {O<:InfiniteMPO}
+    recalculate!(envs, state, operator, state)
+    f = expectation_value(state, operator, envs)
+    isapprox(imag(f), 0; atol=eps(abs(f))^(3 / 4)) || @warn "MPO might not be Hermitian: $f"
 
-"""
-Retract a left-canonical MultilineMPS along Grassmann tangent `g` by distance `alpha`.
-"""
-function retract(x::ManifoldPoint{<:MultilineMPS}, tg, alpha)
-    g = reshape(tg, size(x.state))
-
-    nal = similar(x.state.AL)
-    h = similar(tg)
-    tmap!(h, eachindex(g); scheduler=MPSKit.Defaults.scheduler[]) do i
-        nal[i], th = Grassmann.retract(x.state.AL[i], g[i].Pg, alpha)
-        return PrecGrad(th)
+    A = Core.Compiler.return_type(Grassmann.project, Tuple{eltype(state),eltype(state)})
+    gs = Vector{A}(undef, length(state))
+    tmap!(gs, eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        AC′ = ∂∂AC(i, state, operator, envs) * state.AC[i]
+        g = rmul!(Grassmann.project(AC′, state.AL[i]), -inv(f))
+        return rmul(g, state.C[i]')
     end
-
-    nstate = MPSKit.MultilineMPS(nal, x.state.C[:, end])
-    newpoint = ManifoldPoint(nstate, x.envs)
-
-    return newpoint, h[:]
+    return -log(real(f)), gs
 end
+function fg(state::MultilineMPS, operator::MultilineMPO,
+            envs::MultilineEnvironments=environments(state, operator))
+    @assert length(state) == 1 "not implemented"
+    recalculate!(envs, state, operator, state)
+    f = expectation_value(state, operator, envs)
+    isapprox(imag(f), 0; atol=eps(abs(f))^(3 / 4)) || @warn "MPO might not be Hermitian: $f"
 
-"""
-Retract a left-canonical infinite MPS along Grassmann tangent `g` by distance `alpha`.
-"""
-function retract(x::ManifoldPoint{<:InfiniteMPS}, g, alpha)
-    state = x.state
-    envs = x.envs
-    nal = similar(state.AL)
-    h = similar(g)  # The tangent at the end-point
-
-    tmap!(h, eachindex(g); scheduler=MPSKit.Defaults.scheduler[]) do i
-        nal[i], th = Grassmann.retract(state.AL[i], g[i].Pg, alpha)
-        return PrecGrad(th)
+    A = Core.Compiler.return_type(Grassmann.project, Tuple{eltype(state),eltype(state)})
+    gs = Matrix{A}(undef, size(state))
+    tforeach(eachindex(state); scheduler=MPSKit.Defaults.scheduler[]) do i
+        AC′ = ∂∂AC(i, state, operator, envs) * state.AC[i]
+        g = rmul!(Grassmann.project(AC′, state.AL[i]), -inv(f))
+        gs[i] = rmul(g, state.C[i]')
+        return nothing
     end
-
-    nstate = InfiniteMPS(nal, state.C[end])
-
-    newpoint = ManifoldPoint(nstate, envs)
-
-    return newpoint, h
+    return -log(real(f)), gs
 end
 
 """
-Retract a left-canonical finite MPS along Grassmann tangent `g` by distance `alpha`.
-"""
-function retract(x::ManifoldPoint{<:FiniteMPS}, g, alpha)
-    state = x.state
-    envs = x.envs
+    rho_inv_regularized(C; rtol=eps(real(scalartype(C)))^(3/4))
 
-    y = copy(state)  # The end-point
-    h = similar(g)  # The tangent at the end-point
-    for i in 1:length(g)
-        yal, th = Grassmann.retract(state.AL[i], g[i].Pg, alpha)
-        h[i] = PrecGrad(th)
-        y.AC[i] = (yal, state.C[i])
-    end
-    normalize!(y)
-
-    n_y = ManifoldPoint(y, envs)
-
-    return n_y, h
+Compute the (regularized) inverse of the MPS fixed point `ρ = C * C'`.
+Here we use the Tikhonov regularization, i.e. `inv(ρ) = inv(C * C' + δ²1)`,
+where the regularization parameter is `δ = rtol * norm(C)`.
+"""
+function rho_inv_regularized(C; rtol=eps(real(scalartype(C)))^(3 / 4))
+    U, S, _ = tsvd(C)
+    return U * pinv_tikhonov!!(S; rtol) * U'
 end
 
-"""
-Transport a tangent vector `h` along the retraction from `x` in direction `g` by distance
-`alpha`. `xp` is the end-point of the retraction.
-"""
-function transport!(h, x, g, alpha, xp)
-    tforeach(1:length(h); scheduler=MPSKit.Defaults.scheduler[]) do i
-        return h[i] = PrecGrad(Grassmann.transport!(h[i].Pg, x.state.AL[i], g[i].Pg, alpha,
-                                                    xp.state.AL[i]))
+function pinv_tikhonov!!(S::DiagonalTensorMap{<:Real}; rtol=zero(scalartype(S)))
+    δ² = (rtol * maximum(maximum ∘ last, blocks(S); init=zero(scalartype(S))))^2
+    for (_, b) in blocks(S)
+        b.diag .= inv.(b.diag .^ 2 .+ δ²)
     end
-    return h
+    return S
 end
-
-"""
-Euclidean inner product between two Grassmann tangents of an infinite MPS.
-"""
-function inner(x, g1, g2)
-    return 2 * real(sum(((a, b, c),) -> inner(b, c, a), zip(x.Rhoreg, g1, g2)))
+# TensorKit v0.13 still outputs AbstractTensorMap so define fallback
+function pinv_tikhonov!!(S::AbstractTensorMap{<:Real}; rtol=zero(scalartype(S)))
+    δ² = (rtol * norm(S, Inf))^2
+    return inv(S^2 + δ² * one(S))
 end
 
-"""
-Scale a tangent vector by scalar `alpha`.
-"""
-scale!(g, alpha) = g .* alpha
-
-"""
-Add two tangents vectors, scaling the latter by `alpha`.
-"""
-add!(g1, g2, alpha) = g1 + g2 .* alpha
-
-"""
-Take the L2 Tikhonov regularised of a matrix `m`.
-
-The regularisation parameter is the larger of `delta` (the optional argument that defaults
-to zero) and square root of machine epsilon.
-"""
-function regularize(m, delta=zero(scalartype(m)))
-    U, S, V = tsvd(m)
-
-    #Sreg = real(S*sqrt(one(S) + delta^2*one(S)*norm(S,Inf)^2/S^2));#
-    Sreg = S^2 + (norm(S, Inf) * delta)^2 * one(S)
-
-    Mreg = U * Sreg * U'
-
-    return Mreg
+# utility test function
+function optimtest(ψ, O, envs=environments(ψ, O); alpha=-0.1:0.001:0.1,
+                   retract=retract,
+                   inner=inner)
+    _fg(x) = fg(x, O, envs)
+    return OptimKit.optimtest(_fg, ψ; alpha, retract, inner)
 end
 
 end  # module GrassmannMPS
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 74be9efa5..0d1c922ce 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -26,7 +26,7 @@ function DMRG(; tol=Defaults.tol, maxiter=Defaults.maxiter, eigalg=(;),
 end
 
 function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs=environments(ψ, H))
-    ϵs = map(pos -> calc_galerkin(ψ, pos, envs), 1:length(ψ))
+    ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
     ϵ = maximum(ϵs)
     log = IterLog("DMRG")
 
@@ -39,7 +39,7 @@ function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs=environment
             for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
                 h = ∂∂AC(pos, ψ, H, envs)
                 _, vec = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
-                ϵs[pos] = max(ϵs[pos], calc_galerkin(ψ, pos, envs))
+                ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
                 ψ.AC[pos] = vec
             end
             ϵ = maximum(ϵs)
@@ -93,7 +93,7 @@ function DMRG2(; tol=Defaults.tol, maxiter=Defaults.maxiter, eigalg=(;),
 end
 
 function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environments(ψ, H))
-    ϵs = map(pos -> calc_galerkin(ψ, pos, envs), 1:length(ψ))
+    ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
     ϵ = maximum(ϵs)
     log = IterLog("DMRG2")
 
diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl
index df0c3d15d..56117a86a 100644
--- a/src/algorithms/groundstate/gradient_grassmann.jl
+++ b/src/algorithms/groundstate/gradient_grassmann.jl
@@ -54,16 +54,16 @@ function find_groundstate(ψ::S, H, alg::GradientGrassmann,
         @warn "This is not fully supported - split the mps up in a sum of mps's and optimize seperately"
     normalize!(ψ)
 
-    #optimtest(GrassmannMPS.fg,(ψ,envs);alpha=-0.01:0.001:0.01,retract=GrassmannMPS.retract,inner=GrassmannMPS.inner)
-    x, _, _, _, normgradhistory = optimize(GrassmannMPS.fg,
-                                           GrassmannMPS.ManifoldPoint(ψ, envs),
+    fg(x) = GrassmannMPS.fg(x, H, envs)
+    x, _, _, _, normgradhistory = optimize(fg, ψ,
                                            alg.method;
-                                           (transport!)=GrassmannMPS.transport!,
-                                           retract=GrassmannMPS.retract,
-                                           inner=GrassmannMPS.inner,
-                                           (scale!)=GrassmannMPS.scale!,
-                                           (add!)=GrassmannMPS.add!,
-                                           (finalize!)=alg.finalize!,
+                                           GrassmannMPS.transport!,
+                                           GrassmannMPS.retract,
+                                           GrassmannMPS.inner,
+                                           GrassmannMPS.scale!,
+                                           GrassmannMPS.add!,
+                                           GrassmannMPS.precondition,
+                                           alg.finalize!,
                                            isometrictransport=true)
-    return x.state, x.envs, normgradhistory[end]
+    return x, envs, normgradhistory[end]
 end
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index 207cefd21..7c6575d26 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -23,10 +23,9 @@ function IDMRG1(; tol=Defaults.tol, tol_gauge=Defaults.tolgauge, eigalg=(;),
     return IDMRG1{typeof(eigalg′)}(tol, tol_gauge, eigalg′, maxiter, verbosity)
 end
 
-function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG1, oenvs=environments(ost, H))
-    ϵ::Float64 = calc_galerkin(ost, oenvs)
+function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG1, envs=environments(ost, H))
+    ϵ::Float64 = calc_galerkin(ost, H, ost, envs)
     ψ = copy(ost)
-    envs = IDMRGEnvironments(ost, oenvs)
     log = IterLog("IDMRG")
 
     LoggingExtras.withlevel(; alg.verbosity) do
@@ -40,7 +39,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG1, oenvs=environments(o
                 h = ∂∂AC(pos, ψ, H, envs)
                 _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
                 ψ.AL[pos], ψ.C[pos] = leftorth!(ψ.AC[pos])
-                update_leftenv!(envs, ψ, H, pos + 1)
+                transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
             end
 
             # right to left sweep
@@ -51,7 +50,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG1, oenvs=environments(o
                 ψ.C[pos - 1], temp = rightorth!(_transpose_tail(ψ.AC[pos]))
                 ψ.AR[pos] = _transpose_front(temp)
 
-                update_rightenv!(envs, ψ, H, pos - 1)
+                transfer_rightenv!(envs, ψ, H, ψ, pos - 1)
             end
 
             ϵ = norm(C_current - ψ.C[0])
@@ -69,8 +68,9 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG1, oenvs=environments(o
     end
 
     nst = InfiniteMPS(ψ.AR[1:end]; tol=alg.tol_gauge)
-    nenvs = environments(nst, H; solver=oenvs.solver)
-    return nst, nenvs, ϵ
+    recalculate!(envs, nst, H, nst)
+
+    return nst, envs, ϵ
 end
 
 """
@@ -101,12 +101,11 @@ function IDMRG2(; tol=Defaults.tol, tol_gauge=Defaults.tolgauge, eigalg=(;),
     return IDMRG2{typeof(eigalg′)}(tol, tol_gauge, eigalg′, maxiter, verbosity, trscheme)
 end
 
-function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(ost, H))
+function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(ost, H))
     length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    ϵ::Float64 = calc_galerkin(ost, oenvs)
+    ϵ::Float64 = calc_galerkin(ost, H, ost, envs)
 
     ψ = copy(ost)
-    envs = IDMRGEnvironments(ost, oenvs)
     log = IterLog("IDMRG2")
 
     LoggingExtras.withlevel(; alg.verbosity) do
@@ -129,8 +128,8 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(o
                 ψ.AR[pos + 1] = _transpose_front(ar)
                 ψ.AC[pos + 1] = _transpose_front(c * ar)
 
-                update_leftenv!(envs, ψ, H, pos + 1)
-                update_rightenv!(envs, ψ, H, pos)
+                transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
+                transfer_rightenv!(envs, ψ, H, ψ, pos)
             end
 
             # update the edge
@@ -152,8 +151,8 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(o
             C_current = complex(c)
 
             # update environments
-            update_leftenv!(envs, ψ, H, 1)
-            update_rightenv!(envs, ψ, H, 0)
+            transfer_leftenv!(envs, ψ, H, ψ, 1)
+            transfer_rightenv!(envs, ψ, H, ψ, 0)
 
             # sweep from right to left
             for pos in (length(ψ) - 1):-1:1
@@ -170,8 +169,8 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(o
                 ψ.AR[pos + 1] = _transpose_front(ar)
                 ψ.AC[pos + 1] = _transpose_front(c * ar)
 
-                update_leftenv!(envs, ψ, H, pos + 1)
-                update_rightenv!(envs, ψ, H, pos)
+                transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
+                transfer_rightenv!(envs, ψ, H, ψ, pos)
             end
 
             # update the edge
@@ -188,8 +187,8 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(o
             ψ.AR[1] = _transpose_front(ar)
             ψ.AC[1] = _transpose_front(c * ar)
 
-            update_leftenv!(envs, ψ, H, 1)
-            update_rightenv!(envs, ψ, H, 0)
+            transfer_leftenv!(envs, ψ, H, ψ, 1)
+            transfer_rightenv!(envs, ψ, H, ψ, 0)
 
             # update error
             smallest = infimum(_firstspace(C_current), _firstspace(c))
@@ -210,6 +209,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, oenvs=environments(o
     end
 
     ψ′ = InfiniteMPS(ψ.AR[1:end]; tol=alg.tol_gauge)
-    nenvs = environments(ψ′, H; solver=oenvs.solver)
-    return ψ′, nenvs, ϵ
+    recalculate!(envs, ψ′, H, ψ′)
+
+    return ψ′, envs, ϵ
 end
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 7f96bfae3..504311207 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -28,10 +28,12 @@ end
 
 function find_groundstate(ψ::InfiniteMPS, H, alg::VUMPS, envs=environments(ψ, H))
     # initialization
-    ϵ::Float64 = calc_galerkin(ψ, envs)
-    temp_ACs = similar.(ψ.AC)
     scheduler = Defaults.scheduler[]
     log = IterLog("VUMPS")
+    ϵ::Float64 = calc_galerkin(ψ, H, ψ, envs)
+    temp_ACs = similar.(ψ.AC)
+    alg_environments = updatetol(alg.alg_environments, 0, ϵ)
+    recalculate!(envs, ψ, H, ψ; alg_environments.tol)
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ, sum(expectation_value(ψ, H, envs)))
@@ -45,21 +47,21 @@ function find_groundstate(ψ::InfiniteMPS, H, alg::VUMPS, envs=environments(ψ,
             ψ = InfiniteMPS(temp_ACs, ψ.C[end]; alg_gauge.tol, alg_gauge.maxiter)
 
             alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-            recalculate!(envs, ψ; alg_environments.tol)
+            recalculate!(envs, ψ, H, ψ; alg_environments.tol)
 
             ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ),typeof(envs)}
 
-            ϵ = calc_galerkin(ψ, envs)
+            ϵ = calc_galerkin(ψ, H, ψ, envs)
 
             # breaking conditions
             if ϵ <= alg.tol
-                @infov 2 logfinish!(log, iter, ϵ, sum(expectation_value(ψ, H, envs)))
+                @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
                 break
             end
             if iter == alg.maxiter
-                @warnv 1 logcancel!(log, iter, ϵ, sum(expectation_value(ψ, H, envs)))
+                @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs))
             else
-                @infov 3 logiter!(log, iter, ϵ, sum(expectation_value(ψ, H, envs)))
+                @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs))
             end
         end
     end
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index 6e1f18f2e..1a4816616 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -61,7 +61,7 @@ function propagator(A::AbstractFiniteMPS, z::Number, H::FiniteMPOHamiltonian,
             ϵ = 0.0
 
             for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = ac_proj(i, init, mixedenvs)
+                tos = ac_proj(i, init, A, mixedenvs)
 
                 H_AC = ∂∂AC(i, init, H, h_envs)
                 AC = init.AC[i]
@@ -116,7 +116,7 @@ function propagator(A::AbstractFiniteMPS, z, H::FiniteMPOHamiltonian,
             ϵ = 0.0
 
             for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = ac_proj(i, init, mixedenvs)
+                tos = ac_proj(i, init, A, mixedenvs)
                 H1_AC = ∂∂AC(i, init, H, envs1)
                 H2_AC = ∂∂AC(i, init, H2, envs2)
                 H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1))
@@ -141,7 +141,7 @@ function propagator(A::AbstractFiniteMPS, z, H::FiniteMPOHamiltonian,
         end
     end
 
-    a = dot(ac_proj(1, init, mixedenvs), init.AC[1])
+    a = dot(ac_proj(1, init, A, mixedenvs), init.AC[1])
     cb = leftenv(envs1, 1, A) * TransferMatrix(init.AL, H[1:length(A.AL)], A.AL)
     b = zero(a)
     for i in 1:length(cb)
diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl
index 5712b7696..0540f0b72 100644
--- a/src/algorithms/statmech/gradient_grassmann.jl
+++ b/src/algorithms/statmech/gradient_grassmann.jl
@@ -1,24 +1,17 @@
-function leading_boundary(state::InfiniteMPS, H::DenseMPO, alg::GradientGrassmann,
-                          envs=environments(state, H))
-    (multi, envs, err) = leading_boundary(convert(MultilineMPS, state),
-                                          convert(MultilineMPO, H), alg, envs)
-    state = convert(InfiniteMPS, multi)
-    return (state, envs, err)
-end
-
-function leading_boundary(state::MultilineMPS, H, alg::GradientGrassmann,
-                          envs=environments(state, H))
-    res = optimize(GrassmannMPS.fg,
-                   GrassmannMPS.ManifoldPoint(state, envs),
-                   alg.method;
-                   (transport!)=GrassmannMPS.transport!,
-                   retract=GrassmannMPS.retract,
-                   inner=GrassmannMPS.inner,
-                   (scale!)=GrassmannMPS.scale!,
-                   (add!)=GrassmannMPS.add!,
-                   (finalize!)=alg.finalize!,
-                   #precondition = GrassmannMPS.precondition,
-                   isometrictransport=true)
-    (x, fx, gx, numfg, normgradhistory) = res
-    return x.state, x.envs, normgradhistory[end]
+function leading_boundary(state::MultilineMPS,
+                          operator::MultilineMPO,
+                          alg::GradientGrassmann,
+                          envs::MultilineEnvironments=environments(state, H))
+    fg(x) = GrassmannMPS.fg(x, operator, envs)
+    x, _, _, _, normgradhistory = optimize(fg, state,
+                                           alg.method;
+                                           GrassmannMPS.transport!,
+                                           GrassmannMPS.retract,
+                                           GrassmannMPS.inner,
+                                           GrassmannMPS.scale!,
+                                           GrassmannMPS.add!,
+                                           GrassmannMPS.precondition,
+                                           alg.finalize!,
+                                           isometrictransport=true)
+    return x, envs, normgradhistory[end]
 end
diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl
new file mode 100644
index 000000000..1fc759dd5
--- /dev/null
+++ b/src/algorithms/statmech/leading_boundary.jl
@@ -0,0 +1,33 @@
+
+@doc """
+    leading_boundary(ψ, O, [environments]; kwargs...)
+    leading_boundary(ψ, O, algorithm, environments)
+
+Compute the leading boundary MPS for operator `O` with initial guess `ψ`. If not specified, an
+optimization algorithm will be attempted based on the supplied keywords.
+
+## Arguments
+- `ψ::AbstractMPS`: initial guess
+- `O::AbstractMPO`: operator for which to find the leading_boundary
+- `[environments]`: MPS environment manager
+- `algorithm`: optimization algorithm
+
+## Keywords
+- `tol::Float64`: tolerance for convergence criterium
+- `maxiter::Int`: maximum amount of iterations
+- `verbosity::Int`: display progress information
+""" leading_boundary
+
+# TODO: alg selector
+
+# implementation always in terms of Multiline objects
+function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg,
+                          envs=environments(state, operator))
+    state_multi = convert(MultilineMPS, state)
+    operator_multi = convert(MultilineMPO, operator)
+    envs_multi = Multiline([envs])
+    state_multi′, envs_multi′, err = leading_boundary(state_multi, operator_multi, alg,
+                                                      envs_multi)
+    state′ = convert(InfiniteMPS, state_multi′)
+    return state′, envs, err
+end
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 8116e0869..310c4d7cf 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -26,7 +26,7 @@ end
 
 function leading_boundary(ψ::MultilineMPS, O::MultilineMPO, alg::VOMPS,
                           envs=environments(ψ, O))
-    ϵ::Float64 = calc_galerkin(ψ, envs)
+    ϵ::Float64 = calc_galerkin(ψ, O, ψ, envs)
     temp_ACs = similar.(ψ.AC)
     scheduler = Defaults.scheduler[]
     log = IterLog("VOMPS")
@@ -42,11 +42,11 @@ function leading_boundary(ψ::MultilineMPS, O::MultilineMPO, alg::VOMPS,
             ψ = MultilineMPS(temp_ACs, ψ.C[:, end]; alg_gauge.tol, alg_gauge.maxiter)
 
             alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-            recalculate!(envs, ψ; alg_environments.tol)
+            recalculate!(envs, ψ, O, ψ; alg_environments.tol)
 
             ψ, envs = alg.finalize(iter, ψ, O, envs)::Tuple{typeof(ψ),typeof(envs)}
 
-            ϵ = calc_galerkin(ψ, envs)
+            ϵ = calc_galerkin(ψ, O, ψ, envs)
 
             if ϵ <= alg.tol
                 @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, O, envs))
diff --git a/src/algorithms/statmech/vumps.jl b/src/algorithms/statmech/vumps.jl
index e9baf7d6b..aa0987b35 100644
--- a/src/algorithms/statmech/vumps.jl
+++ b/src/algorithms/statmech/vumps.jl
@@ -18,10 +18,13 @@ function leading_boundary(ψ::InfiniteMPS, H, alg, envs=environments(ψ, H))
 end
 
 function leading_boundary(ψ::MultilineMPS, H, alg::VUMPS, envs=environments(ψ, H))
-    ϵ::Float64 = calc_galerkin(ψ, envs)
-    temp_ACs = similar.(ψ.AC)
-    scheduler = Defaults.scheduler[]
+    # initialization
     log = IterLog("VUMPS")
+    ϵ::Float64 = calc_galerkin(ψ, H, ψ, envs)
+    scheduler = Defaults.scheduler[]
+    temp_ACs = similar.(ψ.AC)
+    alg_environments = updatetol(alg.alg_environments, 0, ϵ)
+    recalculate!(envs, ψ, H; alg_environments.tol)
 
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs))
@@ -34,11 +37,11 @@ function leading_boundary(ψ::MultilineMPS, H, alg::VUMPS, envs=environments(ψ,
             ψ = MultilineMPS(temp_ACs, ψ.C[:, end]; alg_gauge.tol, alg_gauge.maxiter)
 
             alg_environments = updatetol(alg.alg_environments, iter, ϵ)
-            recalculate!(envs, ψ; alg_environments.tol)
+            recalculate!(envs, ψ, H, ψ; alg_environments.tol)
 
             ψ, envs = alg.finalize(iter, ψ, H, envs)::Tuple{typeof(ψ),typeof(envs)}
 
-            ϵ = calc_galerkin(ψ, envs)
+            ϵ = calc_galerkin(ψ, H, ψ, envs)
 
             if ϵ <= alg.tol
                 @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index 09fe6c8f3..715da66ae 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -57,7 +57,7 @@ function timestep(ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP,
         ψ′ = InfiniteMPS(ψ.C[0], temp_ACs; tol=alg.tolgauge, maxiter=alg.gaugemaxiter)
     end
 
-    recalculate!(envs, ψ′)
+    recalculate!(envs, ψ′, H)
     return ψ′, envs
 end
 
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 13214b348..b36276819 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -26,38 +26,36 @@ end
 # end
 
 """
-    calc_galerkin(state, envs)
+    calc_galerkin(above, operator, below, envs)
+    calc_galerkin(pos, above, operator, below, envs)
 
-Calculate the galerkin error.
+Calculate the Galerkin error, which is the error between the solution of the original problem, and the solution of the problem projected on the tangent space.
+Concretely, this is the overlap of the current state with the single-site derivative, projected onto the nullspace of the current state:
+
+```math
+\\epsilon = |VL * (VL' * \\frac{above}{\\partial AC_{pos}})|
+```
 """
-function calc_galerkin(state::Union{InfiniteMPS,FiniteMPS,WindowMPS}, loc, envs)::Float64
-    AC´ = ∂∂AC(loc, state, envs.operator, envs) * state.AC[loc]
+function calc_galerkin(pos::Int, above::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator,
+                       below, envs)
+    AC´ = ∂∂AC(pos, above, operator, envs) * above.AC[pos]
     normalize!(AC´)
-    out = add!(AC´, state.AL[loc] * state.AL[loc]' * AC´, -1)
+    out = add!(AC´, below.AL[pos] * below.AL[pos]' * AC´, -1)
     return norm(out)
 end
-function calc_galerkin(state::Union{InfiniteMPS,FiniteMPS,WindowMPS}, envs)::Float64
-    return maximum([calc_galerkin(state, loc, envs) for loc in 1:length(state)])
+function calc_galerkin(pos::CartesianIndex{2}, above::MultilineMPS, operator::MultilineMPO,
+                       below::MultilineMPS, envs::MultilineEnvironments)
+    row, col = pos.I
+    return calc_galerkin(col, above[row], operator[row], below[row + 1], envs[row])
 end
-function calc_galerkin(state::MultilineMPS, envs::InfiniteMPOEnvironments)::Float64
-    above = isnothing(envs.above) ? state : envs.above
-
-    εs = zeros(Float64, size(state, 1), size(state, 2))
-    for (row, col) in product(1:size(state, 1), 1:size(state, 2))
-        AC´ = ∂∂AC(row, col, state, envs.operator, envs) * above.AC[row, col]
-        normalize!(AC´)
-        out = add!(AC´, state.AL[row + 1, col] * state.AL[row + 1, col]' * AC´, -1)
-        εs[row, col] = norm(out)
-    end
-
-    return maximum(εs[:])
+function calc_galerkin(above::Union{InfiniteMPS,FiniteMPS,WindowMPS}, operator,
+                       below, envs)
+    return maximum(pos -> calc_galerkin(pos, above, operator, below, envs), 1:length(above))
 end
-function calc_galerkin(state::InfiniteMPS, site::Int,
-                       envs::InfiniteMPOHamiltonianEnvironments)
-    AC´ = ∂∂AC(site, state, envs.operator, envs) * state.AC[site]
-    normalize!(AC´)
-    out = add!(AC´, state.AL[site] * (state.AL[site]' * AC´), -1)
-    return norm(out)
+function calc_galerkin(above::MultilineMPS, operator::MultilineMPO, below::MultilineMPS,
+                       envs::MultilineEnvironments)
+    return maximum(pos -> calc_galerkin(pos, above, operator, below, envs),
+                   CartesianIndices(size(above)))
 end
 
 """
@@ -171,11 +169,8 @@ function variance end
 function variance(state::InfiniteMPS, H::InfiniteMPOHamiltonian,
                   envs=environments(state, H))
     e_local = map(1:length(state)) do i
-        return @plansor state.AC[i][3 7; 5] *
-                        leftenv(envs, i, state)[1 2; 3] *
-                        H[i][:, :, :, end][2 4; 7 8] *
-                        rightenv(envs, i, state)[end][5 8; 6] *
-                        conj(state.AC[i][1 4; 6])
+        return contract_mpo_expval(state.AC[i], envs.GLs[i], H[i][:, :, :, end],
+                                   envs.GRs[i][end])
     end
     lattice = physicalspace.(Ref(state), 1:length(state))
     H_renormalized = InfiniteMPOHamiltonian(lattice,
@@ -202,11 +197,9 @@ function variance(state::InfiniteQP, H::InfiniteMPOHamiltonian, envs=environment
     gs = state.left_gs
 
     e_local = map(1:length(state)) do i
-        return @plansor gs.AC[i][3 7; 5] *
-                        leftenv(envs.leftenvs, i, gs)[1 2; 3] *
-                        H[i][:, :, :, end][2 4; 7 8] *
-                        rightenv(envs.rightenvs, i, gs)[end][5 8; 6] *
-                        conj(gs.AC[i][1 4; 6])
+        GL = leftenv(envs, i, gs)
+        GR = rightenv(envs, i, gs)
+        return contract_mpo_expval(gs.AC[i], GL, H[i][:, :, :, end], GR[end])
     end
     lattice = physicalspace.(Ref(gs), 1:length(state))
     H_regularized = H - InfiniteMPOHamiltonian(lattice,
@@ -264,10 +257,10 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O},
         F_left = i == 1 ? cup : F_right
         F_right = i == L ? cup :
                   isomorphism(ST, V_right ← V_wrap' * right_virtualspace(mpo, i))
-        @plansor contractcheck = true output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
-                                                                τ[-3 1; 4 3] *
-                                                                mpo[i][2 -2; 3 5] *
-                                                                conj(F_right[-4; 4 5])
+        @plansor output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
+                                           τ[-3 1; 4 3] *
+                                           mpo[i][2 -2; 3 5] *
+                                           conj(F_right[-4; 4 5])
     end
 
     mpo isa SparseMPO && dropzeros!.(output) # the above process fills sparse mpos with zeros.
diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl
index f7ce69df0..c26f0d76a 100644
--- a/src/environments/abstract_envs.jl
+++ b/src/environments/abstract_envs.jl
@@ -13,7 +13,6 @@ Base.unlock(envs::AbstractMPSEnvironments) = unlock(envs.lock);
 
 # Allocating tensors
 # ------------------
-
 function allocate_GL(bra::AbstractMPS, mpo::AbstractMPO, ket::AbstractMPS, i::Int)
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
     V = left_virtualspace(bra, i) ⊗ left_virtualspace(mpo, i)' ←
@@ -62,46 +61,30 @@ function allocate_GBR(bra::QP, mpo::AbstractMPO, ket::QP, i::Int)
     return TT(undef, V)
 end
 
-# Abstract Infinite Environments
-# ------------------------------
+# Environment algorithms
+# ----------------------
 """
-    AbstractInfiniteEnvironments <: AbstractEnvironments
+    environment_alg(above, operator, below; kwargs...)
 
-Abstract supertype for infinite environment managers.
+Determine an appropriate algorithm for computing the environments, based on the given `kwargs...`.
 """
-abstract type AbstractInfiniteEnvironments <: AbstractMPSEnvironments end
-
-leftenv(envs, pos::CartesianIndex, state) = leftenv(envs, Tuple(pos)..., state)
-rightenv(envs, pos::CartesianIndex, state) = rightenv(envs, Tuple(pos)..., state)
-
-# recalculate logic
-# -----------------
-function check_recalculate!(envs::AbstractInfiniteEnvironments, state)
-    # check if dependency got updated - cheap test to avoid having to lock
-    if !check_dependency(envs, state)
-        # acquire lock and check again (might have updated while waiting)
-        lock(envs) do
-            return check_dependency(envs, state) || recalculate!(envs, state)
-        end
-    end
-    return envs
+function environment_alg(::Union{InfiniteMPS,MultilineMPS},
+                         ::Union{InfiniteMPO,MultilineMPO},
+                         ::Union{InfiniteMPS,MultilineMPS};
+                         tol=Defaults.tol, maxiter=Defaults.maxiter,
+                         krylovdim=Defaults.krylovdim, verbosity=Defaults.VERBOSE_NONE,
+                         eager=true)
+    return Arnoldi(; tol, maxiter, krylovdim, verbosity, eager)
 end
-
-function recalculate!(envs::AbstractInfiniteEnvironments, state; tol=envs.solver.tol)
-    # check if the virtual spaces have changed and reallocate if necessary
-    if !issamespace(envs, state)
-        envs.leftenvs, envs.rightenvs = initialize_environments(state, envs.operator)
-    end
-
-    solver = envs.solver
-    envs.solver = solver.tol == tol ? solver : Accessors.@set solver.tol = tol
-    envs.dependency = state
-
-    @sync begin
-        Threads.@spawn compute_leftenv!(envs)
-        Threads.@spawn compute_rightenv!(envs)
-    end
-    normalize!(envs)
-
-    return envs
+function environment_alg(above, ::InfiniteMPOHamiltonian, below;
+                         tol=Defaults.tol, maxiter=Defaults.maxiter,
+                         krylovdim=Defaults.krylovdim, verbosity=Defaults.VERBOSE_NONE)
+    return GMRES(; tol, maxiter, krylovdim, verbosity)
+end
+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)
+    return GMRES(; tol, maxiter, krylovdim, verbosity)
 end
diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index bf2c2d4fa..2fd437970 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -12,8 +12,8 @@ struct FiniteEnvironments{A,B,C,D} <: AbstractMPSEnvironments
     ldependencies::Vector{C} #the data we used to calculate leftenvs/rightenvs
     rdependencies::Vector{C}
 
-    leftenvs::Vector{D}
-    rightenvs::Vector{D}
+    GLs::Vector{D}
+    GRs::Vector{D}
 end
 
 function environments(below, (operator, above)::Tuple, args...; kwargs...)
@@ -53,8 +53,8 @@ function environments(below::WindowMPS, O::Union{InfiniteMPOHamiltonian,Infinite
                       above=nothing;
                       lenvs=environments(below.left_gs, O),
                       renvs=environments(below.right_gs, O))
-    leftstart = copy(leftenv(lenvs, 1, below.left_gs))
-    rightstart = copy(rightenv(renvs, length(below), below.right_gs))
+    leftstart = copy(lenvs.GLs[1])
+    rightstart = copy(renvs.GRs[end])
 
     return environments(below, O, above, leftstart, rightstart)
 end
@@ -93,13 +93,13 @@ function rightenv(ca::FiniteEnvironments, ind, state)
         #we need to recalculate
         for j in a:-1:(ind + 1)
             above = isnothing(ca.above) ? state.AR[j] : ca.above.AR[j]
-            ca.rightenvs[j] = TransferMatrix(above, ca.operator[j], state.AR[j]) *
-                              ca.rightenvs[j + 1]
+            ca.GRs[j] = TransferMatrix(above, ca.operator[j], state.AR[j]) *
+                        ca.GRs[j + 1]
             ca.rdependencies[j] = state.AR[j]
         end
     end
 
-    return ca.rightenvs[ind + 1]
+    return ca.GRs[ind + 1]
 end
 
 function leftenv(ca::FiniteEnvironments, ind, state)
@@ -109,11 +109,11 @@ function leftenv(ca::FiniteEnvironments, ind, state)
         #we need to recalculate
         for j in a:(ind - 1)
             above = isnothing(ca.above) ? state.AL[j] : ca.above.AL[j]
-            ca.leftenvs[j + 1] = ca.leftenvs[j] *
-                                 TransferMatrix(above, ca.operator[j], state.AL[j])
+            ca.GLs[j + 1] = ca.GLs[j] *
+                            TransferMatrix(above, ca.operator[j], state.AL[j])
             ca.ldependencies[j] = state.AL[j]
         end
     end
 
-    return ca.leftenvs[ind]
+    return ca.GLs[ind]
 end
diff --git a/src/environments/idmrg_envs.jl b/src/environments/idmrg_envs.jl
deleted file mode 100644
index 9a3682832..000000000
--- a/src/environments/idmrg_envs.jl
+++ /dev/null
@@ -1,55 +0,0 @@
-#=
-Idmrg environments are only to be used internally.
-They have to be updated manually, without any kind of checks
-=#
-"""
-    IDMRGEnvironments{O,V} <: AbstractMPSEnvironments
-
-Environment manager for IDMRG
-"""
-struct IDMRGEnvironments{O,V} <: AbstractMPSEnvironments
-    operator::O
-    leftenvs::PeriodicMatrix{V}
-    rightenvs::PeriodicMatrix{V}
-end
-
-function IDMRGEnvironments(ψ::InfiniteMPS, envs::InfiniteMPOHamiltonianEnvironments)
-    check_recalculate!(envs, ψ)
-    L = length(ψ)
-    leftenvs = PeriodicMatrix(reshape(deepcopy(envs.leftenvs), (1, L)))
-    rightenvs = PeriodicMatrix(reshape(deepcopy(envs.rightenvs), (1, L)))
-    return IDMRGEnvironments(envs.operator, leftenvs, rightenvs)
-end
-function IDMRGEnvironments(ψ::Union{InfiniteMPS,MultilineMPS},
-                           envs::InfiniteMPOEnvironments)
-    check_recalculate!(envs, ψ)
-    return IDMRGEnvironments(envs.operator, deepcopy(envs.leftenvs),
-                             deepcopy(envs.rightenvs))
-end
-
-leftenv(envs::IDMRGEnvironments, site::Int) = envs.leftenvs[site]
-leftenv(envs::IDMRGEnvironments, site::Int, ::InfiniteMPS) = envs.leftenvs[site]
-leftenv(envs::IDMRGEnvironments, row::Int, col::Int) = envs.leftenvs[row, col]
-setleftenv!(envs::IDMRGEnvironments, site::Int, GL) = (envs.leftenvs[site] = GL; envs)
-function setleftenv!(envs::IDMRGEnvironments, row::Int, col::Int, GL)
-    envs.leftenvs[row, col] = GL
-    return envs
-end
-
-rightenv(envs::IDMRGEnvironments, site::Int) = envs.rightenvs[site]
-rightenv(envs::IDMRGEnvironments, site::Int, ::InfiniteMPS) = envs.rightenvs[site]
-rightenv(envs::IDMRGEnvironments, row::Int, col::Int) = envs.rightenvs[row, col]
-setrightenv!(envs::IDMRGEnvironments, site::Int, GR) = (envs.rightenvs[site] = GR; envs)
-function setrightenv!(envs::IDMRGEnvironments, row::Int, col::Int, GR)
-    envs.rightenvs[row, col] = GR
-    return envs
-end
-
-function update_leftenv!(envs::IDMRGEnvironments, state, O, site::Int)
-    T = TransferMatrix(state.AL[site - 1], O[site - 1], state.AL[site - 1])
-    return setleftenv!(envs, site, leftenv(envs, site - 1) * T)
-end
-function update_rightenv!(envs::IDMRGEnvironments, state, O, site::Int)
-    T = TransferMatrix(state.AR[site + 1], O[site + 1], state.AR[site + 1])
-    return setrightenv!(envs, site, T * rightenv(envs, site + 1))
-end
diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl
new file mode 100644
index 000000000..f0b0de501
--- /dev/null
+++ b/src/environments/infinite_envs.jl
@@ -0,0 +1,295 @@
+"""
+    InfiniteEnvironments <: AbstractMPSEnvironments
+
+Environments for an infinite MPS-MPO-MPS combination. These solve the corresponding fixedpoint equations:
+```math
+GLs[i] * T_LL[i] = λ GLs[i + 1]
+T_RR[i] * GRs[i] = λ GRs[i - 1]
+```
+where `T_LL` and `T_RR` are the (regularized) transfer matrix operators on a give site for `AL-O-AL` and `AR-O-AR` respectively.
+"""
+struct InfiniteEnvironments{V<:GenericMPSTensor} <: AbstractMPSEnvironments
+    GLs::PeriodicVector{V}
+    GRs::PeriodicVector{V}
+end
+
+Base.length(envs::InfiniteEnvironments) = length(envs.GLs)
+
+leftenv(envs::InfiniteEnvironments, site::Int, state) = envs.GLs[site]
+rightenv(envs::InfiniteEnvironments, site::Int, state) = envs.GRs[site]
+
+function environments(above::InfiniteMPS,
+                      operator::Union{InfiniteMPO,InfiniteMPOHamiltonian},
+                      below::InfiniteMPS=above;
+                      kwargs...)
+    GLs, GRs = initialize_environments(above, operator, below)
+    envs = InfiniteEnvironments(GLs, GRs)
+    return recalculate!(envs, above, operator, below; kwargs...)
+end
+
+function issamespace(envs::InfiniteEnvironments, above::InfiniteMPS,
+                     operator::Union{InfiniteMPO,InfiniteMPOHamiltonian},
+                     below::InfiniteMPS)
+    L = check_length(above, operator, below)
+    for i in 1:L
+        space(envs.GLs[i]) ==
+        (left_virtualspace(below, i) ⊗ left_virtualspace(operator, i)' ←
+         left_virtualspace(above, i)) || return false
+        space(envs.GRs[i]) ==
+        (right_virtualspace(above, i) ⊗ right_virtualspace(operator, i) ←
+         right_virtualspace(below, i)) || return false
+    end
+    return true
+end
+
+function recalculate!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                      operator::Union{InfiniteMPO,InfiniteMPOHamiltonian},
+                      below::InfiniteMPS=above; kwargs...)
+    if !issamespace(envs, above, operator, below)
+        # TODO: in-place initialization?
+        GLs, GRs = initialize_environments(above, operator, below)
+        copy!(envs.GLs, GLs)
+        copy!(envs.GRs, GRs)
+    end
+
+    alg = environment_alg(above, operator, below; kwargs...)
+
+    @sync begin
+        @spawn compute_leftenvs!(envs, above, operator, below, alg)
+        @spawn compute_rightenvs!(envs, above, operator, below, alg)
+    end
+    normalize!(envs, above, operator, below)
+
+    return envs
+end
+
+# InfiniteMPO environments
+# ------------------------
+function initialize_environments(above::InfiniteMPS, operator::InfiniteMPO,
+                                 below::InfiniteMPS=above)
+    L = check_length(above, operator, below)
+    GLs = PeriodicVector([randomize!(allocate_GL(above, operator, below, i)) for i in 1:L])
+    GRs = PeriodicVector([randomize!(allocate_GR(above, operator, below, i)) for i in 1:L])
+    return GLs, GRs
+end
+
+function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                           operator::InfiniteMPO, below::InfiniteMPS, alg)
+    # compute eigenvector
+    T = TransferMatrix(above.AL, operator, below.AL)
+    λ, envs.GLs[1] = fixedpoint(flip(T), envs.GLs[1], :LM, alg)
+    # push through unitcell
+    for i in 2:length(operator)
+        envs.GLs[i] = envs.GLs[i - 1] *
+                      TransferMatrix(above.AL[i - 1], operator[i - 1], below.AL[i - 1])
+    end
+    return λ, envs
+end
+
+function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                            operator::InfiniteMPO, below::InfiniteMPS, alg)
+    # compute eigenvector
+    T = TransferMatrix(above.AR, operator, below.AR)
+    λ, envs.GRs[end] = fixedpoint(T, envs.GRs[end], :LM, alg)
+    # 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]) * envs.GRs[i + 1]
+    end
+    return λ, envs
+end
+
+function TensorKit.normalize!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                              operator::InfiniteMPO, below::InfiniteMPS)
+    for i in 1:length(operator)
+        λ = dot(below.C[i], MPO_∂∂C(envs.GLs[i + 1], envs.GRs[i]) * above.C[i])
+        scale!(envs.GLs[i + 1], inv(λ))
+    end
+    return envs
+end
+
+# InfiniteMPOHamiltonian environments
+# -----------------------------------
+function initialize_environments(above::InfiniteMPS, operator::InfiniteMPOHamiltonian,
+                                 below::InfiniteMPS=above)
+    L = check_length(above, operator, below)
+    GLs = PeriodicVector([allocate_GL(above, operator, below, i) for i in 1:L])
+    GRs = PeriodicVector([allocate_GR(above, operator, below, i) for i in 1:L])
+
+    # GL = (1, 0, 0)
+    GL = first(GLs)
+    for i in 1:length(GL)
+        if i == 1
+            GL[i] = isomorphism(storagetype(GL), space(GL[i]))
+        else
+            fill!(GL[i], zero(scalartype(GL)))
+        end
+    end
+
+    # GR = (0, 0, 1)^T
+    GR = last(GRs)
+    for i in 1:length(GR)
+        if i == length(GR)
+            GR[i] = isomorphism(storagetype(GR), space(GR[i]))
+        else
+            fill!(GR[i], zero(scalartype(GR)))
+        end
+    end
+
+    return GLs, GRs
+end
+
+function compute_leftenvs!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                           operator::InfiniteMPOHamiltonian, below::InfiniteMPS, alg)
+    L = check_length(above, operator, below)
+    GLs = envs.GLs
+    vsize = length(first(GLs))
+
+    @assert above === below "not implemented"
+
+    ρ_left = l_LL(above)
+    ρ_right = r_LL(above)
+
+    # the start element
+    # TODO: check if this is necessary
+    # leftutil = similar(above.AL[1], space(GL[1], 2)[1])
+    # fill_data!(leftutil, one)
+    # @plansor GL[1][1][-1 -2; -3] = ρ_left[-1; -3] * leftutil[-2]
+
+    (L > 1) && left_cyclethrough!(1, GLs, above, operator, below)
+
+    for i in 2:vsize
+        prev = copy(GLs[1][i])
+        zerovector!(GLs[1][i])
+        left_cyclethrough!(i, GLs, above, operator, below)
+
+        if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization
+            T = regularize(TransferMatrix(above.AL, below.AL), ρ_left, ρ_right)
+            GLs[1][i], convhist = linsolve(flip(T), GLs[1][i], prev, alg, 1, -1)
+            convhist.converged == 0 &&
+                @warn "GL$i failed to converge: normres = $(convhist.normres)"
+
+            (L > 1) && left_cyclethrough!(i, GLs, above, operator, below)
+
+            # go through the unitcell, again subtracting fixpoints
+            for site in 1:L
+                @plansor GLs[site][i][-1 -2; -3] -= GLs[site][i][1 -2; 2] *
+                                                    r_LL(above, site - 1)[2; 1] *
+                                                    l_LL(above, site)[-1; -3]
+            end
+
+        else
+            if !isemptylevel(operator, i)
+                diag = map(h -> h[i, 1, 1, i], operator[:])
+                T = TransferMatrix(above.AL, diag, below.AL)
+                GLs[1][i], convhist = linsolve(flip(T), GLs[1][i], prev, alg, 1, -1)
+                convhist.converged == 0 &&
+                    @warn "GL$i failed to converge: normres = $(convhist.normres)"
+            end
+            (L > 1) && left_cyclethrough!(i, GLs, above, operator, below)
+        end
+    end
+
+    return GLs
+end
+
+function left_cyclethrough!(index::Int, GL, above::InfiniteMPS, H::InfiniteMPOHamiltonian,
+                            below::InfiniteMPS=above)
+    # TODO: efficient transfer matrix slicing for large unitcells
+    leftinds = 1:index
+    for site in eachindex(GL)
+        GL[site + 1][index] = GL[site][leftinds] * TransferMatrix(above.AL[site],
+                                                                  H[site][leftinds, 1, 1, index],
+                                                                  below.AL[site])
+    end
+    return GL
+end
+
+function compute_rightenvs!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                            operator::InfiniteMPOHamiltonian, below::InfiniteMPS, alg)
+    L = check_length(above, operator, below)
+    GRs = envs.GRs
+    vsize = length(last(GRs))
+
+    @assert above === below "not implemented"
+
+    ρ_left = l_RR(above)
+    ρ_right = r_RR(above)
+
+    # the start element
+    # TODO: check if this is necessary
+    # rightutil = similar(state.AL[1], space(GR[end], 2)[end])
+    # fill_data!(rightutil, one)
+    # @plansor GR[end][end][-1 -2; -3] = r_RR(state)[-1; -3] * rightutil[-2]
+
+    (L > 1) && right_cyclethrough!(vsize, GRs, above, operator, below) # populate other sites
+
+    for i in (vsize - 1):-1:1
+        prev = copy(GRs[end][i])
+        zerovector!(GRs[end][i])
+        right_cyclethrough!(i, GRs, above, operator, below)
+
+        if isidentitylevel(operator, i) # identity matrices; do the hacky renormalization
+            # subtract fixpoints
+            T = regularize(TransferMatrix(above.AR, below.AR), ρ_left, ρ_right)
+            GRs[end][i], convhist = linsolve(T, GRs[end][i], prev, alg, 1, -1)
+            convhist.converged == 0 &&
+                @warn "GR$i failed to converge: normres = $(convhist.normres)"
+
+            L > 1 && right_cyclethrough!(i, GRs, above, operator, below)
+
+            # go through the unitcell, again subtracting fixpoints
+            for site in 1:L
+                @plansor GRs[site][i][-1 -2; -3] -= GRs[site][i][1 -2; 2] *
+                                                    l_RR(above, site + 1)[2; 1] *
+                                                    r_RR(above, site)[-1; -3]
+            end
+        else
+            if !isemptylevel(operator, i)
+                diag = map(b -> b[i, 1, 1, i], operator[:])
+                T = TransferMatrix(above.AR, diag, below.AR)
+                GRs[end][i], convhist = linsolve(T, GRs[end][i], prev, alg, 1, -1)
+                convhist.converged == 0 &&
+                    @warn "GR$i failed to converge: normres = $(convhist.normres)"
+            end
+
+            (L > 1) && right_cyclethrough!(i, GRs, above, operator, below)
+        end
+    end
+
+    return GRs
+end
+
+function right_cyclethrough!(index::Int, GR, above::InfiniteMPS,
+                             operator::InfiniteMPOHamiltonian,
+                             below::InfiniteMPS)
+    # TODO: efficient transfer matrix slicing for large unitcells
+    for site in reverse(eachindex(GR))
+        rightinds = index:length(GR[site])
+        GR[site - 1][index] = TransferMatrix(above.AR[site],
+                                             operator[site][index, 1, 1, rightinds],
+                                             below.AR[site]) * GR[site][rightinds]
+    end
+    return GR
+end
+
+# no normalization necessary -- for consistant interface
+function TensorKit.normalize!(envs::InfiniteEnvironments, above::InfiniteMPS,
+                              operator::InfiniteMPOHamiltonian, below::InfiniteMPS)
+    return envs
+end
+
+# Transfer operations
+# -------------------
+
+function transfer_leftenv!(envs::InfiniteEnvironments, above, operator, below, site::Int)
+    T = TransferMatrix(above.AL[site - 1], operator[site - 1], below.AL[site - 1])
+    envs.GLs[site] = envs.GLs[site - 1] * T
+    return envs
+end
+
+function transfer_rightenv!(envs::InfiniteEnvironments, above, operator, below, site::Int)
+    T = TransferMatrix(above.AR[site + 1], operator[site + 1], below.AR[site + 1])
+    envs.GRs[site] = T * envs.GRs[site + 1]
+    return envs
+end
diff --git a/src/environments/infinitempo_envs.jl b/src/environments/infinitempo_envs.jl
deleted file mode 100644
index 2ea1d5312..000000000
--- a/src/environments/infinitempo_envs.jl
+++ /dev/null
@@ -1,186 +0,0 @@
-"""
-    InfiniteMPOEnvironments{O<:MultilineMPO,V,S<:MultilineMPS,A} <: AbstractMPSEnvironments
-
-Environment manager for `InfiniteMPO` and its `Multiline` version.
-"""
-mutable struct InfiniteMPOEnvironments{O,V,S<:MultilineMPS,A} <:
-               AbstractInfiniteEnvironments
-    above::Union{S,Nothing}
-    operator::O
-
-    dependency::S
-    solver::A
-
-    leftenvs::PeriodicMatrix{V}
-    rightenvs::PeriodicMatrix{V}
-
-    lock::ReentrantLock
-end
-function InfiniteMPOEnvironments(bra, O, ket, solver, GL, GR)
-    return InfiniteMPOEnvironments(bra, O, ket, solver, GL, GR, ReentrantLock())
-end
-
-# Constructors
-# ------------
-function environments(state::InfiniteMPS, O::InfiniteMPO; kwargs...)
-    return environments(convert(MultilineMPS, state), convert(MultilineMPO, O); kwargs...)
-end
-function environments(below::InfiniteMPS,
-                      (mpo, above)::Tuple{<:InfiniteMPO,<:InfiniteMPS}; kwargs...)
-    return environments(convert(MultilineMPS, below),
-                        (convert(MultilineMPO, mpo), convert(MultilineMPS, above));
-                        kwargs...)
-end
-
-function environments(state::MultilineMPS, mpo::MultilineMPO; solver=Defaults.eigsolver)
-    GL, GR = initialize_environments(state, mpo, state)
-    envs = InfiniteMPOEnvironments(nothing, mpo, state, solver, GL, GR)
-    return recalculate!(envs, state)
-end
-
-function environments(below::MultilineMPS,
-                      (mpo, above)::Tuple{<:MultilineMPO,<:MultilineMPS};
-                      solver=Defaults.eigsolver)
-    GL, GR = initialize_environments(above, mpo, below)
-    envs = InfiniteMPOEnvironments(above, mpo, below, solver, GL, GR)
-    return recalculate!(envs, below)
-end
-
-function initialize_environments(ket::MultilineMPS, operator::MultilineMPO,
-                                 bra::MultilineMPS=ket)
-    # allocate
-    GL = PeriodicArray([allocate_GL(bra[row], operator[row], ket[row], col)
-                        for row in 1:size(ket, 1), col in 1:size(ket, 2)])
-    GR = PeriodicArray([allocate_GR(bra[row], operator[row], ket[row], col)
-                        for row in 1:size(ket, 1), col in 1:size(ket, 2)])
-
-    # initialize: randomize
-    foreach(randomize!, GL)
-    foreach(randomize!, GR)
-
-    return GL, GR
-end
-
-# Getter/Setters
-# --------------
-function leftenv(envs::InfiniteMPOEnvironments, pos::Int, state::InfiniteMPS)
-    check_recalculate!(envs, state)
-    return envs.leftenvs[1, pos]
-end
-function leftenv(envs::InfiniteMPOEnvironments, pos::Int, state::MultilineMPS)
-    check_recalculate!(envs, state)
-    return envs.leftenvs[:, pos]
-end
-function leftenv(envs::InfiniteMPOEnvironments, row::Int, col::Int, state)
-    check_recalculate!(envs, state)
-    return envs.leftenvs[row, col]
-end
-
-function rightenv(envs::InfiniteMPOEnvironments, pos::Int, state::InfiniteMPS)
-    check_recalculate!(envs, state)
-    return envs.rightenvs[1, pos]
-end
-function rightenv(envs::InfiniteMPOEnvironments, pos::Int, state::MultilineMPS)
-    check_recalculate!(envs, state)
-    return envs.rightenvs[:, pos]
-end
-function rightenv(envs::InfiniteMPOEnvironments, row::Int, col::Int, state)
-    check_recalculate!(envs, state)
-    return envs.rightenvs[row, col]
-end
-
-# Utility
-# -------
-function check_dependency(envs::InfiniteMPOEnvironments, state::MultilineMPS)
-    return all(x -> ===(x...), zip(envs.dependency, state))
-end
-
-function issamespace(envs::InfiniteMPOEnvironments, state::MultilineMPS)
-    for row in 1:size(state, 1)
-        newstate = state[row]
-        oldstate = envs.dependency[row]
-        for col in 1:size(state, 2)
-            if left_virtualspace(oldstate, col) != left_virtualspace(newstate, col)
-                return false
-            end
-            if right_virtualspace(oldstate, col) != right_virtualspace(newstate, col)
-                return false
-            end
-        end
-    end
-    return true
-end
-
-# Calculation
-# -----------
-function recalculate!(envs::InfiniteMPOEnvironments, nstate::InfiniteMPS; kwargs...)
-    return recalculate!(envs, convert(MultilineMPS, nstate); kwargs...)
-end
-
-function compute_leftenv!(envs::InfiniteMPOEnvironments)
-    below = envs.dependency
-    above = something(envs.above, below)
-    mpo = envs.operator
-
-    # sanity check
-    numrows, numcols = size(above)
-    @assert size(above) == size(mpo)
-    @assert size(below) == size(mpo)
-
-    @threads for row in 1:numrows
-        T = TransferMatrix(above[row].AL, mpo[row, :], below[row + 1].AL)
-        _, envs.leftenvs[row, 1] = fixedpoint(flip(T), envs.leftenvs[row, 1], :LM,
-                                              envs.solver)
-        # compute rest of unitcell
-        for col in 2:numcols
-            envs.leftenvs[row, col] = envs.leftenvs[row, col - 1] *
-                                      TransferMatrix(above[row].AL[col - 1],
-                                                     mpo[row, col - 1],
-                                                     below[row + 1].AL[col - 1])
-        end
-    end
-
-    return envs
-end
-
-function compute_rightenv!(envs::InfiniteMPOEnvironments)
-    below = envs.dependency
-    above = something(envs.above, below)
-    mpo = envs.operator
-
-    # sanity check
-    numrows, numcols = size(above)
-    @assert size(above) == size(mpo)
-    @assert size(below) == size(mpo)
-
-    @threads for row in 1:numrows
-        T = TransferMatrix(above[row].AR, mpo[row, :], below[row + 1].AR)
-        _, envs.rightenvs[row, end] = fixedpoint(T, envs.rightenvs[row, end], :LM,
-                                                 envs.solver)
-        # compute rest of unitcell
-        for col in (numcols - 1):-1:1
-            envs.rightenvs[row, col] = TransferMatrix(above[row].AR[col + 1],
-                                                      mpo[row, col + 1],
-                                                      below[row].AR[col + 1]) *
-                                       envs.rightenvs[row, col + 1]
-        end
-    end
-
-    return envs
-end
-
-function TensorKit.normalize!(envs::InfiniteMPOEnvironments)
-    below = envs.dependency
-    above = something(envs.above, below)
-
-    for row in 1:size(below, 1)
-        # fix normalization
-        CRs_top, CRs_bot = above[row].C, below[row + 1].C
-        for col in 1:size(below, 2)
-            λ = dot(CRs_bot[col],
-                    MPO_∂∂C(envs.leftenvs[row, col + 1], envs.rightenvs[row, col]) *
-                    CRs_top[col])
-            scale!(envs.leftenvs[row, col + 1], inv(λ))
-        end
-    end
-end
diff --git a/src/environments/infinitempohamiltonian_envs.jl b/src/environments/infinitempohamiltonian_envs.jl
deleted file mode 100644
index 62c789ade..000000000
--- a/src/environments/infinitempohamiltonian_envs.jl
+++ /dev/null
@@ -1,211 +0,0 @@
-"""
-    InfiniteMPOHamiltonianEnvironments{O<:InfiniteMPOHamiltonian,V,S<:InfiniteMPS,A} <: AbstractMPSEnvironments
-
-Environment manager for `InfiniteMPOHamiltonian`.
-"""
-mutable struct InfiniteMPOHamiltonianEnvironments{O,V,S,A} <: AbstractInfiniteEnvironments
-    operator::O
-    dependency::S
-
-    solver::A
-
-    leftenvs::PeriodicVector{V}
-    rightenvs::PeriodicVector{V}
-
-    lock::ReentrantLock
-end
-function InfiniteMPOHamiltonianEnvironments(operator, dependency, solver, leftenvs,
-                                            rightenvs)
-    return InfiniteMPOHamiltonianEnvironments(operator, dependency, solver, leftenvs,
-                                              rightenvs, ReentrantLock)
-end
-
-# Constructors
-# ------------
-function environments(state::InfiniteMPS, H::InfiniteMPOHamiltonian;
-                      solver=Defaults.linearsolver)
-    GL, GR = initialize_environments(state, H)
-    envs = InfiniteMPOHamiltonianEnvironments(H, state, solver, GL, GR, ReentrantLock())
-    return recalculate!(envs, state)
-end
-
-function initialize_environments(state::InfiniteMPS, H::InfiniteMPOHamiltonian)
-    # allocate
-    GL = PeriodicArray([allocate_GL(state, H, state, i) for i in 1:length(state)])
-    GR = PeriodicArray([allocate_GR(state, H, state, i) for i in 1:length(state)])
-
-    # initialize:
-    # GL = (1, 0, 0)
-    for i in 1:length(GL[1])
-        if i == 1
-            GL[1][i] = isomorphism(storagetype(eltype(GL)), space(GL[1][i]))
-        else
-            fill!(GL[1][i], zero(scalartype(GL)))
-        end
-    end
-    # GR = (0, 0, 1)^T
-    for i in 1:length(GR[end])
-        if i == length(GR[end])
-            GR[end][i] = isomorphism(storagetype(eltype(GR)), space(GR[end][i]))
-        else
-            fill!(GR[end][i], zero(scalartype(GR)))
-        end
-    end
-
-    return GL, GR
-end
-
-# Getter/Setters
-# --------------
-function leftenv(envs::InfiniteMPOHamiltonianEnvironments, pos::Int, ψ)
-    check_recalculate!(envs, ψ)
-    return envs.leftenvs[pos]
-end
-
-function rightenv(envs::InfiniteMPOHamiltonianEnvironments, pos::Int, ψ)
-    check_recalculate!(envs, ψ)
-    return envs.rightenvs[pos]
-end
-
-# Utility
-# -------
-function check_dependency(envs::InfiniteMPOHamiltonianEnvironments, state::InfiniteMPS)
-    return envs.dependency === state
-end
-
-function issamespace(envs::InfiniteMPOHamiltonianEnvironments, state::InfiniteMPS)
-    return all(1:length(state)) do i
-        return left_virtualspace(envs.dependency, i) == left_virtualspace(state, i)
-    end
-end
-
-# Calculation
-# -----------
-function compute_leftenv!(envs::InfiniteMPOHamiltonianEnvironments)
-    state = envs.dependency
-    H = envs.operator
-    GL = envs.leftenvs
-    solver = envs.solver
-    L = length(state)
-
-    # the start element
-    # TODO: check if this is necessary
-    leftutil = similar(state.AL[1], space(GL[1], 2)[1])
-    fill_data!(leftutil, one)
-    ρ_left = l_LL(state)
-    @plansor GL[1][1][-1 -2; -3] = ρ_left[-1; -3] * leftutil[-2]
-
-    (L > 1) && left_cyclethrough!(1, GL, H, state)
-
-    for i in 2:length(GL[1])
-        prev = copy(GL[1][i])
-        zerovector!(GL[1][i])
-        left_cyclethrough!(i, GL, H, state)
-
-        if isidentitylevel(H, i) # identity matrices; do the hacky renormalization
-            T = regularize(TransferMatrix(state.AL, state.AL), ρ_left, r_LL(state))
-            GL[1][i], convhist = linsolve(flip(T), GL[1][i], prev, solver, 1, -1)
-            convhist.converged == 0 &&
-                @warn "GL$i failed to converge: normres = $(convhist.normres)"
-
-            (L > 1) && left_cyclethrough!(i, GL, H, state)
-
-            # go through the unitcell, again subtracting fixpoints
-            for site in 1:L
-                @plansor GL[site][i][-1 -2; -3] -= GL[site][i][1 -2; 2] *
-                                                   r_LL(state, site - 1)[2; 1] *
-                                                   l_LL(state, site)[-1; -3]
-            end
-
-        else
-            if !isemptylevel(H, i)
-                diag = map(h -> h[i, 1, 1, i], H[:])
-                T = TransferMatrix(state.AL, diag, state.AL)
-                GL[1][i], convhist = linsolve(flip(T), GL[1][i], prev, solver, 1, -1)
-                convhist.converged == 0 &&
-                    @warn "GL$i failed to converge: normres = $(convhist.normres)"
-            end
-            (L > 1) && left_cyclethrough!(i, GL, H, state)
-        end
-    end
-
-    return GL
-end
-
-function left_cyclethrough!(index::Int, GL, H::InfiniteMPOHamiltonian, state)
-    # TODO: efficient transfer matrix slicing for large unitcells
-    for site in 1:length(GL)
-        leftinds = 1:index
-        GL[site + 1][index] = GL[site][leftinds] * TransferMatrix(state.AL[site],
-                                                                  H[site][leftinds, 1, 1, index],
-                                                                  state.AL[site])
-    end
-    return GL
-end
-
-function compute_rightenv!(envs::InfiniteMPOHamiltonianEnvironments)
-    GR = envs.rightenvs
-    H = envs.operator
-    state = envs.dependency
-    solver = envs.solver
-    L = length(state)
-
-    odim = length(GR[end])
-
-    # the start element
-    rightutil = similar(state.AL[1], space(GR[end], 2)[end])
-    fill_data!(rightutil, one)
-    @plansor GR[end][end][-1 -2; -3] = r_RR(state)[-1; -3] * rightutil[-2]
-
-    (L > 1) && right_cyclethrough!(odim, GR, H, state) # populate other sites
-
-    for i in (odim - 1):-1:1
-        prev = copy(GR[end][i])
-        zerovector!(GR[end][i])
-        right_cyclethrough!(i, GR, H, state)
-
-        if isidentitylevel(H, i) # identity matrices; do the hacky renormalization
-            # subtract fixpoints
-            T = regularize(TransferMatrix(state.AR, state.AR), l_RR(state), r_RR(state))
-            GR[end][i], convhist = linsolve(T, GR[end][i], prev, solver, 1, -1)
-            convhist.converged == 0 &&
-                @warn "GR$i failed to converge: normres = $(convhist.normres)"
-
-            L > 1 && right_cyclethrough!(i, GR, H, state)
-
-            # go through the unitcell, again subtracting fixpoints
-            for site in 1:L
-                @plansor GR[site][i][-1 -2; -3] -= GR[site][i][1 -2; 2] *
-                                                   l_RR(state, site + 1)[2; 1] *
-                                                   r_RR(state, site)[-1; -3]
-            end
-        else
-            if !isemptylevel(H, i)
-                diag = map(b -> b[i, 1, 1, i], H[:])
-                T = TransferMatrix(state.AR, diag, state.AR)
-                GR[end][i], convhist = linsolve(T, GR[end][i], prev, solver, 1, -1)
-                convhist.converged == 0 &&
-                    @warn "GR$i failed to converge: normres = $(convhist.normres)"
-            end
-
-            (L > 1) && right_cyclethrough!(i, GR, H, state)
-        end
-    end
-
-    return GR
-end
-
-function right_cyclethrough!(index::Int, GR, H::InfiniteMPOHamiltonian, state)
-    # TODO: efficient transfer matrix slicing for large unitcells
-    L = length(GR)
-    for site in reverse(1:L)
-        rightinds = index:length(GR[site])
-        GR[site - 1][index] = TransferMatrix(state.AR[site],
-                                             H[site][index, 1, 1, rightinds],
-                                             state.AR[site]) * GR[site][rightinds]
-    end
-    return GR
-end
-
-# no normalization here, but need this for consistent interface
-TensorKit.normalize!(envs::InfiniteMPOHamiltonianEnvironments) = envs
diff --git a/src/environments/multiline_envs.jl b/src/environments/multiline_envs.jl
new file mode 100644
index 000000000..de84fef0e
--- /dev/null
+++ b/src/environments/multiline_envs.jl
@@ -0,0 +1,65 @@
+const MultilineEnvironments{E<:AbstractMPSEnvironments} = Multiline{E}
+
+function environments(above::MultilineMPS, operator::MultilineMPO,
+                      below::MultilineMPS=above; kwargs...)
+    (rows = size(above, 1)) == size(operator, 1) == size(below, 1) ||
+        throw(ArgumentError("Incompatible sizes"))
+    envs = map(1:rows) do row
+        return environments(above[row], operator[row], below[row + 1]; kwargs...)
+    end
+    return Multiline(PeriodicVector(envs))
+end
+
+function recalculate!(envs::MultilineEnvironments, above::MultilineMPS,
+                      operator::MultilineMPO, below::MultilineMPS=above; kwargs...)
+    (rows = size(above, 1)) == size(operator, 1) == size(below, 1) ||
+        throw(ArgumentError("Incompatible sizes"))
+    @threads for row in 1:rows
+        recalculate!(envs[row], above[row], operator[row], below[row + 1]; kwargs...)
+    end
+    return envs
+end
+function recalculate!(envs::MultilineEnvironments, below, (operator, above)::Tuple;
+                      kwargs...)
+    return recalculate!(envs, above, operator, below; kwargs...)
+end
+
+function TensorKit.normalize!(envs::MultilineEnvironments, above, operator, below)
+    for row in 1:size(above, 1)
+        normalize!(envs[row], above[row], operator[row], below[row + 1])
+    end
+    return envs
+end
+function TensorKit.normalize!(envs::MultilineEnvironments, below, (operator, above))
+    for row in 1:size(above, 1)
+        normalize!(envs[row], above[row], operator[row], below[row + 1])
+    end
+    return envs
+end
+
+function leftenv(envs::MultilineEnvironments, col::Int, state)
+    return leftenv.(parent(envs), col, parent(state))
+end
+function rightenv(envs::MultilineEnvironments, col::Int, state)
+    return rightenv.(parent(envs), col, parent(state))
+end
+
+function transfer_leftenv!(envs::MultilineEnvironments, above, operator, below, site::Int)
+    for row in 1:size(above, 1)
+        transfer_leftenv!(envs[row], above[row], operator[row], below[row + 1], site)
+    end
+    return envs
+end
+function transfer_leftenv!(envs::MultilineEnvironments, below, (O, above)::Tuple, site::Int)
+    return transfer_leftenv!(envs, above, O, below, site)
+end
+function transfer_rightenv!(envs::MultilineEnvironments, above, operator, below, site::Int)
+    for row in 1:size(above, 1)
+        transfer_rightenv!(envs[row], above[row], operator[row], below[row + 1], site)
+    end
+    return envs
+end
+function transfer_rightenv!(envs::MultilineEnvironments, below, (O, above)::Tuple,
+                            site::Int)
+    return transfer_rightenv!(envs, above, O, below, site)
+end
diff --git a/src/environments/multiple_envs.jl b/src/environments/multiple_envs.jl
index 6eb2a75b3..97a6f70f0 100644
--- a/src/environments/multiple_envs.jl
+++ b/src/environments/multiple_envs.jl
@@ -1,5 +1,4 @@
-struct MultipleEnvironments{O,C} <: AbstractMPSEnvironments
-    operator::O
+struct MultipleEnvironments{C} <: AbstractMPSEnvironments
     envs::Vector{C}
 end
 
@@ -11,38 +10,15 @@ 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(st, H::LazySum)
-    return MultipleEnvironments(H, map(op -> environments(st, op), H.ops))
+function environments(state, H::LazySum; kwargs...)
+    return MultipleEnvironments(map(Base.Fix1(environments, state), H.ops))
 end
 
-function environments(st::Union{InfiniteMPS,MultilineMPS}, H::LazySum;
-                      solver=Defaults.linearsolver)
-    if !(solver isa Vector)
-        solver = repeat([solver], length(H))
-    end
-    return MultipleEnvironments(H,
-                                map((op, solv) -> environments(st, op; solver=solv),
-                                    H.ops, solver))
-end
-
-# TODO: fix this such that `T(...) isa T`
-function IDMRGEnvironments(ψ::Union{MultilineMPS,InfiniteMPS}, env::MultipleEnvironments)
-    envs = IDMRGEnvironments.(Ref(ψ), env.envs)
-    Hs = getproperty.(env.envs, :operator)
-    return MultipleEnvironments(LazySum(Hs), envs)
-end
-
-#broadcast vs map?
-# function environments(state, H::LinearCombination)
-#     return MultipleEnvironments(H, broadcast(o -> environments(state, o), H.opps))
-# end;
-
 function environments(st::WindowMPS,
                       H::LazySum;
                       lenvs=environments(st.left_gs, H),
                       renvs=environments(st.right_gs, H))
-    return MultipleEnvironments(H,
-                                map((op, sublenv, subrenv) -> environments(st, op;
+    return MultipleEnvironments(map((op, sublenv, subrenv) -> environments(st, op;
                                                                            lenvs=sublenv,
                                                                            renvs=subrenv),
                                     H.ops, lenvs, renvs))
@@ -52,18 +28,12 @@ end
 """
     Recalculate in-place each sub-env in MultipleEnvironments
 """
-function recalculate!(env::MultipleEnvironments, args...; kwargs...)
-    for subenv in env.envs
-        recalculate!(subenv, args...; kwargs...)
-    end
-    return env
-end
-
-function check_recalculate!(env::MultipleEnvironments, args...; kwargs...)
-    for subenv in env.envs
-        check_recalculate!(subenv, args...; kwargs...)
+function recalculate!(envs::MultipleEnvironments, above, operator::LazySum, below=above;
+                      kwargs...)
+    for (subenvs, subO) in zip(envs.envs, operator)
+        recalculate!(subenvs, above, subO, below; kwargs...)
     end
-    return env
+    return envs
 end
 
 #maybe this can be used to provide compatibility with existing code?
@@ -75,18 +45,18 @@ function Base.getproperty(envs::MultipleEnvironments, prop::Symbol)
     end
 end
 
-function update_rightenv!(envs::MultipleEnvironments{<:LazySum,<:IDMRGEnvironments}, st, H,
-                          pos::Int)
-    for (subH, subenv) in zip(H, envs.envs)
-        tm = TransferMatrix(st.AR[pos + 1], subH[pos + 1], st.AR[pos + 1])
-        setrightenv!(subenv, pos, tm * rightenv(subenv, pos + 1))
+function transfer_rightenv!(envs::MultipleEnvironments{<:InfiniteEnvironments},
+                            above, operator, below, pos::Int)
+    for (subH, subenv) in zip(operator, envs.envs)
+        transfer_rightenv!(subenv, above, subH, below, pos)
     end
+    return envs
 end
 
-function update_leftenv!(envs::MultipleEnvironments{<:LazySum,<:IDMRGEnvironments}, st, H,
-                         pos::Int)
-    for (subH, subenv) in zip(H, envs.envs)
-        tm = TransferMatrix(st.AL[pos - 1], subH[pos - 1], st.AL[pos - 1])
-        setleftenv!(subenv, pos, leftenv(subenv, pos - 1) * tm)
+function transfer_leftenv!(envs::MultipleEnvironments{<:InfiniteEnvironments},
+                           above, operator, below, pos::Int)
+    for (subH, subenv) in zip(operator, envs.envs)
+        transfer_leftenv!(subenv, above, subH, below, pos)
     end
+    return envs
 end
diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl
index baaf63169..b48fdd0e2 100644
--- a/src/environments/qp_envs.jl
+++ b/src/environments/qp_envs.jl
@@ -1,17 +1,14 @@
-#=
-nothing fancy - only used internally (and therefore cryptic) - stores some partially contracted things
-seperates out this bit of logic from effective_excitation_hamiltonian (now more readable)
-can also - potentially - partially reuse this in other algorithms
-=#
-struct QPEnv{A,B} <: AbstractMPSEnvironments
-    lBs::PeriodicArray{A,2}
-    rBs::PeriodicArray{A,2}
-
-    lenvs::B
-    renvs::B
-end
-
-struct QuasiparticleEnvironments{A,B} <: AbstractMPSEnvironments
+"""
+    InfiniteQPEnvironments <: AbstractMPSEnvironments
+
+Environments for an infinite QP-MPO-QP combination. These solve the corresponding fixedpoint equations:
+```math
+GLs[i] * T_BL[i] + GBLs[i] * T_RL[i] = GBLs[i + 1]
+T_BR[i] * GRs[i] + T_LR[i] * GBRs[i] = GBRs[i - 1]
+```
+where `T_BL`, `T_BR`, `T_RL` and `T_LR` are the (regularized) transfer matrix operators on a given site for `B-O-AL`, `B-O-AR`, `AR-O-AL` and `AL-O-AR` respectively.
+"""
+struct InfiniteQPEnvironments{A,B} <: AbstractMPSEnvironments
     leftBenvs::PeriodicVector{A}
     rightBenvs::PeriodicVector{A}
 
@@ -19,52 +16,41 @@ struct QuasiparticleEnvironments{A,B} <: AbstractMPSEnvironments
     rightenvs::B
 end
 
-function environments(exci::Union{InfiniteQP,Multiline{<:InfiniteQP}}, H;
-                      solver=Defaults.linearsolver)
-    # Explicitly define optional arguments as these depend on solver,
-    # which needs to come after these arguments.
-    lenvs = environments(exci.left_gs, H; solver=solver)
+Base.length(envs::InfiniteQPEnvironments) = length(envs.leftenvs)
 
-    return environments(exci, H, lenvs; solver=solver)
+function leftenv(envs::InfiniteQPEnvironments, site::Int, state)
+    return leftenv(envs.leftenvs, site, state)
+end
+function rightenv(envs::InfiniteQPEnvironments, site::Int, state)
+    return rightenv(envs.rightenvs, site, state)
 end
 
-function environments(exci::Union{InfiniteQP,Multiline{<:InfiniteQP}}, H, lenvs;
-                      solver=Defaults.linearsolver)
-    # Explicitly define optional arguments as these depend on solver,
-    # which needs to come after these arguments.
-    renvs = exci.trivial ? lenvs : environments(exci.right_gs, H; solver=solver)
+# Explicitly define optional arguments as these depend on kwargs,
+# which needs to come after these arguments.
 
-    return environments(exci, H, lenvs, renvs; solver=solver)
+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 = exci.trivial ? lenvs : environments(exci.right_gs, H; kwargs...)
+    return environments(exci, H, lenvs, renvs; kwargs...)
 end
 
-function gen_exci_lw_rw(left_gs::InfiniteMPS, H::Union{InfiniteMPO,InfiniteMPOHamiltonian},
-                        right_gs,
-                        excileg)
-    B = tensormaptype(spacetype(left_gs), 2, 2, storagetype(eltype(left_gs)))
-    BB = tensormaptype(sumspacetype(spacetype(B)), 2, 2, B)
-
-    GBL = PeriodicVector{BB}(undef, length(left_gs))
-    GBR = PeriodicVector{BB}(undef, length(left_gs))
-
-    for site in 1:length(left_gs)
-        GBL[site] = BB(undef,
-                       left_virtualspace(left_gs, site - 1) ⊗ left_virtualspace(H, site)' ←
-                       excileg' ⊗ left_virtualspace(right_gs, site - 1))
-        fill!(GBL[site], zero(scalartype(B)))
-
-        GBR[site] = BB(undef,
-                       right_virtualspace(left_gs, site)' ⊗ right_virtualspace(H, site)' ←
-                       excileg' ⊗ right_virtualspace(right_gs, site)')
-        fill!(GBR[site], zero(scalartype(B)))
+function environments(qp::MultilineQP, operator::MultilineMPO, lenvs, renvs; kwargs...)
+    (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...)
     end
-
-    return GBL, GBR
+    return Multiline(PeriodicVector(envs))
 end
 
-function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
-                      solver=Defaults.linearsolver)
+function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian,
+                      lenvs, renvs;
+                      kwargs...)
     ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1))
-    # ids = collect(Iterators.filter(x -> isidentitylevel(H, x), 2:(H.odim - 1)))
+    solver = environment_alg(exci, H, exci; kwargs...)
 
     AL = exci.left_gs.AL
     AR = exci.right_gs.AR
@@ -72,7 +58,6 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
     lBs = PeriodicVector([allocate_GBL(exci, H, exci, i) for i in 1:length(exci)])
     rBs = PeriodicVector([allocate_GBR(exci, H, exci, i) for i in 1:length(exci)])
 
-    # lBs, rBs = gen_exci_lw_rw(exci.left_gs, H, exci.right_gs, space(exci[1], 3))
     zerovector!(lBs[1])
     for pos in 1:length(exci)
         lBs[pos + 1] = lBs[pos] * TransferMatrix(AR[pos], H[pos], AL[pos]) /
@@ -159,7 +144,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
         rBs[i - 1] += rB_cur
     end
 
-    return QuasiparticleEnvironments(lBs, rBs, lenvs, renvs)
+    return InfiniteQPEnvironments(lBs, rBs, lenvs, renvs)
 end
 
 function environments(exci::FiniteQP,
@@ -188,150 +173,13 @@ function environments(exci::FiniteQP,
                         rightenv(renvs, pos, exci.right_gs)
     end
 
-    return QuasiparticleEnvironments(lBs, rBs, lenvs, renvs)
-end
-
-function environments(exci::Multiline{<:InfiniteQP},
-                      ham::MultilineMPO,
-                      lenvs,
-                      renvs;
-                      solver=Defaults.linearsolver)
-    exci.trivial ||
-        @warn "there is a phase ambiguity in topologically nontrivial statmech excitations"
-
-    left_gs = exci.left_gs
-    right_gs = exci.right_gs
-
-    exci_space = space(exci[1][1], 3)
-
-    (numrows, numcols) = size(left_gs)
-
-    st = site_type(typeof(left_gs))
-    B_type = tensormaptype(spacetype(st), 2, 2, storagetype(st))
-
-    lBs = PeriodicArray{B_type,2}(undef, size(left_gs, 1), size(left_gs, 2))
-    rBs = PeriodicArray{B_type,2}(undef, size(left_gs, 1), size(left_gs, 2))
-
-    for row in 1:numrows
-        c_lenvs = broadcast(col -> leftenv(lenvs, col, left_gs)[row], 1:numcols)
-        c_renvs = broadcast(col -> rightenv(renvs, col, right_gs)[row], 1:numcols)
-
-        hamrow = ham[row, :]
-
-        left_above = left_gs[row]
-        left_below = left_gs[row + 1]
-        right_above = right_gs[row]
-        right_below = right_gs[row + 1]
-
-        left_renorms = fill(zero(scalartype(B_type)), numcols)
-        right_renorms = fill(zero(scalartype(B_type)), numcols)
-
-        for col in 1:numcols
-            lv = leftenv(lenvs, col, left_gs)[row]
-            rv = rightenv(lenvs, col, left_gs)[row]
-            left_renorms[col] = @plansor lv[1 2; 3] *
-                                         left_above.AC[col][3 4; 5] *
-                                         hamrow[col][2 6; 4 7] *
-                                         rv[5 7; 8] *
-                                         conj(left_below.AC[col][1 6; 8])
-
-            lv = leftenv(renvs, col, right_gs)[row]
-            rv = rightenv(renvs, col, right_gs)[row]
-            right_renorms[col] = @plansor lv[1 2; 3] *
-                                          right_above.AC[col][3 4; 5] *
-                                          hamrow[col][2 6; 4 7] *
-                                          rv[5 7; 8] *
-                                          conj(right_below.AC[col][1 6; 8])
-        end
-
-        left_renorms = left_renorms .^ -1
-        right_renorms = right_renorms .^ -1
-
-        lB_cur = fill_data!(similar(left_below.AL[1],
-                                    left_virtualspace(left_below, 0) *
-                                    _firstspace(hamrow[1])',
-                                    exci_space' * right_virtualspace(right_above, 0)),
-                            zero)
-        rB_cur = fill_data!(similar(left_below.AL[1],
-                                    left_virtualspace(left_below, 0) *
-                                    _firstspace(hamrow[1]),
-                                    exci_space' * right_virtualspace(right_above, 0)),
-                            zero)
-        for col in 1:numcols
-            lB_cur = lB_cur *
-                     TransferMatrix(right_above.AR[col], hamrow[col], left_below.AL[col])
-            lB_cur += c_lenvs[col] *
-                      TransferMatrix(exci[row][col], hamrow[col], left_below.AL[col])
-            lB_cur *= left_renorms[col] * exp(-1im * exci.momentum)
-            lBs[row, col] = lB_cur
-
-            col = numcols - col + 1
-
-            rB_cur = TransferMatrix(left_above.AL[col], hamrow[col], right_below.AR[col]) *
-                     rB_cur
-            rB_cur += TransferMatrix(exci[row][col], hamrow[col], right_below.AR[col]) *
-                      c_renvs[col]
-            rB_cur *= exp(1im * exci.momentum) * right_renorms[col]
-            rBs[row, col] = rB_cur
-        end
-
-        tm_RL = TransferMatrix(right_above.AR, hamrow, left_below.AL)
-        tm_LR = TransferMatrix(left_above.AL, hamrow, right_below.AR)
-
-        if exci.trivial
-            @plansor rvec[-1 -2; -3] := rightenv(lenvs, 0, left_gs)[row][-1 -2; 1] *
-                                        conj(left_below.C[0][-3; 1])
-            @plansor lvec[-1 -2; -3] := leftenv(lenvs, 1, left_gs)[row][-1 -2; 1] *
-                                        left_above.C[0][1; -3]
-
-            tm_RL = regularize(tm_RL, lvec, rvec)
-
-            @plansor rvec[-1 -2; -3] := rightenv(renvs, 0, right_gs)[row][1 -2; -3] *
-                                        right_above.C[0][-1; 1]
-            @plansor lvec[-1 -2; -3] := conj(right_below.C[0][-3; 1]) *
-                                        leftenv(renvs, 1, right_gs)[row][-1 -2; 1]
-
-            tm_LR = regularize(tm_LR, lvec, rvec)
-        end
-
-        lBs[row, end], convhist = linsolve(flip(tm_RL), lB_cur, lB_cur, solver, 1,
-                                           -exp(-1im * numcols * exci.momentum) *
-                                           prod(left_renorms))
-        convhist.converged == 0 &&
-            @warn "GBL[$row] failed to converge: normres = $(convhist.normres)"
-
-        rBs[row, 1], convhist = linsolve(tm_LR, rB_cur, rB_cur, GMRES(), 1,
-                                         -exp(1im * numcols * exci.momentum) *
-                                         prod(right_renorms))
-        convhist.converged == 0 &&
-            @warn "GBR[$row] failed to converge: normres = $(convhist.normres)"
-
-        left_cur = lBs[row, end]
-        right_cur = rBs[row, 1]
-        for col in 1:(numcols - 1)
-            left_cur = left_renorms[col] * left_cur *
-                       TransferMatrix(right_above.AR[col], hamrow[col],
-                                      left_below.AL[col]) * exp(-1im * exci.momentum)
-            lBs[row, col] += left_cur
-
-            col = numcols - col + 1
-            right_cur = TransferMatrix(left_above.AL[col], hamrow[col],
-                                       right_below.AR[col]) * right_cur *
-                        exp(1im * exci.momentum) * right_renorms[col]
-            rBs[row, col] += right_cur
-        end
-    end
-
-    return QPEnv(lBs, rBs, lenvs, renvs)
+    return InfiniteQPEnvironments(lBs, rBs, lenvs, renvs)
 end
 
-function environments(exci::InfiniteQP,
-                      O::InfiniteMPO,
-                      lenvs,
-                      renvs;
-                      solver=Defaults.linearsolver)
+function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...)
     exci.trivial ||
         @warn "there is a phase ambiguity in topologically nontrivial statmech excitations"
+    solver = environment_alg(exci, O, exci; kwargs...)
 
     left_gs = exci.left_gs
     right_gs = exci.right_gs
@@ -415,5 +263,5 @@ function environments(exci::InfiniteQP,
         GBR[col] += right_cur
     end
 
-    return QuasiparticleEnvironments(GBL, GBR, lenvs, renvs)
+    return InfiniteQPEnvironments(GBL, GBR, lenvs, renvs)
 end
diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl
index 416d20774..77b55fbaa 100644
--- a/src/operators/multipliedoperator.jl
+++ b/src/operators/multipliedoperator.jl
@@ -52,3 +52,6 @@ Base.:*(op::UntimedOperator, g::Function) = MultipliedOperator(op.op, t -> g(t)
 function environments(st, x::MultipliedOperator, args...; kwargs...)
     return environments(st, x.op, args...; kwargs...)
 end
+function recalculate!(envs, above, x::MultipliedOperator, below=above; kwargs...)
+    return recalculate!(envs, above, x.op, below; kwargs...)
+end
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index e96ea640b..2072c28be 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -280,6 +280,8 @@ function Base.similar(ψ::FiniteMPS{A,B}) where {A,B}
     return FiniteMPS{A,B}(similar(ψ.ALs), similar(ψ.ARs), similar(ψ.ACs), similar(ψ.Cs))
 end
 
+Base.eachindex(ψ::FiniteMPS) = eachindex(ψ.AL)
+Base.eachindex(l::IndexStyle, ψ::FiniteMPS) = eachindex(l, ψ.AL)
 Base.checkbounds(::Type{Bool}, ψ::FiniteMPS, i::Integer) = 1 <= i <= length(ψ)
 
 Base.@propagate_inbounds function Base.getindex(ψ::FiniteMPS, i::Int)
diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl
index eddb593e4..9c7561836 100644
--- a/src/states/infinitemps.jl
+++ b/src/states/infinitemps.jl
@@ -216,6 +216,14 @@ Base.size(ψ::InfiniteMPS, args...) = size(ψ.AL, args...)
 Base.length(ψ::InfiniteMPS) = length(ψ.AL)
 Base.eltype(ψ::InfiniteMPS) = eltype(ψ.AL)
 Base.copy(ψ::InfiniteMPS) = InfiniteMPS(copy(ψ.AL), copy(ψ.AR), copy(ψ.C), copy(ψ.AC))
+function Base.copy!(ψ::InfiniteMPS, ϕ::InfiniteMPS)
+    copy!.(ψ.AL, ϕ.AL)
+    copy!.(ψ.AR, ϕ.AR)
+    copy!.(ψ.AC, ϕ.AC)
+    copy!.(ψ.C, ϕ.C)
+    return ψ
+end
+
 function Base.repeat(ψ::InfiniteMPS, i::Int)
     return InfiniteMPS(repeat(ψ.AL, i), repeat(ψ.AR, i), repeat(ψ.C, i), repeat(ψ.AC, i))
 end
@@ -227,6 +235,9 @@ function Base.circshift(ψ::InfiniteMPS, n)
                        circshift(ψ.AC, n))
 end
 
+Base.eachindex(ψ::InfiniteMPS) = eachindex(ψ.AL)
+Base.eachindex(l::IndexStyle, ψ::InfiniteMPS) = eachindex(l, ψ.AL)
+
 Base.checkbounds(::Type{Bool}, ψ::InfiniteMPS, i::Integer) = true
 
 site_type(::Type{<:InfiniteMPS{A}}) where {A} = A
diff --git a/src/states/multilinemps.jl b/src/states/multilinemps.jl
index b8e1eeafe..d5ff0ef2e 100644
--- a/src/states/multilinemps.jl
+++ b/src/states/multilinemps.jl
@@ -83,5 +83,8 @@ TensorKit.normalize!(a::MultilineMPS) = (normalize!.(parent(a)); return a)
 Base.convert(::Type{MultilineMPS}, st::InfiniteMPS) = Multiline([st])
 Base.convert(::Type{InfiniteMPS}, st::MultilineMPS) = only(st)
 Base.eltype(t::MultilineMPS) = eltype(t[1])
+Base.copy!(ψ::MultilineMPS, ϕ::MultilineMPS) = (copy!.(parent(ψ), parent(ϕ));
+                                                ψ)
+
 left_virtualspace(t::MultilineMPS, i::Int, j::Int) = left_virtualspace(t[i], j)
 right_virtualspace(t::MultilineMPS, i::Int, j::Int) = right_virtualspace(t[i], j)
diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl
index 01a203bc3..10f36cbe0 100644
--- a/src/states/quasiparticle_state.jl
+++ b/src/states/quasiparticle_state.jl
@@ -174,6 +174,7 @@ end
 const QP{S,T1,T2} = Union{LeftGaugedQP{S,T1,T2},RightGaugedQP{S,T1,T2}}
 const FiniteQP{S<:FiniteMPS,T1,T2} = QP{S,T1,T2}
 const InfiniteQP{S<:InfiniteMPS,T1,T2} = QP{S,T1,T2}
+const MultilineQP{Q<:QP} = Multiline{Q}
 
 TensorKit.spacetype(::Union{QP{S},Type{<:QP{S}}}) where {S} = spacetype(S)
 TensorKit.sectortype(::Union{QP{S},Type{<:QP{S}}}) where {S} = sectortype(S)
@@ -350,7 +351,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S<:FiniteMPS}
     return FiniteMPS(Ls + Rs + Bs; normalize=false)
 end
 
-function Base.getproperty(exci::Multiline{<:InfiniteQP}, s::Symbol)
+function Base.getproperty(exci::MultilineQP, s::Symbol)
     if s == :momentum
         return first(exci.data).momentum
     elseif s == :left_gs
diff --git a/src/utility/multiline.jl b/src/utility/multiline.jl
index 281bfe666..2b645176f 100644
--- a/src/utility/multiline.jl
+++ b/src/utility/multiline.jl
@@ -28,6 +28,7 @@ function Base.axes(m::Multiline, i::Int)
     return i == 1 ? axes(parent(m), 1) :
            i == 2 ? axes(parent(m)[1], 1) : throw(ArgumentError("Invalid index $i"))
 end
+Base.eachindex(m::Multiline) = CartesianIndices(size(m))
 
 Base.getindex(m::Multiline, i::Int) = getindex(parent(m), i)
 Base.setindex!(m::Multiline, v, i::Int) = (setindex!(parent(m), v, i); m)
diff --git a/test/algorithms.jl b/test/algorithms.jl
index 6d85fa129..97a924475 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -690,7 +690,7 @@ end
 end
 
 @testset "approximate" verbose = true begin
-    verbosity = 0
+    verbosity = verbosity_conv
     @testset "mpo * infinite ≈ infinite" begin
         ψ = InfiniteMPS([ℙ^2, ℙ^2], [ℙ^10, ℙ^10])
         H = force_planar(repeat(transverse_field_ising(; g=4), 2))