For a Hubbard hamiltonian, $$\hat{H} = \sum_{i,j,\sigma} t_{i j} \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} + \sum_i U_i \hat{n}_{i \uparrow} \hat{n}_{i \downarrow}, $$
the mean-field solution can be obtained by using:
$ \hat{n}_{i \uparrow} \hat{n}_{i \downarrow} \approx n_{i \uparrow} \hat{n}_{i \downarrow} + n_{i \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{i \downarrow}  $
and solving self-consistenly the resulting one-electron hamiltonian.
In general, for a generic many-body term $\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma},$ the approximation is:
$$\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma} \approx \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} + \langle   \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle  \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} - \langle   \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle \langle  \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle  \\
-\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{k \sigma^\prime} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{l \sigma} - \langle   \hat{c}^\dagger_{j \sigma^\prime}  \hat{c}_{l \sigma}\rangle  \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime} + \langle   \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle \langle  \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime}  \rangle   $$

(we assume the particle number is conserved, so that terms like $\langle \hat{c}^\dagger \hat{c}^\dagger \rangle $ are not important). 
Now, $S_z$ also being a good quantum number, if $\sigma \neq \sigma^\prime$ we will also ignore the terms on the second line, as done for the approximation above for Hubbard model. The Hubbard model on a lattice can be solved easily on mean-field thanks to the fact that on the self-consistency loop we only need the charges, and never non-diagonal charges (other elements of the density matrix) like $\langle   \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle$. 

 However think now of a slightly more complicated model on a lattice. It's like the hamiltonian $\hat{H}$ above but now we add many-body interactions between neighboring sites; something like:
$$ \dfrac{1}{2}\sum_{i,j, \sigma,\sigma^\prime} J_{ij} \hat{n}_{i \sigma} \hat{n}_{j \sigma^\prime}.$$
Now, for the terms with $\sigma \neq \sigma^\prime$ we have the simple approximation $\hat{n}_{i \uparrow} \hat{n}_{j \downarrow} \approx n_{i \uparrow} \hat{n}_{j \downarrow} + n_{j \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{j \downarrow}$. However, for terms like $\hat{n}_{i \sigma} \hat{n}_{j \sigma}$, when doing the mean-field approximation we will unavoidable find terms like those of the second line for the generic equation above; that is, the term:
$$-\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma} - \langle   \hat{c}^\dagger_{j \sigma}  \hat{c}_{i \sigma}\rangle  \hat{c}^\dagger_{i \sigma}\hat{c}_{j \sigma} + \langle   \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \langle  \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma}  \rangle,$$
and here is the problem. We need to include somehow in the self-consistency the non-diagonal charges $n_{ij\sigma}:=\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle. $ 

In order to do self-consistency on the charges we should express the $n_{ij \sigma}$ in terms of the $n_{i \sigma},$ but I don't think there is any analytical way of doing this even when all expected values are taken over Slater deterinants. Am I right? 
I know that I can express sums of non-diagonal charges, $\sum_j n_{i j \sigma}$ comfortably as sums of electron-hole interactions, $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$, which is convenient to write the Hartree-Fock energy of lattice models when the $J_{ij}$ parameters are constant. But in general, an isolated $n_{ij \sigma}$ seems untreatable. 
Two possible solutions come to mind:

1) Finding approximations for $n_{ij \sigma}$ in terms of $n_{k \sigma}$. How can this be done? Can the identity with $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$ (that is exact over Slater determinants) be used?

2) Doing self-consistency in terms of $n_{ij \sigma}$ instead of just the charges. I have tried this last option for things like PPP parametrization on hydrocarbons and I've encountered cases on which, even using a mixing algorithm on all the density matrix, the convergence is quite problematic. Has this be done before and how? (meaning, Lattice Hartree-Fock where self-consistency is done in terms of the full density matrix)