Meaning of orbital occupations in the Hubbard Model I was reading the paper 'A linear response approach to the calculation of the effective interaction parametersin the LDA+U method' by Cococcioni et. al (1)
Quoted verbatim from section 3.A), 

In order to fully define how the approach works the first thing to do is to select the degrees of freedom on which 'HubbardU' will operate and define the corresponding occupation matrix, n I mm' . Although it is usually straight-forward to identify in a given system the atomic levels to be treated in a special way (the 'd' electrons in transition metals and the 'f' ones in the rare earths and actinides series) there is no unique or rigorous way to define occupation of localized atomic levels in a multi-atom system. Equally legitimate choices are projections on normalized atomic orbitals, projections on wannier functions etc.

I do not understand the physical interpretation of the occupation matrix projected on normalized atomic orbitals. Or wannier functions too, for that matter. What does this represent? 
PS: Sorry if I'm missing something very straightforward here.
 A: First of all, consider the concept of occupation numbers for an isolated atom. In this case, atomic occupations numbers are well-defined because we have a natural basis for constructing many-body states (Slater determinants), consisting of bound states of an electron moving in the effective one-electron potential of an atom. These bound states, $\phi_{nlm}(\mathbf{r})$, are localized because their wavefunctions decay exponentially, which implies that we can unambiguously associate them with a given atom. Also, thanks to rotation symmetry of an isolated atom, the states have well-defined orbital quantum numbers $(l,m)$ and are therefore referred to as orbital states.
In a crystalline solid, we encounter two issues:


*

*the atoms are sitting very close to each other, which means that, save for the core states (atomic-like states with rather high binding energies), practically all other states represent non-localized Bloch waves, $\Psi_{\nu\mathbf{k}}(\mathbf{r})$, characterized by band indices $\nu$ and by points $\mathbf{k}$ of the Brillouin zone. Such a wavefunction is usually very extended, having a significant overlap with many different atoms, which makes it meaningless to assign such a state to one particular atom;

*rotation symmetry is broken inside a solid, hence no well-defined quantum numbers $(l, m)$.
Now, we want to set up the Hubbard model on top of our Kohn-Sham DFT (Bloch) states. The very reason we want to do this is an observation that LDA (GGA, meta-GGA, etc.) functionals fail to describe certain classes of solid systems because of the presence of partially occupied localized states that resemble those of an isolated atom. What this means is that some of the Bloch waves can be represented by a linear combination of strongly localized atomic states, only slightly overlapping with each other. For instance, we look at the band structure of NiO and see some suspiciously flat bands around the Fermi level. It turns out that the states corresponding to these bands can be rather accurately represented by linear combinations of $3d$-states of an isolated Ni atom. In this case, one can say that the tight-binding approximation applies well to these states. And this is important because the Hubbard model is conceptually formulated in the tight-binding picture.
So, the conclusion is that if we want to use the Hubbard model (the basis of $+U$ methods), we first need to construct such localized states, $\phi(\mathrm{r})$, that they satisfy several criteria:


*

*Linear combinations of these states must produce our Bloch states exactly. Otherwise stated, the desired states must be formed by linear combinations of Bloch states,
$$
|\phi_{\lambda}\rangle = \sum_{\nu \mathbf{k}} c_{\nu\mathbf{k},\lambda} |\Psi_{\nu \mathbf{k}}\rangle,
$$
where $\lambda$ is a set of indices (usually, $\mathrm{R}$ -- atomic position, $l$, and $m$) characterizing the state.

*The states must be normalized and strongly localized around atoms (e.g., Ni atoms in NiO). That is, the overlap of $\phi_{\mathbf{R},lm}(\mathbf{r})$ and $\phi_{\mathbf{R'},l'm'}(\mathbf{r})$ for any pair of nearest neighbors at $\mathbf{R}$ and $\mathbf{R}'$ must be much smaller than 1.

*The states must have very strong orbital character, implying that
$$
\int \phi_{l',m'}(\mathrm{r}) Y_{l,m}(\hat{\mathrm{r}}) d \hat{\mathrm{r}} \approx 
\phi_{l,m}(r) \delta_{l',l} \delta_{m',m},
$$
where the integration is over angles and the origin is shifted to a respective atomic site. After all, the interaction term ($U$-term) in the Hubbard model usually represents screened Coulomb interaction between states of a spherically symmetric atom.
Once such a set of states is constructed, the local occupation number (density) matrix (per spin channel) is given by (omitting position indices $\mathbf{R}$)
$$
n_{lm,l'm'} = \langle \phi_{lm} | \phi_{l'm'} \rangle = 
   \sum_{\nu\mathbf{k}} f_{\nu\mathbf{k}} c^*_{lm,\nu\mathbf{k}} c_{\nu\mathbf{k},l'm'},
$$
where $f_{\nu\mathbf{k}}$ are Fermi weights cutting out only occupied Bloch states. In particular, $n_{lm,lm}$ yields the occupation number of a state $|\phi_{lm}\rangle$.
The procedure outlined above is generally called projection. Now, what the author of the cited paper refers to are various methods to perform projection.
For example, one way of doing this is to search for such coefficients $c_{\nu\mathbf{k},lm}$ that $\phi_{lm}(\mathbf{r})$ have the lowest possible extent defined by $\int r^2 |\phi_{lm}(\mathbf{r})|^2 d^3 r$. This way, one arrives at maximally localized Wannier functions (MLWF).
Another way of constructing localized states for some ion in a solid is to start from a basis of functions, $\chi_{lm}$, having proper atomic-like character. For instance, one could simply use bound states of a corresponding isolated atom. To construct the desired states, we perform projection by forming linear combinations of Bloch states, with the coefficients defined by the overlap of Bloch states with $\chi_{lm}$:
$$
|\phi_{lm}\rangle = \sum_{\nu \mathbf{k}} |\nu, \mathbf{k}\rangle \langle \nu, \mathbf{k}| \chi_{lm} \rangle.
$$
This method is called projected localized orbitals.
TL;DR: "Occupation matrix projected on [some set of localized functions]" means that in order to use the Hubbard model on top of DFT we transformed our Kohn-Sham states to form a set of localized wavefunctions resembling atomic bound states with well-defined orbital character $(l,m)$. Then, it is a matter of a straightforward calculation to derive the occupation of these localized states from the occupation of Bloch (Kohn-Sham) bands.
