diff --git a/src/docstrings.jl b/src/docstrings.jl
index 9ee2e319..a542393a 100644
--- a/src/docstrings.jl
+++ b/src/docstrings.jl
@@ -2666,12 +2666,13 @@ where `U` is the onsite interaction.
 - `potential`: charge-charge potential to use for both Hartree and Fock. Can be a number or a function of position. Default: `1`
 - `hartree`: charge-charge potential `v_H` for the Hartree mean field. Can be a number or a function of position. Overrides `potential`. Default: `potential`
 - `fock`: charge-charge potential `v_F` for the Fock mean field. Can be a number, a function of position or `nothing`. In the latter case all Fock terms (even onsite) will be dropped. Default: `hartree`
-- `onsite`: charge-charge onsite potential. Overrides both Hartree and Fock potentials for onsite interactions. Default: `hartree(0)`
+- `onsite`: charge-charge onsite potential. Overrides both Hartree and Fock potentials for onsite interactions, unless `fock = nothing` (it then overrides only Hartree, see above). Default: `hartree(0)`
 - `charge`: a number (in single-orbital systems) or a matrix (in multi-orbital systems) representing the charge operator on each site. Default: `I`
 - `nambu::Bool`: specifies whether the model is defined in Nambu space. In such case, `charge` should also be in Nambu space, typically `SA[1 0; 0 -1]` or similar. Default: `false`
 - `namburotation::Bool`: if `nambu == true` and spinful systems, specifies whether the spinor basis is `[c↑, c↓, c↓⁺, -c↑⁺]` (`namburotation = true`) or `[c↑, c↓, c↑⁺, c↓⁺]` (`namburotation = false`). Default: `false`
 - `selector::NamedTuple`: a collection of `hopselector` directives that defines the pairs of sites (`pair_selection` above) that interact through the charge-charge potential. Default: `(; range = 0)` (i.e. only onsite)
-- `selector_hartree::NamedTuple`: same as `selector`, but restricted to the Hartree interaction. Useful to do efficient Hartree-only simulations, by setting `fock = 0` and leaving `selector` at its `(; range =0)` default. Then, having a large range in `selector_hartree` will be cheap. Default: `selector`.
+- `selector_hartree::NamedTuple`: same as `selector`, but restricted to the Hartree interactions. Useful e.g. to do efficient Hartree-only simulations: a long-range `selector_hartree` is cheap in periodic systems, as lattice periodicity is exploited to make the computational cost independent of the range. Default: `selector`.
+- `selector_fock::NamedTuple`: same as `selector_hartree`, but for Fock interactions. Note that `selector_fock` increases the number of entries that need to be computed for the density matrix, so increasing the Fock range is typically not cheap. Default: `selector`.
 
 Any additional keywords `kw` are passed to the `densitymatrix` function used to compute the
 mean field, see above
diff --git a/src/meanfield.jl b/src/meanfield.jl
index de666312..808af7e4 100644
--- a/src/meanfield.jl
+++ b/src/meanfield.jl
@@ -26,7 +26,7 @@ struct ZeroField end
 function meanfield(g::GreenFunction{T,E}, args...;
     potential = Returns(1), hartree = potential, fock = hartree,
     onsite = missing, charge = I, nambu::Bool = false, namburotation = missing,
-    selector::NamedTuple = (; range = 0), selector_hartree = selector, kw...) where {T,E}
+    selector::NamedTuple = (; range = 0), selector_fock = selector, selector_hartree = selector, kw...) where {T,E}
 
     Vh = sanitize_potential(hartree)
     Vf = sanitize_potential(fock)
@@ -40,13 +40,13 @@ 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`")
 
-    gs = g[sitepairs(; selector..., includeonsite = true)]
+    gs = g[sitepairs(; selector_fock..., 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.
-    hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : T(Vf(dr)); selector..., includeonsite = true)
-    hHartree = (Uf == U && Vh == Vf && selector == selector_hartree) ? hFock :
+    hFock = lat |> hopping((r, dr) -> iszero(dr) ? Uf : T(Vf(dr)); selector_fock..., includeonsite = true)
+    hHartree = (Uf == U && Vh == Vf && selector_fock == selector_hartree) ? hFock :
         lat |> hopping((r, dr) -> iszero(dr) ? U : T(Vh(dr)); selector_hartree..., includeonsite = true)
     # this drops zeros
     potHartree = T.(sum(unflat, harmonics(hHartree)))
@@ -66,7 +66,7 @@ function meanfield(g::GreenFunction{T,E}, args...;
 
     encoder, decoder = nambu ? NambuEncoderDecoder(is_nambu_rotated´) : (identity, identity)
     S = typeof(encoder(zero(Q)))
-    output = call!_output(g[sitepairs(; selector..., includeonsite = true, kernel = Q)])
+    output = call!_output(g[sitepairs(; selector_fock..., includeonsite = true, kernel = Q)])
     sparse_enc = similar(output, S)
     output = CompressedOrbitalMatrix(sparse_enc; encoder, decoder, hermitian = true)