diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 29c324a7c..142046e80 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -27,7 +27,7 @@ export entanglementplot, transferplot
 
 # hamiltonian things
 export Cache
-export SparseMPO, MPOHamiltonian, DenseMPO, MPOMultiline
+export SparseMPO, MPOHamiltonian, DenseMPO, MPOMultiline, LazySum
 export ∂C, ∂AC, ∂AC2, environments, expectation_value, effective_excitation_hamiltonian
 export leftenv, rightenv
 
@@ -79,6 +79,7 @@ include("operators/sparsempo/sparsempo.jl")
 include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
 include("operators/mpomultiline.jl")
 include("operators/projection.jl")
+include("operators/lazysum.jl")
 
 include("transfermatrix/transfermatrix.jl")
 include("transfermatrix/transfer.jl")
@@ -90,6 +91,7 @@ include("environments/abstractinfenv.jl")
 include("environments/permpoinfenv.jl")
 include("environments/mpohaminfenv.jl")
 include("environments/qpenv.jl")
+include("environments/multipleenv.jl")
 include("environments/idmrgenv.jl")
 include("environments/lazylincocache.jl")
 
diff --git a/src/algorithms/derivatives.jl b/src/algorithms/derivatives.jl
index 9470efc6e..c7c1ece9a 100644
--- a/src/algorithms/derivatives.jl
+++ b/src/algorithms/derivatives.jl
@@ -21,6 +21,8 @@ struct MPO_∂∂AC2{O,L,R}
     rightenv::R
 end
 
+const DerivativeOperator = Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}
+
 Base.:*(h::Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}, v) = h(v);
 
 (h::MPO_∂∂C)(x) = ∂C(x, h.leftenv, h.rightenv);
@@ -326,3 +328,24 @@ function ∂∂AC2(pos::Int, state, opp::ProjectionOperator, env)
         rightenv(env, pos + 1, state),
     )
 end
+
+(x::LazySum{<:MPSKit.DerivativeOperator})(y) = sum(O -> O(y), x)
+Base.:*(h::LazySum{<:MPSKit.DerivativeOperator}, v) = h(v)
+
+function ∂∂C(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{MPO_∂∂C}(
+        map((op, openv) -> ∂∂C(pos, mps, op, openv), opp.ops, cache.envs)
+    )
+end
+
+function ∂∂AC(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{MPO_∂∂AC}(
+        map((op, openv) -> ∂∂AC(pos, mps, op, openv), opp.ops, cache.envs)
+    )
+end
+
+function ∂∂AC2(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{MPO_∂∂AC2}(
+        map((op, openv) -> ∂∂AC2(pos, mps, op, openv), opp.ops, cache.envs)
+    )
+end
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index 2cea5c151..2d3c6ad3b 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -1,13 +1,30 @@
-#works for general tensors
+"""
+    expectation_value(ψ, O, [location], [environments])
+
+Compute the expectation value of an operator `O` on a state `ψ`. If `location` is given, the
+operator is applied at that location. If `environments` is given, the expectation value
+might be computed more efficiently by re-using previously calculated environments.
+
+# Arguments
+* `ψ::AbstractMPS` : the state on which to compute the expectation value
+* `O` : the operator to compute the expectation value of. This can either be an `AbstractMPO`, a single `AbstractTensorMap` or an array of `AbstractTensorMap`s.
+* `location::Union{Int,AbstractRange{Int}}` : the location at which to apply the operator. Only applicable for operators that act on a subset of all sites.
+* `environments::Cache` : the environments to use for the calculation. If not given, they will be calculated.
+"""
+function expectation_value end
+
+# Local operators
+# ---------------
 function expectation_value(
-    state::Union{InfiniteMPS,WindowMPS,FiniteMPS}, opp::AbstractTensorMap
-)
-    return expectation_value(state, fill(opp, length(state)))
+    ψ::Union{InfiniteMPS,WindowMPS,FiniteMPS}, O::AbstractTensorMap{S,1,1}
+) where {S}
+    return expectation_value(ψ, fill(O, length(ψ)))
 end
 function expectation_value(
-    state::Union{InfiniteMPS,WindowMPS,FiniteMPS}, opps::AbstractArray{<:AbstractTensorMap}
-)
-    map(zip(state.AC, opps)) do (ac, opp)
+    ψ::Union{InfiniteMPS,WindowMPS,FiniteMPS},
+    opps::AbstractArray{<:AbstractTensorMap{S,1,1}},
+) where {S}
+    return map(zip(ψ.AC, opps)) do (ac, opp)
         tr(
             ac' * transpose(
                 opp * transpose(
@@ -22,86 +39,82 @@ function expectation_value(
     end
 end
 
-"""
-calculates the expectation value of op, where op is a plain tensormap where the first index works on site at
-"""
+# Multi-site operators
+# --------------------
+# TODO: replace Vector{MPOTensor} with FiniteMPO
 function expectation_value(
-    state::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, op::AbstractTensorMap, at::Int
-) where {T<:MPSTensor}
-    return expectation_value(state, decompose_localmpo(add_util_leg(op)), at)
+    ψ::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, O::AbstractTensorMap{S,N,N}, at::Int
+) where {S,N,T<:MPSTensor{S}}
+    return expectation_value(ψ, decompose_localmpo(add_util_leg(O)), at)
 end
-
-"""
-calculates the expectation value of op = op1*op2*op3*... (ie an N site operator) starting at site at
-"""
 function expectation_value(
-    state::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}},
-    op::AbstractArray{<:AbstractTensorMap},
+    ψ::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}},
+    O::AbstractArray{<:MPOTensor{S}},
     at::Int,
