Castep

Castep is a plane-wave, pseudopotential density-functional theory program for calculating the groundstate electronic charge density of a periodic system of electrons and nuclei.

The main computational effort is the search for the Kohn-Sham wavefunctions and groundstate electronic charge density for a fixed, 3D-periodic configuration of ions. This requires the solution of a series of one-particle Schrödinger equations:

$$\hat{H}\psi_{bk}(r) = E_{bk}\psi_{bk}(r)$$ (2.1)

where $\hat{H}$ is the Hamiltonian, the index $b$ labels different eigenstates, also known as bands, and the index $k$ labels different sampling points (`k-points') in the reciprocal-space Brillouin zone of the periodic simulation cell. $\{\psi_{bk}(r)\}$ are the Kohn-Sham wavefunctions, the eigenstates of these Schrödinger equations, and $\{E_{bk}\}$ are the corresponding band-energies.

The equations for different k-points are not completely independent, but interact via the electronic charge density $n(r)$

$$n(r) = \sum_{bk} f_{bk} \int \vert \psi_{bk}(r) \vert^2 d^3r$$ (2.2)

where $f_{bk}$ is the occupancy of band $b$ at k-point $k$. For insulating systems all bands below the Fermi energy are fully occupied ($f_{bk}=1$) and all bands above the Fermi energy are unoccupied ($f_{bk}=0$). For metals it is common to introduce a small thermal-like distribution function which smears the occupancies in the vicinity of the Fermi level.

For the 3D periodic systems Castep is designed to simulate, the density shares the same periodicity as the simulation system and so it is natural to express it in a Fourier basis. The basis functions are plane-waves whose wave-vectors are reciprocal lattice vectors of the simulation cell, the so-called G-vectors. The bands themselves need not be quite periodic, but they can always be written as the product of a phase factor $e^{i\vec{k}.\vec{r}}$ and a periodic function. It is these k-points that we sample, approximating the integral in equation 2.2 by a summation over discrete k-vectors drawn uniformly from the first Brillouin zone.

The Hamiltonian consists of some terms which are diagonal, or near-diagonal, in reciprocal space (e.g. the kinetic energy operator) and some terms which are diagonal in real-space (e.g. the ionic potential operator). This necessitates the use of Fourier transforms to transform the data from reciprocal- to real-space and back again.

Diagonalising the Hamiltonian directly is too expensive to be practical, and yields far more eigenstates than the lowest few we require; hence the search for the groundstate is an iterative procedure to compute just the lowest $N_b$ eigenstates, and in Castep this search can proceed in one of two main modes: density mixing (DM) or ensemble density functional theory (EDFT).

In the DM methodology the search proceeds by computing approximate eigenstates for a fixed trial charge density, computing a new density from these eigenstates, and then employing a density mixer to produce a new estimate for the true groundstate density.

Figure 2.1: Diagram showing the basic Castep self-consistent field (SCF) cycle for the density mixing (2.1(a)) and ensemble DFT (2.1(b)) algorithms

[Density mixing (DM) algorithm] [Ensemble density functional algorithm (EDFT)]

In the EDFT method the density used is always that obtained from the Kohn-Sham wavefunctions. A trial step is performed along the search direction and the density and energy recomputed. Castep then uses the data from this sample point, along with the data at the initial point, to fit a quadratic and estimate the so-called self-consistent minimum step, and hence find a better trial wavefunction. For metallic systems the band-occupancies are determined in a similar manner, updating the density and energy at each trial step of an occupancy cycle and always aiming to take steps that lower the total energy.

Because the EDFT method requires the construction of the density many more times than the DM method, it performs many more Fourier transforms of the wavefunction and is significantly slower per SCF cycle. However by fixing the Hamiltonian in the DM wavefunction updates, the DM algorithm has a tendency to overshoot the true eigenstate and it is only the density mixer itself that drives the algorithm to the correct groundstate. For this reason it is common for DM calculations to take many more SCF cycles to converge than corresponding EDFT calculations, though each DM cycle is considerably faster.