diff --git a/src/docstrings.jl b/src/docstrings.jl
index a5547d8c..c8d5d115 100644
--- a/src/docstrings.jl
+++ b/src/docstrings.jl
@@ -2086,7 +2086,7 @@ The `quadgk_opts` are extra keyword arguments (other than `atol`) to pass on to
 
 Currently, the following GreenSolvers implement dedicated densitymatrix algorithms:
 
-- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries are not currently supported. No `opts`.
+- `GS.Schur`: based on numerical integration over Bloch phase. Boundaries and non-empty contacts are not currently supported. Assumes Hermitian Hamiltonian. No `opts`.
 - `GS.Spectrum`: based on summation occupation-weigthed eigenvectors. No `opts`.
 - `GS.KPM`: based on the Chebyshev expansion of the Fermi function. Currently only works for zero temperature and only supports `nothing` contacts (see `attach`). No `opts`.
 
diff --git a/src/meanfield.jl b/src/meanfield.jl
index b5c24ab2..3cab3f15 100644
--- a/src/meanfield.jl
+++ b/src/meanfield.jl
@@ -7,17 +7,16 @@
 #   we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry
 #region
 
-struct MeanField{B,T,C<:CompressedOrbitalMatrix,S<:SparseMatrixCSC,H<:DensityMatrix,F<:DensityMatrix}
+struct MeanField{B,T,C<:CompressedOrbitalMatrix,S<:SparseMatrixCSC,F<:DensityMatrix}
     output::C
     potHartree::S
     potFock::S
-    rhoHartree::H
-    rhoFock::F
+    rho::F
     charge::B
     nambu::Bool
     is_nambu_rotated::Bool
     rowcol_ranges::NTuple{2,Vector{UnitRange{Int}}}
-    onsite_tmp::Vector{Complex{T}}
+    onsite_tmp::Vector{T}
 end
 
 struct ZeroField end