-) where {T<:MPSTensor}
-    firstspace = _firstspace(first(op))
-    (firstspace == oneunit(firstspace) && _lastspace(last(op)) == firstspace') || throw(
+) where {S,T<:MPSTensor{S}}
+    firstspace = _firstspace(first(O))
+    (firstspace == oneunit(firstspace) && _lastspace(last(O)) == firstspace') || throw(
         ArgumentError(
             "localmpo should start and end in a trivial leg, not with $(firstspace)"
         ),
     )
 
-    ut = fill_data!(similar(op[1], firstspace), one)
+    ut = fill_data!(similar(O[1], firstspace), one)
     @plansor v[-1 -2; -3] :=
         isomorphism(
-            storagetype(T),
-            left_virtualspace(state, at - 1),
-            left_virtualspace(state, at - 1),
+            storagetype(T), left_virtualspace(ψ, at - 1), left_virtualspace(ψ, at - 1)
         )[
             -1
             -3
         ] * conj(ut[-2])
     tmp =
-        v * TransferMatrix(
-            state.AL[at:(at + length(op) - 1)], op, state.AL[at:(at + length(op) - 1)]
-        )
+        v * TransferMatrix(ψ.AL[at:(at + length(O) - 1)], O, ψ.AL[at:(at + length(O) - 1)])
     return @plansor tmp[1 2; 3] *
         ut[2] *
-        state.CR[at + length(op) - 1][3; 4] *
-        conj(state.CR[at + length(op) - 1][1; 4])
+        ψ.CR[at + length(O) - 1][3; 4] *
+        conj(ψ.CR[at + length(O) - 1][1; 4])
+end
+
+# MPOHamiltonian
+# --------------
+# TODO: add section in documentation to explain convention
+expectation_value(ψ, envs::Cache) = expectation_value(ψ, envs.opp, envs)
+function expectation_value(ψ, ham::MPOHamiltonian)
+    return expectation_value(ψ, ham, environments(ψ, ham))
 end
 
 """
-    calculates the expectation value for the given operator/hamiltonian
+    expectation_value(ψ::WindowMPS, ham::MPOHAmiltonian, envs) -> vals, tot
+
+TODO
 """
-expectation_value(state, envs::Cache) = expectation_value(state, envs.opp, envs);
-function expectation_value(state, ham::MPOHamiltonian)
-    return expectation_value(state, ham, environments(state, ham))
-end
-function expectation_value(state::WindowMPS, ham::MPOHamiltonian, envs::FinEnv)
-    vals = expectation_value_fimpl(state, ham, envs)
+function expectation_value(ψ::WindowMPS, ham::MPOHamiltonian, envs::FinEnv)
+    vals = expectation_value_fimpl(ψ, ham, envs)
 
     tot = 0.0 + 0im
     for i in 1:(ham.odim), j in 1:(ham.odim)
-        tot += @plansor leftenv(envs, length(state), state)[i][1 2; 3] *
-            state.AC[end][3 4; 5] *
-            rightenv(envs, length(state), state)[j][5 6; 7] *
-            ham[length(state)][i, j][2 8; 4 6] *
-            conj(state.AC[end][1 8; 7])
+        tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] *
+            ψ.AC[end][3 4; 5] *
+            rightenv(envs, length(ψ), ψ)[j][5 6; 7] *
+            ham[length(ψ)][i, j][2 8; 4 6] *
+            conj(ψ.AC[end][1 8; 7])
     end
 
-    return vals, tot / (norm(state.AC[end])^2)
+    return vals, tot / (norm(ψ.AC[end])^2)
 end
 
-function expectation_value(state::FiniteMPS, ham::MPOHamiltonian, envs::FinEnv)
-    return expectation_value_fimpl(state, ham, envs)
+function expectation_value(ψ::FiniteMPS, ham::MPOHamiltonian, envs::FinEnv)
+    return expectation_value_fimpl(ψ, ham, envs)
 end
