Background Theory: Fourier Transforms and Plane-Wave Expansions

Introduction

Density functional theory (DFT) calculations rely on solving the Kohn-Sham equations self-consistently. The wavefunction of the system can be represented in multiple basis sets. For material systems with periodic boundary conditions, a plane-wave basis set is frequently chosen, due to its periodic nature and numerical efficiency as discussed below. This strategy is adopted by many popular DFT packages in solid-state physics, including e.g. ABINIT, CASTEP, Quantum ESPRESSO, and VASP.

Plane-wave expansion

The general form of a plane wave is:

P(\mathbf{r}) = P_0 e^{i\mathbf{q}\cdot \mathbf{r}}.

The name “plane wave” originates from the fact that the wavefront (where the wave assumes the same value) is a plane perpendicular to the propagation direction $\mathbf{q}$ .

In solid-state physics, where one typically studies periodic crystalline solids, Bloch’s theorem tells us that any wavefunction can be written as $\psi_{n\mathbf{k}}(r) = u_{n\mathbf{k}}(r) e^{i\mathbf{k}\cdot \mathbf{r}}$ , where $u_{n\mathbf{k}}(r)$ is a function with the same periodicity of the system.

Since this function is periodic, we can then use a plane-wave basis set to represent $u_{n\mathbf{k}}(r)$ :

u_{n\mathbf{k}}(r) = \sum_{\mathbf{G}} c_i(\mathbf G) e^{i\mathbf{G\cdot r}}

where the sum is over all reciprocal lattice vectors $\mathbf G$ . Similarly, the charge density (that also has the periodicity of the system) is represented in a plane-wave basis set.

The main advantage of this basis set is the possibility to use fast Fourier transform (FFT) algorithms to convert between the reciprocal-space representation (the coefficient of the plane waves) and the real-space representation of the wavefunctions on a grid. This is appealing because the operation of the potential on $\psi$ is very simple in real space (since for local potentials it is simply an element-wise multiplication), while it is very efficient to compute the operation of the kinetic operator in reciprocal space, where the second derivative also becomes simply an element-wise multiplication. One can then pass between the two representations using FFTs, when needed.

Since the number of plane waves (or, equivalently, of $\mathbf G$ vectors) determines the size of the numerical problem to solve, it is important to find a balance between using enough plane waves to accurately represent the plane waves, but at the same time not too many, to avoid to make an overly expensive calculation. This notebook helps in visualizing the effect of the basis set size on the functions that we want to represent.

Fourier series and Fast fourier transform

In real space, a Fourier series is a linear combination of cosine and sine functions used to represent a smooth periodic function defined in a certain range. For a one dimensional function $f(x)$ defined on $[-\pi, \pi]$ , the Fourier series is given by:

f(x) = \frac {A_0} {2} + \sum_{k=1}^{\infty}(A_k \cos kx + B_k \sin kx)

where one can prove that the coefficients are given by:

A_k = \frac {1}{\pi} \int_{-\pi}^{\pi} f(x) \cos(kx)dx

and

B_k = \frac {1}{\pi} \int_{-\pi}^{\pi} f(x) \sin(kx)dx.

Note that you can interpret the Fourier coefficients as the inner product of $f(x)$ and a (normalized) cosine or sine function. Geometrically, this can be interpreted as the projection of the function onto a unit basis vector of the Hilbert space spanned by the cosine and sine functions. Thus, the Fourier series is essentially the representation of a function using an infinite basis set composed of trigonometric functions.

In complex space, a Fourier series shares the same form as a plane-wave expansion:

f(x)=\sum_{-\infty}^{\infty}C_k e^{ikx} = \sum_{-\infty}^{\infty}(\alpha_k + i\beta_k)(\cos kx + i\sin kx).

Note that, for both the real and complex case, the basis set is composed of trigonometric functions with increasing frequencies, all multiples of the base frequency, determined by the periodicity in real space.

If all frequencies are considered, under appropriate conditions the representation is exact. However, in practical simulations one always considers a finite number of basis vectors:

\hat{f_k}= \sum_{n=0}^{N-1} f_n e^{-ikn/N} \quad k=0,1,...,N-1 \quad(1)

where $f_k$ is the $k-$ th sampling of the objective function and $\hat{f_k}$ is the set of Fourier coefficients, with which the original function can be constructed by inverse Fourier transform

f_k = \frac 1 N \sum_{n=0}^{N-1} \hat{f_n} e^{ikn/N}.

A naive implementation of a Discrete Fourier Transform would require $O(N^2)$ operations to compute. This becomes clear if we express equation (1) in matrix form:

\begin{pmatrix}\hat{f_0} \\ \hat{f_1} \\ \hat{f_2} \\ \vdots \\ \hat{f_{N-1}} \end{pmatrix} =\begin{pmatrix} 1 & 1 & 1 &\dots &1 \\ 1 & w_N & w_N^2 & \dots & w_N^{N-1}\\ 1 & w_N^2 & w_N^4 & \dots & w_N^{2(N-1)}\\ \vdots & & & & \vdots \\ 1 & w_N^{N-1} & w_N^{2(N-1)} & \dots & w_N^{(N-1)^2} \end{pmatrix} \begin{pmatrix} f_0 \\ f_1 \\ f_2 \\ \vdots \\ f_{N-1} \end{pmatrix}

where $w_N = e^{- i /N}$ as the unit frequency component. The Fast Fourier transform (FFT) algorithm, instead, allows us to find $\hat{f_k}$ with only $O(N\log N)$ complexity, making it applicable to large systems. One of the most common FFT is Cooley-Tukey method [Math. Comp. 19, 297–301 (1965)], which uses a divide-and-conquer approach by recursively breaking down the matrix multiplication into two smaller parts. This method is used in many computational tools including the ones that support this notebook.