diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl
index f905c5405..3c564cde2 100644
--- a/src/algorithms/approximate/approximate.jl
+++ b/src/algorithms/approximate/approximate.jl
@@ -32,14 +32,13 @@ function approximate(
         envs = environments(ψ, toapprox)
     )
     envs′ = Multiline([envs])
-    multi, envs = approximate(
+    multi, envs, δ = approximate(
         convert(MultilineMPS, ψ),
         (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])),
-        algorithm,
-        envs′
+        algorithm, envs′
     )
     ψ = convert(InfiniteMPS, multi)
-    return ψ, envs
+    return ψ, envs, δ
 end
 
 # dispatch to in-place method
@@ -56,12 +55,11 @@ function approximate(
         algorithm::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox)
     )
     envs′ = Multiline([envs])
-    multi, envs = approximate(
+    multi, envs, δ = approximate(
         convert(MultilineMPS, ψ),
         (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])),
-        algorithm,
-        envs′
+        algorithm, envs′
     )
     ψ = convert(InfiniteMPS, multi)
-    return ψ, envs
+    return ψ, envs, δ
 end
diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl
index 7ad26e706..095e260f6 100644
--- a/src/algorithms/changebonds/svdcut.jl
+++ b/src/algorithms/changebonds/svdcut.jl
@@ -29,9 +29,9 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::SvdCut; normalize::Bool = true
     for i in (length(ψ) - 1):-1:1
         U, S, V, = tsvd(ψ.C[i]; trunc = alg.trscheme, alg = alg.alg_svd)
         AL′ = ψ.AL[i] * U
-        ψ.AC[i] = (AL′, complex(S))
+        ψ.AC[i] = (AL′, S)
         AR′ = _transpose_front(V * _transpose_tail(ψ.AR[i + 1]))
-        ψ.AC[i + 1] = (complex(S), AR′)
+        ψ.AC[i + 1] = (S, AR′)
     end
     return normalize ? normalize!(ψ) : ψ
 end
@@ -87,7 +87,7 @@ function changebonds(ψ::MultilineMPS, alg::SvdCut)
     return Multiline(map(x -> changebonds(x, alg), ψ.data))
 end
 function changebonds(ψ::InfiniteMPS, alg::SvdCut)
-    copied = complex.(ψ.AL)
+    copied = copy.(ψ.AL)
     ncr = ψ.C[1]
 
     for i in 1:length(ψ)
@@ -103,7 +103,7 @@ function changebonds(ψ::InfiniteMPS, alg::SvdCut)
     ψ = if space(ncr, 1) != space(copied[1], 1)
         InfiniteMPS(copied)
     else
-        C₀ = ncr isa TensorMap ? complex(ncr) : TensorMap(complex(ncr))
+        C₀ = ncr isa TensorMap ? ncr : TensorMap(ncr)
         InfiniteMPS(copied, C₀)
     end
     return normalize!(ψ)
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index e8ac68738..dfce5bc35 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -37,53 +37,42 @@ function expectation_value end
 function expectation_value(ψ::AbstractMPS, (inds, O)::Pair)
     @boundscheck foreach(Base.Fix1(Base.checkbounds, ψ), inds)
 
-    sites, local_mpo = instantiate_operator(physicalspace(ψ), inds => O)
-    @assert _firstspace(first(local_mpo)) == oneunit(_firstspace(first(local_mpo))) ==
-        dual(_lastspace(last(local_mpo)))
-    for (site, o) in zip(sites, local_mpo)
-        if o isa MPOTensor
-            physicalspace(ψ, site) == physicalspace(o) ||
-                throw(SpaceMismatch("physical space does not match at site $site"))
-        end
+    # special cases that can be handled more efficiently
+    if length(inds) == 1
+        return local_expectation_value1(ψ, inds[1], O)
+    elseif length(inds) == 2 && eltype(inds) <: Integer && inds[1] + 1 == inds[2]
+        return local_expectation_value2(ψ, inds[1], O)
     end
 
-    Ut = fill_data!(similar(local_mpo[1], _firstspace(first(local_mpo))), one)
-
-    # some special cases that avoid using transfer matrices
-    if length(sites) == 1
-        AC = ψ.AC[sites[1]]
-        E = @plansor conj(AC[4 5; 6]) * conj(Ut[1]) * local_mpo[1][1 5; 3 2] * Ut[2] *
-            AC[4 3; 6]
-    elseif length(sites) == 2 && (sites[1] + 1 == sites[2])
-        AC = ψ.AC[sites[1]]
-        AR = ψ.AR[sites[2]]
-        O1, O2 = local_mpo
-        E = @plansor conj(AC[4 5; 10]) * conj(Ut[1]) * O1[1 5; 3 8] * AC[4 3; 6] *
-            conj(AR[10 9; 11]) * Ut[2] * O2[8 9; 7 2] * AR[6 7; 11]
-    else
-        # generic case: write as Vl * T^N * Vr
-        # left side
-        T = storagetype(site_type(ψ))
-        @plansor Vl[-1 -2; -3] := isomorphism(
-            T, left_virtualspace(ψ, sites[1]), left_virtualspace(ψ, sites[1])
-        )[-1; -3] *
-            conj(Ut[-2])
-
-        # middle
-        M = foldl(zip(sites, local_mpo); init = Vl) do v, (site, o)
-            if o isa Number
-                return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o)
-            else
-                return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site])
-            end
-        end
+    # generic case: convert to MPO and write as Vl * T^N * Vr
+    sites, local_mpo = instantiate_operator(ψ, inds => O)
+
+    # left side
+    Vl = insertrightunit(l_LL(ψ, sites[1]), 1; dual = true)
 