-function expectation_value_fimpl(
-    state::AbstractFiniteMPS, ham::MPOHamiltonian, envs::FinEnv
-)
-    ens = zeros(scalartype(state), length(state))
-    for i in 1:length(state), (j, k) in keys(ham[i])
+function expectation_value_fimpl(ψ::AbstractFiniteMPS, ham::MPOHamiltonian, envs::FinEnv)
+    ens = zeros(scalartype(ψ), length(ψ))
+    for i in 1:length(ψ), (j, k) in keys(ham[i])
         !((j == 1 && k != 1) || (k == ham.odim && j != ham.odim)) && continue
 
-        cur = @plansor leftenv(envs, i, state)[j][1 2; 3] *
-            state.AC[i][3 7; 5] *
-            rightenv(envs, i, state)[k][5 8; 6] *
-            conj(state.AC[i][1 4; 6]) *
+        cur = @plansor leftenv(envs, i, ψ)[j][1 2; 3] *
+            ψ.AC[i][3 7; 5] *
+            rightenv(envs, i, ψ)[k][5 8; 6] *
+            conj(ψ.AC[i][1 4; 6]) *
             ham[i][j, k][2 4; 7 8]
         if !(j == 1 && k == ham.odim)
             cur /= 2
@@ -110,7 +123,7 @@ function expectation_value_fimpl(
         ens[i] += cur
     end
 
-    n = norm(state.AC[end])^2
+    n = norm(ψ.AC[end])^2
     return ens ./ n
 end
 
@@ -164,28 +177,46 @@ function expectation_value(
     return tot
 end
 
-function expectation_value(st::InfiniteMPS, mpo::DenseMPO)
-    return expectation_value(convert(MPSMultiline, st), convert(MPOMultiline, mpo))
-end;
-function expectation_value(st::MPSMultiline, mpo::MPOMultiline)
-    return expectation_value(st, environments(st, mpo))
-end;
-function expectation_value(st::InfiniteMPS, ca::PerMPOInfEnv)
-    return expectation_value(convert(MPSMultiline, st), ca)
-end;
-function expectation_value(st::MPSMultiline, opp::MPOMultiline, ca::PerMPOInfEnv)
-    retval = PeriodicArray{scalartype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2))
-    for (i, j) in product(1:size(st, 1), 1:size(st, 2))
-        retval[i, j] = @plansor leftenv(ca, i, j, st)[1 2; 3] *
-            opp[i, j][2 4; 6 5] *
-            st.AC[i, j][3 6; 7] *
-            rightenv(ca, i, j, st)[7 5; 8] *
-            conj(st.AC[i + 1, j][1 4; 8])
+# Transfer matrices
+# -----------------
+function expectation_value(ψ::InfiniteMPS, mpo::DenseMPO)
+    return expectation_value(convert(MPSMultiline, ψ), convert(MPOMultiline, mpo))
+end
+function expectation_value(ψ::MPSMultiline, mpo::MPOMultiline)
+    return expectation_value(ψ, environments(ψ, mpo))
+end
+function expectation_value(ψ::InfiniteMPS, ca::PerMPOInfEnv)
+    return expectation_value(convert(MPSMultiline, ψ), ca)
+end
+function expectation_value(ψ::MPSMultiline, O::MPOMultiline, ca::PerMPOInfEnv)
+    retval = PeriodicMatrix{scalartype(ψ)}(undef, size(ψ, 1), size(ψ, 2))
+    for (i, j) in product(1:size(ψ, 1), 1:size(ψ, 2))
+        retval[i, j] = @plansor leftenv(ca, i, j, ψ)[1 2; 3] *
+            O[i, j][2 4; 6 5] *
+            ψ.AC[i, j][3 6; 7] *
+            rightenv(ca, i, j, ψ)[7 5; 8] *
+            conj(ψ.AC[i + 1, j][1 4; 8])
     end
     return retval
 end
 
-expectation_value(state::FiniteQP, opp) = expectation_value(convert(FiniteMPS, state), opp)
-function expectation_value(state::FiniteQP, opp::MPOHamiltonian)
-    return expectation_value(convert(FiniteMPS, state), opp)
+expectation_value(ψ::FiniteQP, O) = expectation_value(convert(FiniteMPS, ψ), O)
+function expectation_value(ψ::FiniteQP, O::MPOHamiltonian)
+    return expectation_value(convert(FiniteMPS, ψ), O)
+end
+
+# define expectation_value for LazySum
+function expectation_value(ψ, ops::LazySum, at::Int)
+    return sum(op -> expectation_value(ψ, op, at), ops)
+end
+function expectation_value(Ψ, ops::LazySum, envs::MultipleEnvironments=environments(Ψ, ops))
+    return sum(((op, env),) -> expectation_value(Ψ, op, env), zip(ops.ops, envs))
+end
+
+# for now we also have LinearCombination
+function expectation_value(Ψ, H::LinearCombination, envs::LazyLincoCache=environments(Ψ, H))
+    return return sum(
+        ((c, op, env),) -> c * expectation_value(Ψ, op, env),
+        zip(H.coeffs, H.opps, envs.envs),
+    )
 end
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index 214d45b58..c22d25286 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -35,8 +35,7 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG1, oenvs=environments
             Ψ.AC[pos] = vecs[1]
             Ψ.AL[pos], Ψ.CR[pos] = leftorth(vecs[1])
 
-            tm = TransferMatrix(Ψ.AL[pos], ham[pos], Ψ.AL[pos])
-            setleftenv!(envs, pos + 1, leftenv(envs, pos) * tm)
+            update_leftenv!(envs, Ψ, ham, pos + 1)
         end
 
         for pos in length(Ψ):-1:1
@@ -47,8 +46,7 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG1, oenvs=environments
             Ψ.CR[pos - 1], temp = rightorth(_transpose_tail(vecs[1]))
             Ψ.AR[pos] = _transpose_front(temp)
 
-            tm = TransferMatrix(Ψ.AR[pos], ham[pos], Ψ.AR[pos])
-            setrightenv!(envs, pos - 1, tm * rightenv(envs, pos))
+            update_rightenv!(envs, Ψ, ham, pos - 1)
         end
 
         delta = norm(curc - Ψ.CR[0])
@@ -110,10 +108,8 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG2, oenvs=environments
             st.AR[pos + 1] = _transpose_front(ar)
             st.AC[pos + 1] = _transpose_front(c * ar)
 
-            tm_L = TransferMatrix(st.AL[pos], ham[pos], st.AL[pos])
-            setleftenv!(envs, pos + 1, leftenv(envs, pos) * tm_L)
-            tm_R = TransferMatrix(st.AR[pos + 1], ham[pos + 1], st.AR[pos + 1])
-            setrightenv!(envs, pos, tm_R * rightenv(envs, pos + 1))
+            update_leftenv!(envs, st, ham, pos + 1)
+            update_rightenv!(envs, st, ham, pos)
         end
 
         # update the edge
@@ -135,10 +131,8 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG2, oenvs=environments
         curc = complex(c)
 
         # update environments
-        tm_L = TransferMatrix(st.AL[0], ham[0], st.AL[0])
-        setleftenv!(envs, 1, leftenv(envs, 0) * tm_L)
-        tm_R = TransferMatrix(st.AR[1], ham[1], st.AR[1])
-        setrightenv!(envs, 0, tm_R * rightenv(envs, 1))
+        update_leftenv!(envs, st, ham, 1)
+        update_rightenv!(envs, st, ham, 0)
 
         # sweep from right to left
         for pos in (length(st) - 1):-1:1
@@ -155,10 +149,8 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG2, oenvs=environments
             st.AR[pos + 1] = _transpose_front(ar)
             st.AC[pos + 1] = _transpose_front(c * ar)
 
-            tm_L = TransferMatrix(st.AL[pos], ham[pos], st.AL[pos])
-            setleftenv!(envs, pos + 1, leftenv(envs, pos) * tm_L)
-            tm_R = TransferMatrix(st.AR[pos + 1], ham[pos + 1], st.AR[pos + 1])
-            setrightenv!(envs, pos, tm_R * rightenv(envs, pos + 1))
+            update_leftenv!(envs, st, ham, pos + 1)
+            update_rightenv!(envs, st, ham, pos)
         end
 
         # update the edge
@@ -178,10 +170,8 @@ function find_groundstate(ost::InfiniteMPS, ham, alg::IDMRG2, oenvs=environments
         st.AR[1] = _transpose_front(ar)
         st.AC[1] = _transpose_front(c * ar)
 
-        tm_L = TransferMatrix(st.AL[0], ham[0], st.AL[0])
-        setleftenv!(envs, 1, leftenv(envs, 0) * tm_L)
-        tm_R = TransferMatrix(st.AR[1], ham[1], st.AR[1])
-        setrightenv!(envs, 0, tm_R * rightenv(envs, 1))
+        update_leftenv!(envs, st, ham, 1)
+        update_rightenv!(envs, st, ham, 0)
 
         # update error
         smallest = infimum(_firstspace(curc), _firstspace(c))
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 7762b03ca..54030736e 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -190,6 +190,13 @@ function variance(state::InfiniteQP, ham::MPOHamiltonian, envs=environments(stat
     )
 end
 
+function variance(Ψ, ham::LazySum, envs=environments(Ψ, sum(ham)))
+    envs isa MultipleEnvironments && throw(
+        ArgumentError("The environment cannot be Lazy i.e. use environments of sum(H)")
+    )
+    return variance(Ψ, sum(ham), envs)
+end
+
 """
 You can impose periodic boundary conditions on an mpo-hamiltonian (for a given size)
 That creates a new mpo-hamiltonian with larger bond dimension
diff --git a/src/environments/idmrgenv.jl b/src/environments/idmrgenv.jl
index ae452ed93..a62354037 100644
--- a/src/environments/idmrgenv.jl
+++ b/src/environments/idmrgenv.jl
@@ -9,14 +9,14 @@ struct IDMRGEnv{H,V}
     rw::PeriodicArray{V,2}
 end
 
-function IDMRGEnv(state::Union{MPSMultiline,InfiniteMPS}, env)
-    state === env.dependency || recalculate!(env, state)
+function IDMRGEnv(Ψ::Union{MPSMultiline,InfiniteMPS}, env)
+    Ψ === env.dependency || recalculate!(env, Ψ)
     return IDMRGEnv(env.opp, deepcopy(env.lw), deepcopy(env.rw))
 end
 
 leftenv(envs::IDMRGEnv, pos::Int) = envs.lw[:, pos];
 leftenv(envs::IDMRGEnv, row::Int, col::Int) = envs.lw[row, col];
-leftenv(envs::IDMRGEnv, pos::Int, state::InfiniteMPS) = envs.lw[:, pos];
+leftenv(envs::IDMRGEnv, pos::Int, Ψ::InfiniteMPS) = envs.lw[:, pos];
 function setleftenv!(envs::IDMRGEnv, pos, lw)
     return envs.lw[:, pos] = lw[:]
 end
@@ -26,10 +26,49 @@ end
 
 rightenv(envs::IDMRGEnv, row::Int, col::Int) = envs.rw[row, col];
 rightenv(envs::IDMRGEnv, pos::Int) = envs.rw[:, pos];
-rightenv(envs::IDMRGEnv, pos::Int, state::InfiniteMPS) = envs.rw[:, pos];
+rightenv(envs::IDMRGEnv, pos::Int, Ψ::InfiniteMPS) = envs.rw[:, pos];
 function setrightenv!(envs::IDMRGEnv, pos, rw)
     return envs.rw[:, pos] = rw[:]
 end
 function setrightenv!(envs::IDMRGEnv, row, col, val)
     return envs.rw[row, col] = val
 end
+
+# For MultipleEnvironments
+
+function IDMRGEnv(Ψ::Union{MPSMultiline,InfiniteMPS}, env::MultipleEnvironments)
+    tmp = map(env.envs) do subenv
+        Ψ === subenv.dependency || recalculate!(subenv, Ψ)
+        (subenv.opp, IDMRGEnv(subenv.opp, deepcopy(subenv.lw), deepcopy(subenv.rw)))
+    end
+    hams, envs = collect.(zip(tmp...))
+    return MultipleEnvironments(LazySum(hams), envs)
+end
+
+function update_rightenv!(
+    envs::MultipleEnvironments{<:LazySum,<:IDMRGEnv}, st, ham, pos::Int
+)
+    for (subham, subenv) in zip(ham, envs.envs)
+        tm = TransferMatrix(st.AR[pos + 1], subham[pos + 1], st.AR[pos + 1])
+        setrightenv!(subenv, pos, tm * rightenv(subenv, pos + 1))
+    end
+end
+
+function update_leftenv!(
+    envs::MultipleEnvironments{<:LazySum,<:IDMRGEnv}, st, ham, pos::Int
+)
+    for (subham, subenv) in zip(ham, envs.envs)
+        tm = TransferMatrix(st.AL[pos - 1], subham[pos - 1], st.AL[pos - 1])
+        setleftenv!(subenv, pos, leftenv(subenv, pos - 1) * tm)
+    end
+end
+
+function update_rightenv!(envs::IDMRGEnv, st, ham, pos::Int)
+    tm = TransferMatrix(st.AR[pos + 1], ham[pos + 1], st.AR[pos + 1])
+    return setrightenv!(envs, pos, tm * rightenv(envs, pos + 1))
+end
+
+function update_leftenv!(envs::IDMRGEnv, st, ham, pos::Int)
+    tm = TransferMatrix(st.AL[pos - 1], ham[pos - 1], st.AL[pos - 1])
+    return setleftenv!(envs, pos, leftenv(envs, pos - 1) * tm)
+end
diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl
new file mode 100644
index 000000000..ebb4e9309
--- /dev/null
+++ b/src/environments/multipleenv.jl
@@ -0,0 +1,69 @@
+struct MultipleEnvironments{O,C} <: Cache
+    opp::O
+    envs::Vector{C}
+end
+
+Base.size(x::MultipleEnvironments) = size(x.envs)
+Base.getindex(x::MultipleEnvironments, i) = x.envs[i]
+Base.length(x::MultipleEnvironments) = prod(size(x))
+
+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, ham::LazySum)
+    return MultipleEnvironments(ham, map(op -> environments(st, op), ham.ops))
+end
+
+function environments(
+    st::Union{InfiniteMPS,MPSMultiline}, ham::LazySum; solver=Defaults.linearsolver
+)
+    if !(solver isa Vector)
+        solver = repeat([solver], length(ham))
+    end
+    return MultipleEnvironments(
+        ham, map((op, solv) -> environments(st, op; solver=solv), ham.ops, solver)
+    )
+end
+
+#broadcast vs map?
+# function environments(state, ham::LinearCombination)
+#     return MultipleEnvironments(ham, broadcast(o -> environments(state, o), ham.opps))
+# end;
+
+function environments(
+    st::WindowMPS,
+    ham::LazySum;
+    lenvs=environments(st.left_gs, ham),
+    renvs=environments(st.right_gs, ham),
+)
+    return MultipleEnvironments(
+        ham,
+        map(
+            (op, sublenv, subrenv) -> environments(st, op; lenvs=sublenv, renvs=subrenv),
+            ham.ops,
+            lenvs,
+            renvs,
+        ),
+    )
+end
+
+# we need to define how to recalculate
+"""
+    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
+
+#maybe this can be used to provide compatibility with existing code?
+function Base.getproperty(envs::MultipleEnvironments, prop::Symbol)
+    if prop === :solver
+        return map(env -> env.solver, envs)
+    else
+        return getfield(envs, prop)
+    end
+end
diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl
new file mode 100644
index 000000000..36f3a6a1a
--- /dev/null
+++ b/src/operators/lazysum.jl
@@ -0,0 +1,72 @@
+"""
+    LazySum{O} <: AbstractVector{O}
+
+Type that represents a lazy sum i.e explicit summation is only done when needed. 
+This type is basically an `AbstractVector` with some extra functionality to calculate things efficiently.
+
+## Fields
+- ops -- Vector of summable objects
+
+---
+
+## Constructors
+    LazySum(x::Vector)
+
+"""
+struct LazySum{O} <: AbstractVector{O}
+    ops::Vector{O}
+end
+
+Base.size(x::LazySum) = size(x.ops)
+Base.getindex(x::LazySum, i) = x.ops[i]
+#iteration and summation gets automatically implementend thanks to subtyping
+
+Base.length(x::LazySum) = prod(size(x))
+Base.similar(x::LazySum, ::Type{S}, dims::Dims) where {S} = LazySum(similar(x.ops, S, dims))
+
+# Holy traits
+struct TimeDependent end
+struct NotTimeDependent end
+
+TimeDependence(x) = NotTimeDependent()
+TimeDependence(x::LazySum) = istimed(x) ? TimeDependent() : NotTimeDependent()
+
+istimed(::TimeDependent) = true
+istimed(::NotTimeDependent) = false
+istimed(x) = istimed(TimeDependence(x))
+istimed(x::LazySum) = any(istimed, x)
+
+# constructors
+LazySum(x) = LazySum([x])
+LazySum(f::Function, x) = LazySum(map(y -> f(y), x))
+
+# Internal use only, works always
+_eval_at(x, args...) = x # -> this is what you should define for your custom structs inside a LazySum
+
+# wrapper around _eval_at
+eval_at(x, args...) = eval_at(TimeDependence(x), x, args...)
+eval_at(::TimeDependent, x::LazySum, t::Number) = LazySum(O -> _eval_at(O, t), x)
+function eval_at(::TimeDependent, x::LazySum)
+    throw(
+        ArgumentError(
+            "attempting to evaluate time-dependent LazySum without specifiying a time"
+        ),
+    )
+end
+eval_at(::NotTimeDependent, x::LazySum) = sum(O -> _eval_at(O), x)
+function eval_at(::NotTimeDependent, x::LazySum, t::Number)
+    throw(ArgumentError("attempting to evaluate time-independent LazySum at time"))
+end
+
+# For users
+(x::LazySum)(t::Number) = eval_at(x, t) # using (t) should return NotTimeDependent LazySum
+
+# we define the addition for LazySum and we do the rest with this
+function Base.:+(SumOfOps1::LazySum, SumOfOps2::LazySum)
+    return LazySum([SumOfOps1..., SumOfOps2...])
+end
+
+Base.:+(op1::LazySum, op2) = op1 + LazySum(op2)
+Base.:+(op1, op2::LazySum) = LazySum(op1) + op2
+
+Base.repeat(x::LazySum, args...) = LazySum(repeat.(x, args...))
diff --git a/src/states/mpsmultiline.jl b/src/states/mpsmultiline.jl
index a25eac80b..8d81f746d 100644
--- a/src/states/mpsmultiline.jl
+++ b/src/states/mpsmultiline.jl
@@ -71,6 +71,7 @@ site_type(::Type{Multiline{S}}) where {S} = site_type(S)
 bond_type(::Type{Multiline{S}}) where {S} = bond_type(S)
 site_type(st::Multiline) = site_type(typeof(st))
 bond_type(st::Multiline) = bond_type(typeof(st))
+VectorInterface.scalartype(::Multiline{T}) where {T} = scalartype(T)
 
 function TensorKit.dot(a::MPSMultiline, b::MPSMultiline; kwargs...)
     return sum(dot.(parent(a), parent(b); kwargs...))
diff --git a/test/algorithms.jl b/test/algorithms.jl
index 2c79badcd..c77e0778d 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -18,12 +18,12 @@ include("setup.jl")
         GradientGrassmann(; tol=tol, verbosity=verbosity),
     ]
 
-    H = force_planar(transverse_field_ising(; g=1.1))
+    H1 = force_planar(transverse_field_ising(; g=1.1))
 
     @testset "Infinite $i" for (i, alg) in enumerate(infinite_algs)
         L = alg isa IDMRG2 ? 2 : 1
         ψ₀ = repeat(InfiniteMPS([ℙ^2], [ℙ^16]), L)
-        H = repeat(force_planar(transverse_field_ising(; g=1.1)), L)
+        H = repeat(H1, L)
 
         v₀ = variance(ψ₀, H)
         ψ, envs, δ = find_groundstate(ψ₀, H, alg)
@@ -33,15 +33,35 @@ include("setup.jl")
         @test v₀ > v && v < 1e-2 # energy variance should be low
     end
 
+    Hlazy1 = LazySum([3 * H1, -1 * H1, 5.557 * H1])
+
+    @testset "LazySum Infinite $i" for (i, alg) in enumerate(infinite_algs)
+        L = alg isa IDMRG2 ? 2 : 1
+        ψ₀ = repeat(InfiniteMPS([ℙ^2], [ℙ^16]), L)
+        Hlazy = repeat(Hlazy1, L)
+
+        v₀ = variance(ψ₀, Hlazy)
+        ψ, envs, δ = find_groundstate(ψ₀, Hlazy, alg)
+        v = variance(ψ, Hlazy)
+
+        @test sum(δ) < 1e-3
+        @test v₀ > v && v < 1e-2 # energy variance should be low
+
+        ψ_nolazy, envs_nolazy, _ = find_groundstate(ψ₀, sum(Hlazy), alg)
+        @test expectation_value(ψ, Hlazy, envs) ≈
+            expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06
+    end
+
     finite_algs = [
         DMRG(; verbose=verbosity > 0),
         DMRG2(; verbose=verbosity > 0, trscheme=truncdim(10)),
         GradientGrassmann(; tol=tol, verbosity=verbosity),
     ]
 
+    H = force_planar(transverse_field_ising(; g=1.1))
+
     @testset "Finite $i" for (i, alg) in enumerate(finite_algs)
         ψ₀ = FiniteMPS(rand, ComplexF64, 10, ℙ^2, ℙ^10)
-        H = force_planar(transverse_field_ising(; g=1.1))
 
         v₀ = variance(ψ₀, H)
         ψ, envs, δ = find_groundstate(ψ₀, H, alg)
@@ -50,6 +70,23 @@ include("setup.jl")
         @test sum(δ) < 100 * tol
         @test v₀ > v && v < 1e-2 # energy variance should be low
     end
+
+    Hlazy = LazySum([H, H, H])
+
+    @testset "LazySum Finite $i" for (i, alg) in enumerate(finite_algs)
+        ψ₀ = FiniteMPS(rand, ComplexF64, 10, ℙ^2, ℙ^10)
+
+        v₀ = variance(ψ₀, Hlazy)
+        ψ, envs, δ = find_groundstate(ψ₀, Hlazy, alg)
+        v = variance(ψ, Hlazy)
+
+        @test sum(δ) < 100 * tol
+        @test v₀ > v && v < 1e-2 # energy variance should be low
+
+        ψ_nolazy, envs_nolazy, _ = find_groundstate(ψ₀, sum(Hlazy), alg)
+        @test expectation_value(ψ, Hlazy, envs) ≈
+            expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06
+    end
 end
 
 @testset "timestep" verbose = true begin
@@ -66,6 +103,14 @@ end
         @test sum(E₀) ≈ sum(E) atol = 1e-2
     end
 
+    Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H])
+
+    @testset "Finite LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs
+        ψ, envs = timestep(ψ₀, Hlazy, dt, alg)
+        E = expectation_value(ψ, Hlazy, envs)
+        @test (3 + 1.55 - 0.1) * sum(E₀) ≈ sum(E) atol = 1e-2
+    end
+
     H = repeat(force_planar(heisenberg_XXX(; spin=1)), 2)
     ψ₀ = InfiniteMPS([ℙ^3, ℙ^3], [ℙ^50, ℙ^50])
     E₀ = expectation_value(ψ₀, H)