@@ -41,10 +40,8 @@ function meanfield(g::GreenFunction{T,E}, args...;
     isempty(boundaries(g)) || argerror("meanfield does not currently support systems with boundaries")
     isfinite(U) || argerror("Onsite potential must be finite, consider setting `onsite`")
 
-    gsHartree = g[diagonal(; cells = 0, kernel = Q)]
-    rhoHartree = densitymatrix(gsHartree, args...; kw...)
-    gsFock = g[sitepairs(; selector..., includeonsite = true)]
-    rhoFock = densitymatrix(gsFock, args...; kw...)
+    gs = g[sitepairs(; selector..., includeonsite = true)]
+    rho = densitymatrix(gs, args...; kw...)
 
     lat = lattice(hamiltonian(g))
     # The sparse structure of hFock will be inherited by the evaluated mean field. Need onsite.
@@ -52,19 +49,19 @@ function meanfield(g::GreenFunction{T,E}, args...;
     hHartree = (Uf == U && Vh === Vf) ? hFock :
         lat |> hopping((r, dr) -> iszero(dr) ? U : T(Vh(dr)); selector..., includeonsite = true)
 
-    potHartree = sum(unflat, harmonics(hHartree))
+    potHartree = T.(sum(unflat, harmonics(hHartree)))
 
-    oaxes = orbaxes(call!_output(gsFock))
+    oaxes = orbaxes(call!_output(gs))
     rowcol_ranges = collect.(orbranges.(oaxes))
-    onsite_tmp = similar(diag(parent(call!_output(gsHartree))))
+    onsite_tmp = Vector{T}(undef, length(last(rowcol_ranges)))
 
-    # build potFock with identical axes as the output of rhoFock
+    # build potFock with identical axes as the output of rho
     cells_fock = cells(first(oaxes))
     hFock_slice = hFock[(; cells = cells_fock), (; cells = 0)]
 
     # this is important for the fast orbrange-based implementation of MeanField evaluation
-    check_cell_order(hFock_slice, rhoFock)
-    potFock = parent(hFock_slice)
+    check_cell_order(hFock_slice, rho)
+    potFock = T.(parent(hFock_slice))
 
     encoder, decoder = nambu ? NambuEncoderDecoder(is_nambu_rotated´) : (identity, identity)
     S = typeof(encoder(zero(Q)))
@@ -72,7 +69,7 @@ function meanfield(g::GreenFunction{T,E}, args...;
     sparse_enc = similar(output, S)
     output = CompressedOrbitalMatrix(sparse_enc; encoder, decoder, hermitian = true)
 
-    return MeanField(output, potHartree, potFock, rhoHartree, rhoFock, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp)
+    return MeanField(output, potHartree, potFock, rho, Q, nambu, is_nambu_rotated´, rowcol_ranges, onsite_tmp)
 end
 
 sanitize_potential(x::Real) = Returns(x)
@@ -91,6 +88,7 @@ sanitize_nambu_rotated(is_nambu_rotated::Bool) = is_nambu_rotated
 
 function sanitize_charge(charge, B, nambu, is_rotated)
     Q = sanitize_charge(charge, B)
+    ishermitian(Q) || argerror("Charge $Q should be Hermitian")
     nambu && check_nambu(Q, is_rotated)
     return Q
 end
@@ -137,9 +135,9 @@ function NambuEncoderDecoder(is_nambu_rotated)
     return encoder, decoder
 end
 
-function check_cell_order(hFock_slice, rhoFock)
+function check_cell_order(hFock_slice, rho)
     opot = first(orbaxes(hFock_slice))
-    orho = first(orbaxes(call!_output(rhoFock.gs)))
+    orho = first(orbaxes(call!_output(rho.gs)))
     cells(opot) == cells(orho) || internalerror("meanfield: Cell order mismatch between potential and density matrix")
     return nothing
 end
@@ -164,15 +162,15 @@ function call!(m::MeanField{B}, args...; chopsmall = true, params...) where {B}
     Q, hartree_pot, fock_pot = m.charge, m.onsite_tmp, m.potFock
     rowrngs, colrngs = m.rowcol_ranges
     check_zero_mu(m, args...)
-    trρQ = maybe_nambufy_traces!(diag(parent(m.rhoHartree(args...; params...))), m)
+    ρ = m.rho(args...; params...)
+    trρQ = maybe_nambufy_traces!(diag_real_tr_rho_Q(ρ, Q), m)
     mul!(hartree_pot, m.potHartree, trρQ)
-    ρFock = m.rhoFock(args...; params...)
     meanfield = m.output
     meanfield_parent = parent(meanfield)
     fill!(nonzeros(meanfield_parent), zero(eltype(meanfield_parent)))
     if chopsmall
         hartree_pot .= Quantica.chopsmall.(hartree_pot)  # this is a Vector
-        nzs = nonzeros(parent(ρFock))
+        nzs = nonzeros(parent(ρ))
         nzs .= Quantica.chopsmall.(nzs)
     end
     rows, cols, nzs = rowvals(fock_pot), axes(fock_pot, 2), nonzeros(fock_pot)
@@ -183,7 +181,7 @@ function call!(m::MeanField{B}, args...; chopsmall = true, params...) where {B}
             row = rows[ptr]
             row < col && ishermitian(meanfield) && continue   # skip upper triangle
             vij = nzs[ptr]
-            ρij = view(ρFock, rowrngs[row], colrngs[col])
+            ρij = view(ρ, rowrngs[row], colrngs[col])
             vQρijQ = vij * Q * sanitize_block(B, ρij) * Q
             if row == col
                 meanfield_parent[row, col] = encoder(meanfield)(viiQ - vQρijQ)
@@ -209,6 +207,9 @@ function maybe_nambufy_traces!(traces, m::MeanField{B}) where {B}
     return traces
 end
 
+diag_real_tr_rho_Q(ρ, Q) =
+    [real(unsafe_trace_prod(Q, view(ρ, rng, rng))) for rng in siteindexdict(orbaxes(ρ, 2))]
+
 hole_id(::Type{<:SMatrix{2,2}}) = SA[0 0; 0 1]
 hole_id(::Type{<:SMatrix{4,4}}) = SA[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]
 hole_id(S) = nambu_dim_error(S)
diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl
index 0af7eaaf..de3174c1 100644
--- a/src/solvers/green/schur.jl
+++ b/src/solvers/green/schur.jl
@@ -689,6 +689,7 @@ end
 
 # use computed Fermi points to integrate spectral function by segments
 # returns an AbstractMatrix
+# we don't use Integrator, because that is meant for integrals over energy, not momentum
 function (s::DensityMatrixSchurSolver)(µ, kBT; params...)
     g = parent(s.gs)
     λs = propagating_eigvals(g, µ + 0im, 1e-2; params...)
@@ -709,7 +710,7 @@ function fermi_h!(result, s, ϕ, µ, β = 0; params...)
     bs = blockstructure(h)
     # Similar to spectrum(h, ϕ; params...), but less work (no sort! or sanitization)
     copy!(s.hmat, call!(h, ϕ; params...))  # sparse to dense
-    ϵs, psis = eigen!(s.hmat)
+    ϵs, psis = eigen!(Hermitian(s.hmat))
     # special-casing β = Inf with views turns out to be slower
     fs = (@. ϵs = fermi(ϵs - µ, β))
     fpsis = (s.psis .= psis .* transpose(fs))
diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl
index cdc1c599..2e61df47 100644
--- a/test/test_greenfunction.jl
+++ b/test/test_greenfunction.jl
@@ -566,10 +566,12 @@ end
 @testset "meanfield" begin
     oh = LP.honeycomb() |> hamiltonian(hopping(SA[1 im; -im -1]) - onsite(SA[1 0; 0 -1], sublats = :A), orbitals = 2) |> supercell((1,-1))
     g = oh |> greenfunction
+    Q = SA[0 1; im 0]  # nonhermitian charge not allowed
+    @test_throws ArgumentError meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q)
     Q = SA[0 -im; im 0]
     m = meanfield(g; selector = (; range = 1), potential = 2, fock = 1.5, charge = Q)
     Φ = m(0.2, 0.3)
-    ρ12 = m.rhoFock(0.2, 0.3)[sites(1), sites(2)]
+    ρ12 = m.rho(0.2, 0.3)[sites(1), sites(2)]
     @test Φ[sites(1), sites(2)] ≈ -1.5 * Q * ρ12 * Q
 
     # nambu
@@ -579,7 +581,7 @@ end
     m = meanfield(g; selector = (; range = 1), nambu = true, hartree = r -> 1/(1+norm(r)), fock = 1.5, charge = Q)
     @test_throws ArgumentError m(0.2, 0.3)  # µ cannot be nonzero
     Φ = m(0, 0.3)
-    ρ11 = m.rhoFock(0, 0.3)[sites(1), sites(1)] - SA[0 0; 0 1]
+    ρ11 = m.rho(0, 0.3)[sites(1), sites(1)] - SA[0 0; 0 1]
     fock = -1.5 * Q * ρ11 * Q
     hartree = 0.5*Q*(tr(Q*ρ11)*(1+1/2+1/2))
     @test Φ[sites(1), sites(1)] ≈ hartree + fock