-        # right side
-        E = @plansor M[1 2; 3] * Ut[2] * ψ.C[sites[end]][3; 4] *
-            conj(ψ.C[sites[end]][1; 4])
+    # middle
+    M = foldl(zip(sites, local_mpo); init = Vl) do v, (site, o)
+        if o isa Number
+            return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o)
+        else
+            return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site])
+        end
     end
 
-    return E / norm(ψ)^2
+    # right side
+    E = @plansor removeunit(M, 2)[1; 2] * ψ.C[sites[end]][2; 3] * conj(ψ.C[sites[end]][1; 3])
+    return E / dot(ψ, ψ)
+end
+
+function local_expectation_value1(ψ::AbstractMPS, site, O)
+    E = contract_mpo_expval1(ψ.AC[site], O, ψ.AC[site])
+    return E / dot(ψ, ψ)
+end
+function local_expectation_value2(ψ::AbstractMPS, site, O)
+    AC = ψ.AC[site]
+    AR = ψ.AR[site + 1]
+    E = contract_mpo_expval2(AC, AR, O, AC, AR)
+    return E / dot(ψ, ψ)
 end
 
 # MPOHamiltonian
@@ -93,6 +82,40 @@ function contract_mpo_expval(
     )
     return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6])
 end
+# generic fallback
+function contract_mpo_expval(AC, GL, O, GR, ACbar = AC)
+    return dot(ACbar, MPO_AC_Hamiltonian(GL, O, GR) * AC)
+end
+
+function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor = AC)
+    numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator"))
+    return @plansor conj(ACbar[2 3; 4]) * O[3; 1] * AC[2 1; 4]
+end
+function contract_mpo_expval1(
+        AC::GenericMPSTensor{S, 3}, O::AbstractTensorMap{<:Any, S},
+        ACbar::GenericMPSTensor{S, 3} = AC
+    ) where {S}
+    numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator"))
+    return @plansor conj(ACbar[2 3 4; 5]) * O[3; 1] * AC[2 1 4; 5]
+end
+
+function contract_mpo_expval2(
+        A1::MPSTensor, A2::MPSTensor, O::AbstractTensorMap,
+        A1bar::MPSTensor = A1, A2bar::MPSTensor = A2
+    )
+    numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator"))
+    return @plansor conj(A1bar[4 5; 6]) * conj(A2bar[6 7; 8]) * O[5 7; 2 3] * A1[4 2; 1] *
+        A2[1 3; 8]
+end
+function contract_mpo_expval2(
+        A1::GenericMPSTensor{S, 3}, A2::GenericMPSTensor{S, 3},
+        O::AbstractTensorMap{<:Any, S},
+        A1bar::GenericMPSTensor{S, 3} = A1, A2bar::GenericMPSTensor{S, 3} = A2
+    ) where {S}
+    numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator"))
+    return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ[9 6; 1 2] *
+        τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7]
+end
 
 function expectation_value(
         ψ::FiniteMPS, H::FiniteMPOHamiltonian,
@@ -177,3 +200,12 @@ function expectation_value(
     n = norm(ψ.AC[end])^2
     return sum(ens) / (n * length(ψ))
 end
+
+# Density matrices
+# ----------------
+function expectation_value(ρ::FiniteMPO, args...)
+    return expectation_value(convert(FiniteMPS, ρ), args...)
+end
+function expectation_value(ρ::InfiniteMPO, args...)
+    return expectation_value(convert(InfiniteMPS, ρ), args...)
+end
diff --git a/src/algorithms/timestep/integrators.jl b/src/algorithms/timestep/integrators.jl
index 7d19591d8..7b2983c34 100644
--- a/src/algorithms/timestep/integrators.jl
+++ b/src/algorithms/timestep/integrators.jl
@@ -17,8 +17,12 @@ function integrate end
 _eval_t(f, t::Number) = Base.Fix2(f, t)
 _eval_x(f, x) = Base.Fix1(f, x)
 
-function integrate(f, y₀, t::Number, dt::Number, alg::Union{Arnoldi, Lanczos})
-    y, convhist = exponentiate(_eval_t(f, t + dt / 2), -1im * dt, y₀, alg)
+function integrate(
+        f, y₀, t::Number, dt::Number, alg::Union{Arnoldi, Lanczos};
+        imaginary_evolution::Bool = false
+    )
+    dt′ = imaginary_evolution ? -dt : -1im * dt
+    y, convhist = exponentiate(_eval_t(f, t + dt / 2), dt′, y₀, alg)
     convhist.converged == 0 &&
         @warn "integrator failed to converge: normres = $(convhist.normres)"
     return y
diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl
index cbad1ee13..d9cd70247 100644
--- a/src/algorithms/timestep/taylorcluster.jl
+++ b/src/algorithms/timestep/taylorcluster.jl
@@ -28,10 +28,11 @@ First order Taylor expansion for a time-evolution MPO.
 const WI = TaylorCluster(; N = 1, extension = false, compression = false)
 
 function make_time_mpo(
-        H::MPOHamiltonian, dt::Number, alg::TaylorCluster; tol = eps(real(scalartype(H)))
+        H::MPOHamiltonian, dt::Number, alg::TaylorCluster;
+        tol = eps(real(scalartype(H))), imaginary_evolution::Bool = false
     )
     N = alg.N
-    τ = -1im * dt
+    τ = imaginary_evolution ? -dt : -1im * dt
 
     # start with H^N
     H_n = H^N
@@ -180,7 +181,7 @@ end
 # Hack to treat FiniteMPOhamiltonians as Infinite
 function make_time_mpo(
         H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster;
-        tol = eps(real(scalartype(H)))
+        tol = eps(real(scalartype(H))), imaginary_evolution::Bool = false
     )
     H′ = copy(parent(H))
 
@@ -201,7 +202,7 @@ function make_time_mpo(
     H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1]
     H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end]
 
-    mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol)
+    mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol, imaginary_evolution)
 
     # Impose boundary conditions
     mpo_fin = open_boundary_conditions(mpo, length(H))
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index b711d022c..e372c4fbd 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -26,11 +26,15 @@ $(TYPEDFIELDS)
 end
 
 function timestep(
-        ψ_::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP,
-        envs::AbstractMPSEnvironments = environments(ψ_, H);
-        leftorthflag = true
+        ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP,
+        envs::AbstractMPSEnvironments = environments(ψ, H);
+        leftorthflag = true, imaginary_evolution::Bool = false
     )