@@ -75,6 +120,14 @@ end
         E = expectation_value(ψ, H, envs)
         @test sum(E₀) ≈ sum(E) atol = 1e-2
     end
+
+    Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H])
+
+    @testset "Infinite LazySum TDVP" begin
+        ψ, envs = timestep(ψ₀, Hlazy, dt, TDVP())
+        E = expectation_value(ψ, Hlazy, envs)
+        @test (3 + 1.55 - 0.1) * sum(E₀) ≈ sum(E) atol = 1e-2
+    end
 end
 
 @testset "leading_boundary" verbose = true begin
diff --git a/test/operators.jl b/test/operators.jl
index f1946f6c0..80411d7ed 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -71,6 +71,23 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1//2 => 10, 3//2 => 5, 5//2
     @test real(sum(expectation_value(ts2, h4))) >= 0
 end
 
+@testset "General LazySum of $(eltype(Os))" for Os in (
+    rand(ComplexF64, rand(1:10)),
+    map(i -> TensorMap(rand, ComplexF64, ℂ^13, ℂ^7), 1:rand(1:10)),
+    map(i -> TensorMap(rand, ComplexF64, ℂ^1 ⊗ ℂ^2, ℂ^3 ⊗ ℂ^4), 1:rand(1:10)),
+)
+    LazyOs = LazySum(Os)
+
+    #test user interface
+    summed = sum(Os)
+
+    @test sum(LazyOs) ≈ summed atol = 1 - 08
+
+    LazyOs_added = +(LazyOs, Os...)
+
+    @test sum(LazyOs_added) ≈ 2 * summed atol = 1 - 08
+end
+
 @testset "DenseMPO" for ham in (transverse_field_ising(), heisenberg_XXX(; spin=1))
     physical_space = ham.pspaces[1]
     ou = oneunit(physical_space)
