Integrating over all orbitals

Computing expectation values of observables O over all Kohn-Sham orbitals is an important concept relevant to most ab initio calculations. Typically, we can evaluate these expectation values as integral over all orbitals

$\langle O \rangle = \sum_n \frac{1}{\Omega_{\mathrm{BZ}}} \int_{\Omega_{\mathrm{BZ}}} O_{n\bold{k}} \, \Theta(\epsilon_{\mathrm{F}} - \epsilon_{n\bold{k}}) \, d^{3}k,$

where Ω_BZ is the volume of the Brillouin zone, O_nk is the expectation value of the observable with a single Kohn-Sham orbital, and Θ is the Heaviside step function and limits the integral to orbitals with eigenvalues ϵ_nk below the Fermi energy ϵ_F.

When evaluating the integral above numerically, we need to address two concerns: (i) We need to discretize the k-point integral since we do not know the analytic expression of the observable. (ii) A sharp function cutoff like the Heaviside function is numerically unstable so we need robust smearing methods.

Integrating over k points

Discretizing the k-point integral involves replacing the continuous integral over the Brillouin zone by a k-point mesh.

$\frac{1}{\Omega_{\mathrm{BZ}}} \int_{\Omega_{\mathrm{BZ}}} \to \sum_{\bold{k}} w_{\bold{k}}.$

A k-point mesh consists of k-point coordinates and associated weights w_k. In general, any mesh could be chosen but for periodic boundaries the optimal one are equidistant grids. In VASP, we select this sampling by a KPOINTS file or the KSPACING tag.

We can improve the integrals further because the crystal exhibits certain symmetries. Often we can deduce the value O_nk from a different symmetry-equivalent k point. VASP automatically analyzes the symmetry of the crystal and reduces the k point mesh to the irreducible Brillouin zone. The weights of each irreducible k point measure how many equivalent k points exist in the reducible Brillouin zone.

Integrating near the Fermi energy

Occupation of different smearing techniques near the Fermi energy ϵ_F. The energy is measured in units of the smearing σ (SIGMA). The Methfessel-Paxton method (cyan) is closer to the step function than the Gaussian distribution (blue) but also has nonmonotonous features. The smearing Fermi-Dirac distribution (purple) corresponds to a temperature but is much broader at the same value of σ.

We want to replace the Heaviside step function by a smooth equivalent to make the integral numerically stable. Otherwise small changes in the energies would toggle the inclusion of a sample O_nk close to the Fermi energy.

$\sum_{n\bold{k}} \Theta(\epsilon_{\mathrm{F}} - \epsilon_{n\bold{k}}) \ldots \to \sum_{n\bold{k}} f_\sigma(\epsilon_{n\bold{k}}) \ldots$

Here, the occupations f_σ(ϵ_nk) approach 1 for energies far below the Fermi energy ϵ_F and 0 for energies far above it. The parameter σ determines how wide the broadened step function is.

In the figure, we illustrate the different smearing methods implemented in VASP. The Fermi-Dirac smearing (ISMEAR = -1) uses σ as the temperature in the Fermi-Dirac distribution

$f_\sigma(\epsilon) = \Bigl[\exp(\frac{\epsilon-\epsilon_{\mathrm{F}}}{\sigma})+1\Bigr]^{-1}~.$

One can also use the complementary error function (ISMEAR = 0) which results from integrating a Gaussian distribution as the occupation function. This leads to a narrower edge than the Fermi-Dirac distribution for the same σ and a faster approach to the asymptotic behavior

$f_\sigma(\epsilon) = \frac{1}{2} \text{erfc}\Bigl[\frac{\epsilon-\epsilon_{\mathrm{F}}}{\sigma}\Bigr]~.$

Methfessel and Paxton developed higher order approximations to the step function (ISMEAR > 0).

$f_\sigma(\epsilon) = \frac{1}{2} \text{erfc}\Bigl[\frac{\epsilon-\epsilon_{\mathrm{F}}}{\sigma}\Bigr] + \exp\Bigl[-\bigl(\frac{\epsilon-\epsilon_{\mathrm{F}}}{\sigma}\bigr)^2\Bigr] \sum_{i=1}^{n} \frac{(-1)^i}{4^i i! \sqrt\pi} H_{2i-1}\Bigl[\frac{\epsilon-\epsilon_{\mathrm{F}}}{\sigma}\Bigr]~.$