-    ψ = complex(ψ_)
+    # convert state to complex if necessary
+    if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt))
+        return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution)
+    end
+
     temp_ACs = similar(ψ.AC)
     temp_Cs = similar(ψ.C)
 
@@ -38,24 +42,30 @@ function timestep(
     if scheduler isa SerialScheduler
         temp_ACs = tmap!(temp_ACs, 1:length(ψ); scheduler) do loc
             Hac = AC_hamiltonian(loc, ψ, H, ψ, envs)
-            return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator)
+            return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator; imaginary_evolution)
         end
         temp_Cs = tmap!(temp_Cs, 1:length(ψ); scheduler) do loc
             Hc = C_hamiltonian(loc, ψ, H, ψ, envs)
-            return integrate(Hc, ψ.C[loc], t, dt, alg.integrator)
+            return integrate(Hc, ψ.C[loc], t, dt, alg.integrator; imaginary_evolution)
         end
     else
         @sync begin
             Threads.@spawn begin
                 temp_ACs = tmap!(temp_ACs, 1:length(ψ); scheduler) do loc
                     Hac = AC_hamiltonian(loc, ψ, H, ψ, envs)
-                    return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator)
+                    return integrate(
+                        Hac, ψ.AC[loc], t, dt, alg.integrator;
+                        imaginary_evolution
+                    )
                 end
             end
             Threads.@spawn begin
                 temp_Cs = tmap!(temp_Cs, 1:length(ψ); scheduler) do loc
                     Hc = C_hamiltonian(loc, ψ, H, ψ, envs)
-                    return integrate(Hc, ψ.C[loc], t, dt, alg.integrator)
+                    return integrate(
+                        Hc, ψ.C[loc], t, dt, alg.integrator;
+                        imaginary_evolution
+                    )
                 end
             end
         end
@@ -76,34 +86,47 @@ end
 
 function timestep!(
         ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP,
-        envs::AbstractMPSEnvironments = environments(ψ, H)
+        envs::AbstractMPSEnvironments = environments(ψ, H);
+        imaginary_evolution::Bool = false
     )
 
     # sweep left to right
     for i in 1:(length(ψ) - 1)
         Hac = AC_hamiltonian(i, ψ, H, ψ, envs)
-        ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator)
+        ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution)
 
         Hc = C_hamiltonian(i, ψ, H, ψ, envs)
-        ψ.C[i] = integrate(Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator)
+        ψ.C[i] = integrate(
+            Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator;
+            imaginary_evolution
+        )
     end
 
     # edge case
     Hac = AC_hamiltonian(length(ψ), ψ, H, ψ, envs)
-    ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator)
+    ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator; imaginary_evolution)
 
     # sweep right to left
     for i in length(ψ):-1:2
         Hac = AC_hamiltonian(i, ψ, H, ψ, envs)
-        ψ.AC[i] = integrate(Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator)
+        ψ.AC[i] = integrate(
+            Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator;
+            imaginary_evolution
+        )
 
         Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs)
-        ψ.C[i - 1] = integrate(Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator)
+        ψ.C[i - 1] = integrate(
+            Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator;
+            imaginary_evolution
+        )
     end
 
     # edge case
     Hac = AC_hamiltonian(1, ψ, H, ψ, envs)
-    ψ.AC[1] = integrate(Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator)
+    ψ.AC[1] = integrate(
+        Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator;
+        imaginary_evolution
+    )
 
     return ψ, envs
 end
@@ -143,14 +166,15 @@ end
 
 function timestep!(
         ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2,
-        envs::AbstractMPSEnvironments = environments(ψ, H)
+        envs::AbstractMPSEnvironments = environments(ψ, H);
+        imaginary_evolution::Bool = false
     )
 
     # sweep left to right
     for i in 1:(length(ψ) - 1)
         ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])
         Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs)
