Why do we specifically need an SCF-generated charge density for calculating the density of states (DOS) and band structure?

Question

Why do we specifically need an SCF-generated charge density for calculating the density of states (DOS) and band structure?

Answer 1

Background

At the heart of a conventional DFT simulation are the Kohn-Sham equations, $$ \left(-\frac{\hbar^2}{2m}\nabla^2 + \hat{V}\left[\rho\right]\right)\psi_{bk} = E_{bk}\psi_{bk} \tag{1} $$ where $\psi_{bk}$ is the $b$th Kohn-Sham single-particle state (band, for a material) at point $k$ in the Brillouin zone, with corresponding eigenenergy $E_{bk}$. $\hat{V}\left[\rho\right]$ is the potential operator, which is a functional of the Kohn-Sham density, $$ \rho\left(\mathbf{r}\right) = \sum_{bk^\prime}f_{bk^\prime}\left\vert\psi_{bk^\prime}\left(\mathbf{r}\right)\right\vert^2, \tag{2} $$ where $f_{bk^\prime}$ is the band occupancy at the reciprocal space point $k^\prime$, $\mathbf{r}$ is the position in 3D space. I have used a different k-point label, $k^\prime$, deliberately here, because the density requires us to integrate over the whole Brillouin zone, so the set of k-points in this expression must form an integration grid suitable for numerical integration (not necessarily uniform, though that is the common choice).

There is no requirement for the k-points labelled $k$ in Eq. $(1)$ to be the same set as the $k^\prime$ points used in Eq. $(2)$, but complications arise when the sets are different, and these complications are directly relevant to the question.

Self-Consistent Field calculation

Solving Eq. $(1)$ requires us to know $\hat{V}\left[\rho\right]$, and this in turn requires us to know $\rho(\mathbf{r})$. The problem is that $\rho(\mathbf{r})$ requires us to know $\psi_{bk}$ (from Eq. $(2)$), which is what we were trying to solve in the first place!

The most common solution to this conundrum is to use the self-consistent field (SCF) method, in which Eq. $(1)$ is solved for an initial guess "input" density $\rho^\mathrm{in}$, or potential $\hat{V}$, giving an estimate of the states $\psi_{bk}$. These states can then be used in Eq. $(2)$ to provide an "output" density $\rho^\mathrm{out}$ and, hence, a new potential. Eq. $(1)$ is now solved with this new potential, which gives new estimates for the states, which gives a new density... and so on.

If this iterative method converges, we eventually end up with the same "output" density as we put in, i.e. at some point $\rho^\mathrm{out}=\rho^\mathrm{in}=\rho$, which means that equations $(1)$ and $(2)$ are satisfied simultaneously, and the states $\psi_{bk}$, the density $\rho$ and the potential $\hat{V}\left[\rho\right]$ are all consistent -- we have reached self-consistency.

(NB in practice it is usually necessary to construct the new density (or potential) for the next iteration in a more sophisticated way, mixing together many previous inputs and outputs in a method known as density- or potential-mixing.)

Bandstructure calculations

When we wish to compute the bandstructure for a particular system, all we do is to solve Eq. $(1)$ at the appropriate k-points; for a bandstructure calculation, this would be a set of k-points along the system's high-symmetry directions.

In order to solve Eq. $(1)$ for our bandstructure $k$-points, we need to know the potential, and that means we need to know the density $\rho$. In principle, we could use the SCF method again and compute $\rho$ from Eq. $(2)$, except that the $k$-points we use for the bandstructure are extremely atypical, and do not cover large regions of the Brillouin zone, so they make an extremely poor choice for an integration grid.

The solution is to compute the density using a Brillouin zone integration grid ($k^\prime$ points) with the SCF method, and then use that density in Eq. $(1)$ and solve the equation for the states and eigenenergies at the bandstructure points $k$.

Density of States calculations

In contrast to bandstructure calculations, a density of states (DOS) calculation does require a set of Brillouin zone sampling points which form an integration grid. In fact, typical DOS calculations use a much finer integration grid than the SCF set, so wouldn't it actually be better to use this for the density as well?

The reason we use the SCF density here is much more pragmatic: a typical DOS Brillouin zone integration grid is much more densely sampled than we need, to represent the integrated density well. It is more efficient computationally to compute the density with an SCF optimisation on a coarser grid of $k^\prime$-points, and then "fix" the density and solve Eq. $(1)$ on the finer $k$-point grid for the DOS non-self-consistently, i.e. without constantly recomputing the density and potential using that finer grid.