Here, n is the order of the expansion and H_j is the j-th Hermite polynomial. The first order Methfessel-Paxton smearing is shown in the figure. This method leads to an even narrower distribution but introduces a nonmonotonous behavior that can lead to problems in semiconductors and insulators.

A consequence of these broadening techniques is that the total energy is no longer variational (or minimal). It is necessary to replace the total energy by some generalized free energy

$F = E - \sum_{n\bold{k}} w_{\bold{k}} \sigma S[f_\sigma(\epsilon_{n\bold{k}})].$

For the Fermi-Dirac statistics, we might interpret this as the free energy of the electrons at some finite temperature σ = k_BT. There is no straightforward interpretation of the free energy in the case of Gaussian or Methfessel-Paxton smearing. Despite this, it is possible to obtain an accurate extrapolation for σ → 0 from results at finite σ using the formula

$E_0 = E(\sigma \to 0) = \frac{1}{2} (F + E)~.$

E₀ is a meaningful physical quantity for the ground state energy of the system. Importantly, the calculated forces are the derivatives of the free energy F and not of E₀. Nonetheless, the difference of the forces is generally small and acceptable it a suitable σ is used.

When we consider E₀ as our target property, the smearing methods serve as a mathematical tool to obtain faster convergence with respect to the number of k-points. Generally, the Gaussian broadening requires more careful tuning of the width σ compared to the Methfessel-Paxton method. If σ is too large, the energy E(σ → 0) will converge to the wrong value even for an infinite k-point mesh. If σ is too small, we require a much denser k-point mesh and a significantly larger computational cost. With the Methfessel-Paxton method the sharper edge usually averts the necessity of tuning σ. However, since the occupation function is nonmonotonous, it is not suitable to describe systems with a bandgap.

Tetrahedron method

Four k points forming a single tetrahedron with ordered eigenvalues ϵ₁, ϵ₂, ϵ₃, and ϵ₄. The lines show the occupation of the orbitals with changing Fermi energy ϵ_F where the darkest line corresponds to the lowest and the brightest line to the highest eigenvalue. For the Gaussian smearing (blue) every k point is considered individually leading to half filling when the Fermi energy reaches the eigenvalue. The tetrahedron method (red) does not fill any orbital if the Fermi energy is outside of the bounds of the tetrahedron.

The tetrahedron method is an alternative approach to address the sharp edge of the Heaviside step function. Instead of considering each k point individually, we triangulate the k-point mesh, i.e., we split it into as many tetrahedra as necessary to cover the whole Brillouin zone. We use a linear interpolation of the band energies ϵ_nk within each tetrahedron the band energies. Blöchl derived correction terms to cancel the linearization errors of the tetrahedron method.

With this interpolation, we can solve integral over the Brillouin zone analytically considering each tetrahedron individually. It does not require a choice of a width σ like the broadening methods. The figure illustrates the difference between broadening and interpolation method. A broadening method like the Gaussian smearing considers every k point individually. Therefore, it will start filling the orbital as soon as the Fermi energy reaches the width σ of the broadening. In the tetrahedron method, the k points of the tetrahedron only get occupied once the Fermi energy exceeds the lowest eigenvalue. Similarly, it is completely filled once the Fermi energy exceeds the maximum value. As a consequence, the occupations of the different k points are much closer to each other.

Overall, the tetrahedron method yields very accurate occupations with minimal user input. It is very well suited to obtain accurate integral (e.g. total energy) and band onsets (e.g. density of state). The main drawback is that the Blöchel's correction of the linearization errors is not variational with respect to the partial occupancies. Therefore the calculated forces might be wrong by a few percent. If accurate forces are required we recommend a finite temperature method.

Determining the Fermi energy

One important example of integrals over all orbitals is the calculation of the Fermi energy. In this case, the observable is identity and the sum of all occupations should be equal to the number of electrons N_e

$\sum_{n\bold{k}} f_\sigma(\epsilon_{n\bold{k}}) = N_{\mathrm{e}}~.$

This leads to a straightforward interval-bisection algorithm to compute the Fermi energy. We guess bounds for the Fermi energy and then compute the sum of all occupations in the middle of the interval. If this results in a number larger than the number of electrons, we replace the upper bound otherwise we replace the lower one.

Note that this algorithm is not deterministic for systems with a bandgap, where any Fermi energy in the gap would be fine. VASP achieves more consistent results for these system by a good initial guess for the Fermi energy (see EFERMI). This potentially leads to an early exit of the bisection algorithm so that only for metals many iterations of bisection are considered.

Related tags and sections

KPOINTS, ISMEAR, SIGMA, Smearing technique

---