-        ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator)
+        ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution)
 
         nal, nc, nar = tsvd!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd)
         ψ.AC[i] = (nal, complex(nc))
@@ -158,7 +182,10 @@ function timestep!(
 
         if i != (length(ψ) - 1)
             Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs)
-            ψ.AC[i + 1] = integrate(Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator)
+            ψ.AC[i + 1] = integrate(
+                Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator;
+                imaginary_evolution
+            )
         end
     end
 
@@ -166,7 +193,7 @@ function timestep!(
     for i in length(ψ):-1:2
         ac2 = _transpose_front(ψ.AL[i - 1]) * _transpose_tail(ψ.AC[i])
         Hac2 = AC2_hamiltonian(i - 1, ψ, H, ψ, envs)
-        ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator)
+        ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution)
 
         nal, nc, nar = tsvd!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd)
         ψ.AC[i - 1] = (nal, complex(nc))
@@ -174,20 +201,23 @@ function timestep!(
 
         if i != 2
             Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs)
-            ψ.AC[i - 1] = integrate(Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator)
+            ψ.AC[i - 1] = integrate(
+                Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator;
+                imaginary_evolution
+            )
         end
     end
 
     return ψ, envs
 end
 
-#copying version
+# copying version
 function timestep(
         ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,
         alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...;
-        kwargs...
+        imaginary_evolution::Bool = false, kwargs...
     )
-    isreal = scalartype(ψ) <: Real
+    isreal = (scalartype(ψ) <: Real && !imaginary_evolution)
     ψ′ = isreal ? complex(ψ) : copy(ψ)
     if length(envs) != 0 && isreal
         @warn "Currently cannot reuse real environments for complex evolution"
@@ -198,5 +228,5 @@ function timestep(
         @assert length(envs) == 0 "Invalid signature"
         envs′ = environments(ψ′, H)
     end
-    return timestep!(ψ′, H, time, timestep, alg, envs′; kwargs...)
+    return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, kwargs...)
 end
diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl
index a9b06eeda..a0918ad9f 100644
--- a/src/algorithms/timestep/time_evolve.jl
+++ b/src/algorithms/timestep/time_evolve.jl
@@ -5,19 +5,29 @@
 Time-evolve the initial state `ψ₀` with Hamiltonian `H` over a given time span by stepping
 through each of the time points obtained by iterating t_span.
 
-# Arguments
+## Arguments
+
 - `ψ₀::AbstractMPS`: initial state
 - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent).
 - `t_span::AbstractVector{<:Number}`: time points over which the time evolution is stepped
 - `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref).
 - `[envs]`: MPS environment manager
+
+## Keyword Arguments
+
+- `verbosity::Int=0`: verbosity level for logging
+- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step
+    instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
+    function to compute the ground state of a Hamiltonian, or to compute finite-temperature
+    properties of a system.
 """
 function time_evolve end, function time_evolve! end
 
 for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!))
     @eval function $time_evolve(
-            ψ, H, t_span::AbstractVector{<:Number}, alg, envs = environments(ψ, H);
-            verbosity::Int = 0
+            ψ, H, t_span::AbstractVector{<:Number}, alg,
+            envs = environments(ψ, H);
+            verbosity::Int = 0, imaginary_evolution::Bool = false
         )
         log = IterLog("TDVP")
         LoggingExtras.withlevel(; verbosity) do
@@ -26,7 +36,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time
                 t = t_span[iter]
                 dt = t_span[iter + 1] - t
 
-                ψ, envs = $timestep(ψ, H, t, dt, alg, envs)
+                ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution)
                 ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)}
 
                 @infov 3 logiter!(log, iter, 0, t)
@@ -44,18 +54,33 @@ end
 Time-step the state `ψ₀` with Hamiltonian `H` over a given time step `dt` at time `t`,
 solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``.
 
-# Arguments
+## Arguments
+
 - `ψ₀::AbstractMPS`: initial state
 - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent).
 - `t::Number`: starting time of time-step
 - `dt::Number`: time-step magnitude
 - `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref).
 - `[envs]`: MPS environment manager
+
+## Keyword Arguments
+
+- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step
+    instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
+    function to compute the ground state of a Hamiltonian, or to compute finite-temperature
+    properties of a system.
 """
 function timestep end, function timestep! end
 
 @doc """
-    make_time_mpo(H::MPOHamiltonian, dt::Number, alg) -> O::MPO
+    make_time_mpo(H::MPOHamiltonian, dt::Number, alg; kwargs...) -> O::MPO
 
 Construct an `MPO` that approximates ``\\exp(-iHdt)``.
