Crystals

Now consider a crystal. The Kohn-Sham hamiltonian has the same periodicity of the crystal, because both the external potential and the density $\rho$ have this periodicity. The fact that we are dealing with electrons that are not interacting amongst themselves, but only with a external effective field that has the same periodicity of the crystal, means that the one electron states have the form

$$\psi_{\bf k}({\bf r}) = e^{i {\bf k}\cdot {\bf r}} u_{\bf k} ({\bf r}),$$ (8.43)

where $u_{\bf k} ({\bf r})$ is a function that has the same periodicity of the lattice and ${\bf k}$ is a vector in the BZ. This is known as the Bloch theorem (see Appendix ). There is one such state for every vector in the BZ. In fact, for every vector ${\bf k}$ the hamiltonian has an infinite number of states, usually labelled with a letter $n$, representing different bands. The eigenvectors and eigenvalues of  are therefore usually labelled $\psi_{{\bf k},n}$ and $\epsilon_{{\bf k},n}$. To obtain the charge density one should in principle include all states in the BZ which are occupied by electrons. In an infinite crystal there would be $N/2$ occupied states for every k-point in the BZ, where $N$ is the number of electrons per primitive cell, and so the energy is:

$$E = 2\sum_{n=1}^{N/2} \frac{1}{\Omega}\int \epsilon_{{\bf k},n} d... ... \frac{\delta E_{xc}}{\delta \rho ({\bf r})}\rho ({\bf r})d^3 {\bf r} + E_{xc},$$ (8.44)

where $\Omega = \frac{(2\pi)^3}{V}$, and the charge density is obtained as:

$$\rho({\bf r}) = 2\sum_{n=1}^{N/2} \frac{1}{\Omega}\int \psi_{{\bf k},n}^*({\bf r})\psi_{{\bf k},n}({\bf r}) d^3{\bf k}.$$ (8.45)

In practice, the integrals over the BZ are approximated with finite sums over some appropriately chosen grid of k-points, for example a grid of uniformly distributed points ^8.10, or even using a single point, as originally proposed by Baldereschi ^8.11. The choice and the density of points in the BZ determines the accuracy of the calculations. However, in the case of metals the number of occupied bands in different parts of the BZ can differ, because electrons fill states up to the Fermi energy $\epsilon_F$ where there is a discontinuity (at $T=0$) in the occupation numbers. This discontinuity defines a Fermi surface, which can cause problems in the procedure to drive the KS equations to self-consistency. This is due to the fact that the Fermi surface moves as the self-consistency procedure is progressed from one step to the next, and so the number of occupied states at k-points that are close to the Fermi surface changes from one step to the next, which brings discontinuities in the charge density. A solution to this problem is to replace the discontinuity in the occupation numbers with a continuous variation, an electronic smearing, akin to introducing finite temperature in the system. The Hohenberg-Kohn formulation of DFT is only valid for the GS, however, it was later shown by Mermin ^8.12 that it can be extended to finite temperature, as long as one replaces the electronic energy $E$ with the free energy $F = E - T S$, where the electronic entropy is:

$$S = -k_{\rm B}2\sum_{n=1}^{\infty} \frac{1}{\Omega} \int [f_{{\bf k},n} \ln f_{{\bf k},n} + (1-f_{{\bf k},n})\ln (1-f_{{\bf k},n})]d^3{\bf k},$$ (8.46)

and the occupation numbers $f_{{\bf k},n}$ are defined as in . The charge density then becomes:

$$\rho({\bf r}) = 2\sum_{n=1}^{\infty} \frac{1}{\Omega}\int f_{{\bf k},n} \psi_{{\bf k},n}^*({\bf r})\psi_{{\bf k},n}({\bf r}) d^3{\bf k},$$ (8.47)

and the sum over the KS eigenvalues in  is replaced with:

$$2\sum_{n=1}^{\infty} \frac{1}{\Omega}\int f_{{\bf k},n} \epsilon_{{\bf k},n} d^3{\bf k}.$$ (8.48)

The removal of the discontinuity in the occupation numbers near the Fermi surface helps with convergence, but of course provides a result at finite temperature. If one is interested in a zero temperature calculation one needs to extrapolate to zero temperature, for example by performing calculations at several temperatures. To minimise the finite temperature bias and speed up the extrapolation procedure Methfessel an Paxton proposed to replace the standard Fermi smearing  with a modified one ^8.13.

Here comes an important point. Imagine that instead of using the primitive cell one decided to perform the calculations with a supercell, say a $2 \times 2 \times 2$ for example. We have of course complete freedom of doing that, as any supercell which is periodically repeated to tile the entire space can be used to describe the infinite crystal equally well. As we have seen in the in the previous chapter, if we double the lattice vectors then the reciprocal vectors are halved, which means that the BZ corresponding to this supercell is half the size in each direction. This means that in order to maintain the same density of points in the BZ we need to half the grid divisions in each direction, compared to the divisions used for the primitive cell. This principle is very important, especially for the calculations of phonon frequencies using supercells. We mentioned in the previous chapter that phonon frequencies at some special points in the BZ are exact. For example, for a $2 \times 2 \times 2$ supercell they are exact at zone boundary. A $4 \times 4 \times 4$ supercell has exact phonons at zone boundary, and in addition also to points in the middle of the BZ. This means that the zone boundary phonon frequencies should be the same as computed with either the $2 \times 2 \times 2$ or the $4 \times 4 \times 4$ supercells. This is only true, however, if the electronic structure calculations are performed with the same level of accuracy, which is only achieved if the same set of points in reciprocal space are used. This is realised if, for example, one centres both grids at $\Gamma$ and halves the number of divisions for the $4 \times 4 \times 4$ supercell compared to the $2 \times 2 \times 2$ supercell calculation.