@@ -81,3 +98,67 @@ end
 
     @test abs(dot(W * (W * ts), (W * W) * ts)) ≈ 1.0 atol = 1e-10
 end
+
+pspaces = (ℙ^4, Rep[U₁](0 => 2), Rep[SU₂](1 => 1, 2 => 1))
+vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 => 10, 3 => 5, 5 => 1))
+
+@testset "LazySum of (effective) Hamiltonian $(sectortype(pspace))" for (pspace, Dspace) in
+                                                                        zip(
+    pspaces, vspaces
+)
+    n = TensorMap(rand, ComplexF64, pspace, pspace)
+    n += n'
+    nn = TensorMap(rand, ComplexF64, pspace * pspace, pspace * pspace)
+    nn += nn'
+    nnn = TensorMap(rand, ComplexF64, pspace * pspace * pspace, pspace * pspace * pspace)
+    nnn += nnn'
+
+    H1 = repeat(MPOHamiltonian(n), 2)
+    H2 = repeat(MPOHamiltonian(nn), 2)
+    H3 = repeat(MPOHamiltonian(nnn), 2)
+    Hs = [H1, H2, H3]
+    summedH = LazySum(Hs)
+
+    Ψs = [
+        FiniteMPS(rand, ComplexF64, rand(3:2:20), pspace, Dspace),
+        InfiniteMPS([
+            TensorMap(rand, ComplexF64, Dspace * pspace, Dspace),
+            TensorMap(rand, ComplexF64, Dspace * pspace, Dspace),
+        ]),
+    ]
+
+    @testset "LazySum $(Ψ isa FiniteMPS ? "F" : "Inf")initeMPS" for Ψ in Ψs
+        Envs = map(H -> environments(Ψ, H), Hs)
+        summedEnvs = environments(Ψ, summedH)
+
+        expval = sum(zip(Hs, Envs)) do (H, Env)
+            expectation_value(Ψ, H, Env)
+        end
+        expval1 = expectation_value(Ψ, sum(summedH))
+        expval2 = expectation_value(Ψ, summedH, summedEnvs)
+        expval3 = expectation_value(Ψ, summedH)
+        @test expval ≈ expval1
+        @test expval ≈ expval2
+        @test expval ≈ expval3
+
+        # test derivatives
+        summedhct = MPSKit.∂∂C(1, Ψ, summedH, summedEnvs)
+        sum1 = sum(zip(Hs, Envs)) do (H, env)
+            MPSKit.∂∂C(1, Ψ, H, env)(Ψ.CR[1])
+        end
+        @test summedhct(Ψ.CR[1]) ≈ sum1
+
+        summedhct = MPSKit.∂∂AC(1, Ψ, summedH, summedEnvs)
+        sum2 = sum(zip(Hs, Envs)) do (H, env)
+            MPSKit.∂∂AC(1, Ψ, H, env)(Ψ.AC[1])
+        end
+        @test summedhct(Ψ.AC[1]) ≈ sum2
+
+        v = MPSKit._transpose_front(Ψ.AC[1]) * MPSKit._transpose_tail(Ψ.AR[2])
+        summedhct = MPSKit.∂∂AC2(1, Ψ, summedH, summedEnvs)
+        sum3 = sum(zip(Hs, Envs)) do (H, env)
+            MPSKit.∂∂AC2(1, Ψ, H, env)(v)
+        end
+        @test summedhct(v) ≈ sum3
+    end
+end