+
+## Keyword Arguments
+
+- `imaginary_evolution::Bool=false`: if true, the time evolution operator is constructed
+    with an imaginary time step instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``).
+    This can be useful for using this function to compute the ground state of a Hamiltonian,
+    or to compute finite-temperature properties of a system.
 """ make_time_mpo
diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl
index cd02352d5..3cf65cad6 100644
--- a/src/algorithms/timestep/wii.jl
+++ b/src/algorithms/timestep/wii.jl
@@ -19,13 +19,17 @@ $(TYPEDFIELDS)
     maxiter::Int = Defaults.maxiter
 end
 
-function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII)
+function make_time_mpo(
+        H::InfiniteMPOHamiltonian, dt::Number, alg::WII;
+        imaginary_evolution::Bool = false
+    )
     WA = H.A
     WB = H.B
     WC = H.C
     WD = H.D
 
-    δ = dt * (-1im)
+    # needs to be complex if negative because of square root
+    δ = imaginary_evolution ? -complex(dt) : (-1im * dt)
     exp_alg = Arnoldi(; alg.tol, alg.maxiter)
 
     Wnew = map(1:length(H)) do i
@@ -83,7 +87,10 @@ function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII)
 end
 
 # Hack to treat FiniteMPOhamiltonians as Infinite
-function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII)
+function make_time_mpo(
+        H::FiniteMPOHamiltonian, dt::Number, alg::WII;
+        imaginary_evolution::Bool = false
+    )
     H′ = copy(parent(H))
 
     V_left = left_virtualspace(H[1])
@@ -103,9 +110,9 @@ function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII)
     H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1]
     H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end]
 
-    mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol)
+    mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; imaginary_evolution)
 
     # Impose boundary conditions
     mpo_fin = open_boundary_conditions(mpo, length(H))
-    return remove_orphans!(mpo_fin; tol)
+    return remove_orphans!(mpo_fin; alg.tol)
 end
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index 7fcc7c6e8..bb7932414 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -194,6 +194,26 @@ function SparseBlockTensorMap(W::JordanMPOTensor)
     return W′
 end
 
+for f in (:real, :complex)
+    @eval function Base.$f(W::JordanMPOTensor)
+        E = $f(scalartype(W))
+        W′ = JordanMPOTensor{E}(undef, space(W))
+        for (I, v) in nonzero_pairs(W.A)
+            W′.A[I] = $f(v)
+        end
+        for (I, v) in nonzero_pairs(W.B)
+            W′.B[I] = $f(v)
+        end
+        for (I, v) in nonzero_pairs(W.C)
+            W′.C[I] = $f(v)
+        end
+        for (I, v) in nonzero_pairs(W.D)
+            W′.D[I] = $f(v)
+        end
+        return W′
+    end
+end
+
 # Indexing
 # --------
 
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index 8bc220f58..c152c578a 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -319,13 +319,13 @@ function TensorKit.dot(
     randomize!(ρ₀)
 
     val, = fixedpoint(
-        TransferMatrix(ket.AL, parent(mpo), bra.AL), ρ₀, :LM; ishermitian,
-        krylovdim, kwargs...
+        flip(TransferMatrix(ket.AL, parent(mpo), bra.AL)), ρ₀, :LM;
+        ishermitian, krylovdim, kwargs...
     )
     return val
 end
 
-function TensorKit.dot(mpo₁::FiniteMPO{TO}, mpo₂::FiniteMPO{TO}) where {TO <: MPOTensor}
+function TensorKit.dot(mpo₁::FiniteMPO{<:MPOTensor}, mpo₂::FiniteMPO{<:MPOTensor})
     length(mpo₁) == length(mpo₂) || throw(ArgumentError("dimension mismatch"))
     N = length(mpo₁)
     Nhalf = N ÷ 2
@@ -344,7 +344,7 @@ function TensorKit.dot(mpo₁::FiniteMPO{TO}, mpo₂::FiniteMPO{TO}) where {TO <
     return @plansor ρ_left[1; 2] * ρ_right[2; 1]
 end
 
-function LinearAlgebra.tr(mpo::MPO)
+function LinearAlgebra.tr(mpo::FiniteMPO)
     N = length(mpo)
     # @assert BraidingStyle(sectortype(mpo)) isa SymmetricBraiding
     Nhalf = N ÷ 2
@@ -367,6 +367,14 @@ function LinearAlgebra.tr(mpo::MPO)
 
     return @plansor ρ_left[(); 1] * ρ_right[1; ()]
 end
+function LinearAlgebra.tr(mpo::InfiniteMPO; ishermitian = false, kwargs...)
+    T = prod(parent(mpo)) do O
+        @plansor Tᵢ[-1; -2] := O[-1 3; 1 2] * τ[1 2; -2 3]
+    end
+    v₀ = randn(scalartype(mpo), domain(T))
+    val, = fixedpoint(Base.Fix1(*, T), v₀, :LM; kwargs...)
+    return val
+end
 
 function Base.isapprox(
         mpo₁::MPO, mpo₂::MPO;
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 0727f4ca3..476983658 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -299,11 +299,15 @@ function _split_mpoham_types(W::Matrix)::Type{<:MPOTensor}
 end
 
 """
+    instantiate_operator(state, O::Pair)
     instantiate_operator(lattice::AbstractArray{<:VectorSpace}, O::Pair)
 
