Fast Fourier Transforms

A function given as a finite linear combination of plane waves can also be defined as a set of functional values on a equally spaced grid in real space. The sampling theorem (see e.g. Ref. [23]) gives the maximal grid spacing that still allows to hold the same information as the expansion coefficients of the plane waves. The real space sampling points ${\bf R}$ are defined

$${\bf R} = {\bf h} \: {\bf N} {\bf q} \enspace ,$$ (223)

where ${\bf N}$ is a diagonal matrix with the entries $1/N_{s}$ and ${\bf q}$ is a vector of integers ranging from 0 to $N_{s} - 1$ (s = x, y, z). To fulfill the sampling theorem $N_{s}$ has to be bigger than $2$   max$({\bf g}_{s}) + 1$ . To be able to use fast Fourier techniques, $N_{s}$ must be decomposable into small prime numbers (typically 2, 3, and 5). In applications the smallest number $N_{s}$ that fulfills the above requirements is chosen.

A periodic function can be calculated at the real space grid points

$$f({\bf R}) = \sum_{\bf G} f({\bf G}) \exp \! \left[ i {\bf G} \cdot {\bf R} \right]$$ (224)

$$= \sum_{\bf g} f({\bf G}) \exp \! \left[ 2 \pi i \left( ({\bf h}^{\rm t})^{-1} {\bf g}\right) \cdot \left({\bf h} {\bf N} {\bf q}\right) \right]$$ (225)

$$= \sum_{\bf g} f({\bf G}) \exp \! \left[ {2 \pi \over N_{\rm x}} i ... ...\exp \! \left[ {2 \pi \over N_{\rm z}} i g_{\rm z} q_{\rm z} \right] \enspace .$$ (226)

The function $f({\bf G})$ is zero outside the cutoff region and the sum over ${\bf g}$ can be extended over all indices in the cube $-{\bf g}^{\rm max}_{\rm s} \ldots {\bf g}^{\rm max}_{\rm s}$ . The functions $f({\bf R})$ and $f({\bf G})$ are related by three-dimensional Fourier transforms

$$f({\bf R}) = {\rm FT}^{-1} \left[ f({\bf G}) \right]$$ (227)

$$f({\bf G}) = {\rm FT} \left[ f({\bf R}) \right] \enspace .$$ (228)

The Fourier transforms are defined by

$$\left[ {\rm FT}^{-1} \left[ f({\bf G}) \right] \right]_{uvw} = \sum_{j=0}^{N_{\rm x}-1} \sum_{k=0}^{N_{\rm y}-1} \sum_{l=0}^{N_{\rm z}-1} f_{jkl}^{\bf G}$$

$$\exp \! \left[ i {2 \pi \over N_{\rm x}} j u \right] \exp \... ...}} k v \right] \exp \! \left[ i {2 \pi \over N_{\rm z}} l w \right]$$ (229)

$$\left[ {\rm FT} \left[ f({\bf R}) \right] \right]_{jkl} = \sum_{u=0}^{N_{\rm x}-1} \sum_{v=0}^{N_{\rm y}-1} \sum_{w=0}^{N_{\rm z}-1} f_{uvw}^{\bf R}$$

$$\exp \! \left[ -i {2 \pi \over N_{\rm x}} j u \right] \exp ... ... \right] \exp \! \left[ -i {2 \pi \over N_{\rm z}} l w \right] \enspace ,$$ (230)

where the appropriate mappings of ${\bf q}$ and ${\bf g}$ to the indices

$$[u,v,w] = {\bf q}$$ (231)

$$\{j,k,l \} = {\bf g}_s {\bf g}_s \geq 0$$ (232)

$$\{j,k,l \} = N_s + {\bf g}_s {\bf g}_s < 0$$ (233)

have to be used. From Eqs. (229) and (230) it can be seen, that the calculation of the three-dimensional Fourier transforms can be performed by a series of one dimensional Fourier transforms. The number of transforms in each direction is $N_{\rm x} N_{\rm y}$ , $N_{\rm x} N_{\rm z}$ , and $N_{\rm y} N_{\rm z}$ respectively. Assuming that the one-dimensional transforms are performed within the fast Fourier transform framework, the number of operations per transform of length $n$ is approximately $5 n \log{n}$ . This leads to an estimate for the number of operations for the full three-dimensional transform of $5 N \log{N}$ , where $N = N_{\rm x} N_{\rm y} N_{\rm z}$ .