-Instantiate a local operator `O` on a lattice `lattice` as a vector of MPO tensors, and a
-vector of linear site indices.
+Instantiate a local operator `O` for a `state` or `lattice` as a vector of MPO tensors, and
+a vector of linear site indices.
 """
+function instantiate_operator(state::AbstractMPS, O::Pair)
+    return instantiate_operator(physicalspace(state), O)
+end
 function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O)::Pair)
     inds = inds′ isa Int ? tuple(inds′) : inds′
     mpo = O isa FiniteMPO ? O : FiniteMPO(O)
@@ -320,7 +324,10 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O
     i = 1
     for j in first(indices):last(indices)
         if j == indices[i]
-            @assert space(operators[i], 2) == lattice[j] "operator does not fit into the given Hilbert space: $(space(operators[i], 2)) ≠ $(lattice[j])"
+            # TODO: fix this check for density matrices
+            if !(eltype(lattice) <: ProductSpace) && physicalspace(operators[i]) != lattice[j]
+                throw(SpaceMismatch("physical space does not match at site $j"))
+            end
             push!(local_mpo, operators[i])
             i += 1
         else
@@ -629,19 +636,8 @@ end
 # Base.complex(H::MPOHamiltonian) = MPOHamiltonian(map(complex, parent(H)))
 function Base.complex(H::MPOHamiltonian)
     scalartype(H) <: Complex && return H
-    Ws = map(parent(H)) do W
-        W′ = jordanmpotensortype(spacetype(W), complex(scalartype(W)))
-        W′[1] = W[1]
-        W′[end] = W[end]
-        for (I, v) in nonzero_pairs(W)
-            if v isa BraidingTensor
-                W′[I] = BraidingTensor{scalartype(W′)}(space(v), v.adjoint)
-            else
-                W′[I] = complex(v)
-            end
-        end
-    end
-    return MPOHamiltonian(H)
+    Ws = map(complex, parent(H))
+    return MPOHamiltonian(Ws)
 end
 
 function Base.similar(H::MPOHamiltonian, ::Type{O}, L::Int) where {O <: MPOTensor}
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index 1300ceaa4..cc877603e 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -626,9 +626,11 @@ TensorKit.normalize(ψ::FiniteMPS) = normalize!(copy(ψ))
 Fixedpoints
 ===========================================================================================#
 
-function r_RR(ψ::FiniteMPS{T}) where {T}
-    return isomorphism(storagetype(T), domain(ψ.AR[end]), domain(ψ.AR[end]))
+function r_RR(ψ::FiniteMPS, site::Int = length(ψ))
+    Vr = right_virtualspace(ψ.AR[site])
+    return isomorphism(storagetype(site_type(ψ)), Vr ← Vr)
 end
-function l_LL(ψ::FiniteMPS{T}) where {T}
-    return isomorphism(storagetype(T), space(ψ.AL[1], 1), space(ψ.AL[1], 1))
+function l_LL(ψ::FiniteMPS, site::Int = 1)
+    Vl = left_virtualspace(ψ.AL[site])
+    return isomorphism(storagetype(site_type(ψ)), Vl ← Vl)
 end
diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl
index b92bd9c40..c9a4cf906 100644
--- a/src/states/infinitemps.jl
+++ b/src/states/infinitemps.jl
@@ -193,8 +193,10 @@ function InfiniteMPS(AL::AbstractVector{<:GenericMPSTensor}, C₀::MPSBondTensor
     # initialize tensor storage
     AC = similar.(AL)
     AR = similar.(AL)
-    C = similar(AR, typeof(C₀))
-    ψ = InfiniteMPS{eltype(AL), eltype(C)}(AL, AR, C, AC)
+    T = TensorOperations.promote_contract(scalartype(AL), scalartype(C₀))
+    TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false)
+    C = similar(AR, TC)
+    ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC)
 
     # gaugefix the MPS
     gaugefix!(ψ, AL, C₀; order = :R, kwargs...)
@@ -209,8 +211,10 @@ function InfiniteMPS(C₀::MPSBondTensor, AR::AbstractVector{<:GenericMPSTensor}
     # initialize tensor storage
     AC = similar.(AR)
     AL = similar.(AR)
-    C = similar(AR, typeof(C₀))
-    ψ = InfiniteMPS{eltype(AL), eltype(C)}(AL, AR, C, AC)
+    T = TensorOperations.promote_contract(eltype(AR), eltype(C₀))
+    TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false)
+    C = similar(AR, TC)
+    ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC)
 
     # gaugefix the MPS
     gaugefix!(ψ, AR, C₀; order = :L, kwargs...)
@@ -237,14 +241,19 @@ Base.size(ψ::InfiniteMPS, args...) = size(ψ.AL, args...)
 Base.length(ψ::InfiniteMPS) = length(ψ.AL)
 Base.eltype(ψ::InfiniteMPS) = eltype(typeof(ψ))
 Base.eltype(::Type{<:InfiniteMPS{A}}) where {A} = A
+
 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)
+    ψ.AL .= _copy!!.(ψ.AL, ϕ.AL)
+    ψ.AR .= _copy!!.(ψ.AR, ϕ.AR)
+    ψ.AC .= _copy!!.(ψ.AC, ϕ.AC)
+    ψ.C .= _copy!!.(ψ.C, ϕ.C)
     return ψ
 end
+# possible in-place copy
+function _copy!!(dst::AbstractTensorMap, src::AbstractTensorMap)
+    return space(dst) == space(src) ? copy!(dst, src) : copy(src)
+end
 
 function Base.complex(ψ::InfiniteMPS)
     scalartype(ψ) <: Complex && return ψ
diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl
index c3abe3325..fd810e95d 100644
--- a/src/states/orthoview.jl
+++ b/src/states/orthoview.jl
@@ -160,10 +160,18 @@ function Base.setindex!(
 
     a, b = vec
     return if isa(a, MPSBondTensor) #c/ar
-        setindex!(v.parent.Cs, a, i)
+        if !(scalartype(parent(v)) <: Real) && (scalartype(a) <: Real)
+            setindex!(v.parent.Cs, complex(a), i)
+        else
+            setindex!(v.parent.Cs, a, i)
+        end
         setindex!(v.parent.ARs, b, i)
     elseif isa(b, MPSBondTensor) #al/c
-        setindex!(v.parent.Cs, b, i + 1)
+        if !(scalartype(parent(v)) <: Real) && (scalartype(b) <: Real)
+            setindex!(v.parent.Cs, complex(b), i + 1)
+        else
+            setindex!(v.parent.Cs, b, i + 1)
+        end
         setindex!(v.parent.ALs, a, i)
     else
         throw(ArgumentError("invalid value types"))
diff --git a/src/transfermatrix/transfer.jl b/src/transfermatrix/transfer.jl
index 15cc2237a..aa3d10ed7 100644
--- a/src/transfermatrix/transfer.jl
+++ b/src/transfermatrix/transfer.jl
@@ -81,6 +81,21 @@ function transfer_right(v::MPSTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) wher
     check_unambiguous_braiding(space(v, 2))
     return @plansor v[-1 -2; -3] := A[-1 2; 1] * τ[-2 4; 2 3] * conj(Ab[-3 4; 5]) * v[1 3; 5]
 end
+# similar for Matrix Product Density Operators
+function transfer_left(
+        v::MPSTensor{S}, A::GenericMPSTensor{S, 3}, Ab::GenericMPSTensor{S, 3}
+    ) where {S}
+    check_unambiguous_braiding(space(v, 2))
+    return @plansor v[-1 -2; -3] ≔ v[1 3; 8] * A[8 7 6; -3] * τ[3 2; 7 5] * τ[5 4; 6 -2] *
+        conj(Ab[1 2 4; -1])
+end
+function transfer_right(
+        v::MPSTensor{S}, A::GenericMPSTensor{S, 3}, Ab::GenericMPSTensor{S, 3}
+    ) where {S}
+    check_unambiguous_braiding(space(v, 2))
+    return @plansor v[-1 -2; -3] ≔ A[-1 4 2; 1] * τ[-2 6; 4 5] * τ[5 7; 2 3] *
+        conj(Ab[-3 6 7; 8]) * v[1 3; 8]
+end
 
 # the transfer operation with a utility leg in both the domain and codomain is also ill defined - only due to the codomain utility space
 function transfer_left(v::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S}
diff --git a/test/algorithms.jl b/test/algorithms.jl
index b009c121b..1589b586d 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -950,36 +950,94 @@ module TestAlgorithms
     end
 
     @testset "Finite temperature methods" begin
-        L = 6
-        H = transverse_field_ising(; L)
-        trscheme = truncdim(20)
-        verbosity = 1
-        beta = 0.1
+        imaginary_evolution = true
+        @testset "Finite-size" begin
+            L = 6
+            H = transverse_field_ising(; L)
+            trscheme = truncdim(20)
+            verbosity = 1
+            beta = 0.1
+
+            # exact diagonalization
+            H_dense = convert(TensorMap, H)
+            Z_dense_1 = tr(exp(-beta * H_dense))^(1 / L)
+            Z_dense_2 = tr(exp(-2beta * H_dense))^(1 / L)
+
+            # taylor cluster
+            rho_taylor_1 = make_time_mpo(H, beta, TaylorCluster(; N = 2); imaginary_evolution)
+            Z_taylor_1 = tr(rho_taylor_1)^(1 / L)
+            @test Z_taylor_1 ≈ Z_dense_1 atol = 1.0e-2
+            Z_taylor_2 = real(dot(rho_taylor_1, rho_taylor_1))^(1 / L)
+            @test Z_taylor_2 ≈ Z_dense_2 atol = 1.0e-2
+
+            E_x_taylor = @constinferred expectation_value(rho_taylor_1, 1 => S_x())
+            E_xx_taylor = @constinferred expectation_value(rho_taylor_1, (1, 2) => S_xx())
+
+            # WII
+            rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution)
+            Z_wii = tr(rho_wii)^(1 / L)
+            @test Z_wii ≈ Z_dense_1 atol = 1.0e-2
+            @test expectation_value(rho_wii, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2
+            @test expectation_value(rho_wii, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2
+
+            # MPO multiplication
+            rho_mps = convert(FiniteMPS, rho_taylor_1)
+            rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity))
+            Z_mpomul = tr(convert(FiniteMPO, rho_mps))^(1 / L)
+            @test Z_mpomul ≈ Z_dense_2 atol = 1.0e-2
+
+            # TDVP
+            rho_0 = MPSKit.infinite_temperature_density_matrix(H)
+            rho_0_mps = convert(FiniteMPS, rho_0)
+            rho_mps, = timestep(rho_0_mps, H, 0.0, beta, TDVP2(; trscheme); imaginary_evolution)
+            Z_tdvp = real(dot(rho_mps, rho_mps))^(1 / L)
+            @test Z_tdvp ≈ Z_dense_2 atol = 1.0e-2
+
+            @test expectation_value(rho_0_mps, 1 => S_x()) ≈ 0
+            @test expectation_value(rho_0_mps, (1, 2) => S_xx()) ≈ 0
+            @test expectation_value(rho_mps, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2
+            @test expectation_value(rho_mps, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2
+        end
 
-        # exact diagonalization
-        H_dense = convert(TensorMap, H)
-        Z_dense_1 = tr(exp(-beta * H_dense))^(1 / L)
-        Z_dense_2 = tr(exp(-2beta * H_dense))^(1 / L)
-
-        # taylor cluster
-        rho_taylor_1 = make_time_mpo(H, -im * beta, TaylorCluster(; N = 2))
-        Z_taylor_1 = tr(rho_taylor_1)^(1 / L)
-        @test Z_taylor_1 ≈ Z_dense_1 atol = 1.0e-2
-        Z_taylor_2 = real(dot(rho_taylor_1, rho_taylor_1))^(1 / L)
-        @test Z_taylor_2 ≈ Z_dense_2 atol = 1.0e-2
-
-        # MPO multiplication
-        rho_mps = convert(FiniteMPS, rho_taylor_1)
-        rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity))
-        Z_mpomul = tr(convert(FiniteMPO, rho_mps))^(1 / L)
-        @test Z_mpomul ≈ Z_dense_2 atol = 1.0e-2
-
-        # TDVP
-        rho_0 = MPSKit.infinite_temperature_density_matrix(H)
-        rho_0_mps = convert(FiniteMPS, rho_0)
-        rho_mps, = timestep(rho_0_mps, H, 0.0, im * beta, TDVP2(; trscheme))
-        Z_tdvp = real(dot(rho_mps, rho_mps))^(1 / L)
-        @test Z_tdvp ≈ Z_dense_2 atol = 1.0e-2
+        @testset "Infinite-size" begin
+            H = transverse_field_ising()
+            trscheme = truncdim(20)
+            verbosity = 1
+            beta = 0.1
+
+            # taylor cluster
+            alg_taylor = TaylorCluster(; N = 2)
+            rho_taylor = make_time_mpo(H, beta, alg_taylor; imaginary_evolution)
+            rho_taylor_2 = make_time_mpo(H, 2beta, alg_taylor; imaginary_evolution)
+            E_taylor = @constinferred expectation_value(rho_taylor, H)
+            E_taylor2 = @constinferred expectation_value(rho_taylor_2, H)
+
+            E_z_taylor = @constinferred expectation_value(rho_taylor, 1 => S_z())
+            @test E_z_taylor ≈ 0 atol = 1.0e-4 # no spontaneous symmetry breaking at finite T
+
+            # WII
+            rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution)
+            @test expectation_value(rho_wii, H) ≈ E_taylor atol = 1.0e-2
+            @test expectation_value(rho_wii, 1 => S_z()) ≈ E_z_taylor atol = 1.0e-2
+
+            # MPO multiplication
+            rho_mps = convert(InfiniteMPS, rho_taylor)
+            rho_mps2, = approximate(rho_mps, (rho_taylor, rho_mps), IDMRG(; verbosity))
+            E_mpomul = expectation_value(rho_mps2, H)
+            @test E_mpomul ≈ E_taylor2 atol = 1.0e-2
+
+            # TDVP
+            rho_0 = MPSKit.infinite_temperature_density_matrix(H)
+            rho_0_mps, = changebonds(convert(InfiniteMPS, rho_0), H, OptimalExpand(; trscheme = truncdim(20)))
+            rho_mps_tdvp, = timestep(rho_0_mps, H, 0.0, beta, TDVP(); imaginary_evolution)
+            E_tdvp = expectation_value(rho_mps_tdvp, H)
+            @test E_tdvp ≈ E_taylor atol = 1.0e-2
+
+            num_vals = 2
+            vals_taylor = @constinferred(transfer_spectrum(convert(InfiniteMPS, rho_taylor); num_vals))
+            vals_mps = @constinferred(transfer_spectrum(rho_mps; num_vals))
+            @test vals_taylor[1:num_vals] ≈ vals_mps[1:num_vals]
+        end
     end
 
     @testset "Sector conventions" begin
diff --git a/test/operators.jl b/test/operators.jl
index 3ea1ae003..61e26c67d 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -124,7 +124,7 @@ module TestOperators
     @testset "Finite MPOHamiltonian" begin
         L = 3
         T = ComplexF64
-        for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1))
+        for T in (Float64, ComplexF64), V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1))
             lattice = fill(V, L)
             O₁ = randn(T, V, V)
             O₁ += O₁'
@@ -136,6 +136,20 @@ module TestOperators
             H2 = FiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:(L - 1))
             H3 = FiniteMPOHamiltonian(lattice, 1 => O₁, (2, 3) => O₂, (1, 3) => O₂)
 
+            @test scalartype(H1) == scalartype(H2) == scalartype(H3) == T
+            if !(T <: Complex)
+                for H in (H1, H2, H3)
+                    Hc = @constinferred complex(H)
+                    @test scalartype(Hc) == complex(T)
+                    # cannot define `real(H)`, so only test elementwise
+                    for (Wc, W) in zip(parent(Hc), parent(H))
+                        Wr = @constinferred real(Wc)
+                        @test scalartype(Wr) == T
+                        @test Wr ≈ W atol = 1.0e-6
+                    end
+                end
+            end
+
             # check if constructor works by converting back to tensormap
             H1_tm = convert(TensorMap, H1)
             operators = vcat(fill(E, L - 1), O₁)