PSCF v1.1
Appendix: Periodic Functions and Fourier Series

Appendix: Self-Consistent Field Theory (Prev)         Appendix: Space Group Symmetry (Next)

This page describes the different numerical representations used in the pscf_pc and pscf_pg programs to describe fields (functions of position) that have the periodicity of a Bravais lattice. In a subsequent appendix, we also consider a type of symmetry-adapted Fourier expansion appropriate for periodic fields that are known to be invariant under all symmetry operations of a specified space group.

Below, we give succinct descriptions of Bravais lattices and reciprocal lattices, emphasizing aspects that are relevant to the PSCF software. We also briefly review the use of Fourier series for continuous periodic fields, the discrete Fourier transform (DFT) of periodic fields that are defined on the nodes of a mesh, and the notion of wavevector aliasing. Readers who have not previously encountered some of these topics may find it useful to also consult more detailed treatments. For a general overview of topics related to periodic structures and crystal symmetry, we recommend the book Structure of Materials: An Introduction to Crystallography, Diffraction and Symmetry by Marc De Graef and Michael E. McHenry. For a discussion of Fourier methods and a summary of the numerical benefits of solving differential equations in reciprocal space (sometimes called "spectral accuracy"), we recommend the first four chapters of the book Spectral Methods in MATLAB by Lloyd N. Trefethen.

Throughout these notes, let \( D \) denote the number of directions in which the fields of interest are taken to be periodic, which is the dimensionality of the space in which these fields are defined. Allowed values of \( D \) are \( D = 1 \) for lamellar structures, \( D = 2 \) for structures such as those involving two-dimensionally periodic arrangements of cylinders (e.g., hexagonally packed cylinders), and \( D = 3 \) for crystals with full three-dimensional periodicity (e.g., arrangements of spheres and 3D network structures).

Bravais Lattice

Consider a D-dimensional crystal with a Bravais lattice defined by lattice basis vectors \( {\bf a}_{0}, \ldots, {\bf a}_{D-1} \). The Bravais lattice is defined to be the infinite set containing every vector \( {\bf R} \) that can be expressed as a sum

\[ {\bf R} = \sum_{\alpha=0}^{D-1} m_{\alpha} {\bf a}_{\alpha} \]

with arbitrary integer values of the indices \( m_{0},\ldots,m_{D-1} \). A function \( f({\bf r}) \) of position \( {\bf r} \) is said to have the periodicity of this lattice if it is invariant under translation by any Bravais lattice vector, i.e., if

\[ f({\bf r} + {\bf R}) = f({\bf r}) \]

for any D-dimensional position vector \( {\bf r} \) and any translation vector \( {\bf R} \) that belongs to the Bravais lattice.

Reciprocal Lattice

Let \( {\bf b}_{0}, \ldots, {\bf b}_{D-1} \) denote \( D \) reciprocal lattice basis vectors that are defined such that

\[ {\bf a}_{\alpha} \cdot {\bf b}_{\beta} = 2\pi \delta_{\alpha \beta} \]

for \( \alpha, \beta \in [0,..,D-1] \), where \( \delta_{\alpha \beta}\) is a Kronecker delta function. The reciprocal lattice associated with a Bravais lattice is defined to be the infinite set containing every wavevector \( {\bf G} \) that can be expressed as a sum

\[ {\bf G}(n_{0}, \ldots, n_{D-1}) = \sum_{\alpha=0}^{D-1} n_{\alpha} {\bf b}_{\alpha} \]

with integer values of \( n_{0}, \ldots, n_{D-1} \). In what follows, we will sometimes refer to the integer indices \( n_{0}, \ldots, n_{D-1}\) that appear in this expansion of a reciprocal lattice vector as Miller indices for that wavevector.

Any continuous function \( f({\bf r}) \) of position \( {\bf r} \) with the periodicity of a Bravais lattice may be expanded as a Fourier series

\[ f({\bf r}) = \sum_{{\bf G}} \tilde{f}({\bf G}) e^{i {\bf G} \cdot {\bf r}} \quad, \]

in which \( {\bf G} \) denotes a reciprocal lattice vector, the sum is taken over the infinite set of all reciprocal lattice vectors, and \( \tilde{f}({\bf G}) \) is a complex Fourier coefficient associated with wavevector \( {\bf G} \). The Fourier coefficient \( \tilde{f}({\bf G}) \) is equal to the value of the integral

\[ \tilde{f}({\bf G}) = \frac{1}{V_{\rm cell}} \int\limits_{{\rm cell}} d^{D}r \; e^{-i {\bf G} \cdot {\bf r}} f({\bf r}) \quad, \]

in which the integral is a D-dimensional integral over one unit cell, and \( V_{\rm cell} \) is the generalized volume of the unit cell (i.e., the unit cell length for D=1, area for D=2, or volume for D=3).

The type of Fourier series defined above can be used to expand either real- or complex-valued periodic functions of position. If \( f \) is a real-valued function, then its Fourier coefficients must satisfy the condition

\[ \tilde{f}(-{\bf G}) = \tilde{f}^{*}({\bf G}) \]

for any reciprocal lattice vector \( {\bf G} \). Here and hereafter, we use \( C^{*} \) to denote the complex conjugate of any complex number \( C \). Note that this constraint requires that that \( \tilde{f}({\bf G}) \) be real in the specal case \( {\bf G} = 0 \).

Reduced Coordinates

It is convenient to introduce dimensionless reduced coordinates \( r_{0}, \ldots, r_{D-1} \) associated with any Cartesian position vector \( {\bf r} \) within a crystal, which are defined such that

\[ {\bf r}(r_0, \ldots, r_{D-1}) = \sum_{\alpha = 0}^{D-1} {\bf a}_{\alpha} r_{\alpha} \quad. \]

The reduced coordinates are thus the components of \( {\bf r} \) in a basis of Bravais lattice basis vectors. The primary unit cell of a crystal is defined to be a region in which \( 0 \leq r_{\alpha} < 1 \) for all \( \alpha = 0, \ldots, D-1 \). This corresponds to a unit domain in the space of reduced coordinates that is give by the domain \( [0,1) \) for D=1, a unit square for D=2 or a unit cube for D=3. This maps onto a unit cell in Cartesian coordinates that is a parallelogram for D=2 or a parallelepiped for D=3, and for which the vectors connecting neighboring corners are Bravais lattice basis vectors.

Suppose that \( f({\bf r}) \) is a periodic function of the Cartesian position vector \( {\bf r} \) and that \( F(r_{0}, \ldots, r_{D-1}) \) is an equivalent function of the reduced coordinates \( r_{0}, \ldots, r_{D-1} \), defined such that

\[ F(r_{0}, \ldots, r_{D-1}) = f({\bf r}(r_{0},\ldots,r_{D-1})) \]

for any list of reduced coordinates. The requirement that a function \( f({\bf r})\) be periodic under translations by Bravais lattice vector is equivalent to a requirement that \( F \) be invariant under changes of any reduced coordinate by addition of arbitrary integer. That is, for any periodic function \( F \) of the reduced coordinates,

\[ F(r_{0} + k_{0}, \ldots, r_{D-1} + k_{D-1} ) = F(r_{0}, \ldots, r_{D-1} ) \]

for arbitrary integer values of \( k_{0}, \ldots, k_{D-1} \).

Using the above expressions for \( {\bf r}(r_{0}, \ldots, r_{D-1}) \) and \( {\bf G}(n_{0}, \ldots, n_{D-1}) \) and the bi-orthogonality relation for Bravais and reciprocal basis vectors, it is straightforward to show that

\[ \exp \left ( i{\bf G} \cdot {\bf r} \right ) = \exp \left ( 2\pi i \sum_{\alpha = 0}^{D-1} n_{\alpha} r_{\alpha} \right ) \]

for any choice of reciprocal lattice vector \( {\bf G}(n_{0}, \ldots, n_{D-1}) \) and any position vector \( {\bf r}(r_{0}, \ldots, r_{D-1}) \). Upon substituting this into the general Fourier series expansion for a periodic function \( f({\bf r})\) given above, we find that the corresponding function \( F(r_{1}, \ldots, r_{D-1}) \) of the reduced coordinates may be expressed as an infinite series

\[ F(r_{0}, \ldots, r_{D-1}) = \sum_{n_{0} = -\infty}^{\infty} \cdots \sum_{n_{D-1} = -\infty}^{\infty} \tilde{F}(n_{0},\ldots,n_{D-1}) \exp \left ( 2\pi i \sum_{\alpha = 0}^{D-1} n_{\alpha} r_{\alpha} \right ) \]

in which the sum is taken over all integer values of each of the Miller indices, and in which \( \tilde{F}(n_{0}, \ldots, n_{D-1}) \) is a Fourier amplitude defined such that

\[ \tilde{f}({\bf G}(n_{0},\ldots,n_{D-1})) \equiv \tilde{F}(n_{0}, \ldots, n_{D-1}) \]

for the wavevector \( {\bf G}(n_{0}, \ldots, n_{D-1}) \) with integer components \(n_{0}, \ldots, n_{D-1}\). This form of the Fourier series is also a special case of the general Fourier expansion of a periodic function, as applied here to a function \( F(r_{0}, \ldots, r_{D-1}) \) that has the periodicity of a generalized cubic lattice with orthogonal Bravais lattice basis vectors of unit magnitude.

Mesh Discretization

All fields in the pscf_pc and pscf_pg programs are represented internally using values of fields on a regular D-dimensional mesh in the space of the reduced coordinates \( r_{0}, \ldots, r_{D-1} \). Let \( M_{\alpha} \) denote the number of mesh points per unit cell along direction \( \alpha \) for each integer direction index \( \alpha = 0, \ldots, D-1 \). The values of \( M_{0}, \ldots, M_{D-1} \) are specified in the parameter file for a pscf_pc or pscf_pg program as values for elements of the parameter "mesh", whose value is a vector with integer-valued elements. Nodes on this mesh correspond to points at which each reduced coordinate \( r_{\alpha} \) with \( \alpha \in [0, \ldots, D-1]\) has a value

\[ r_{\alpha} = m_{\alpha}/M_{\alpha} \]

for some integer \( m_{\alpha} \in [0,\ldots,M_{\alpha}-1] \) for nodes within the primary unit cell.

Discrete Fourier Transform (DFT)

Let \( F \) denote a function that is defined only on the nodes of such a regular mesh, where \( F(m_{0}, \ldots, m_{D-1}) \) denotes the value of \( F \) on the mesh node with integer indices \( m_{0}, \ldots, m_{D-1} \).

Forward DFT

Let \( \tilde{F} \) denote the discrete Fourier transform (DFT) of \( F \), which we define here by a sum

\[ \tilde{F} (n_0, \ldots, n_{D-1}) = \frac{1}{M_{\rm tot}} \sum_{m_{0}=0}^{M_{0}-1} \cdots \sum_{m_{D-1}=0}^{M_{D-1}-1} F(m_{0},\ldots,m_{D-1}) \exp \left ( -2\pi i \sum_{\alpha = 0}^{D-1} n_{\alpha} m_{\alpha}/M_{\alpha} \right ) \]

in which \( M_{\rm tot} \) denotes the total number of mesh points in the primary unit cell of the mesh, given by the product

\[ M_{\rm tot} = M_{0} \cdots M_{D-1} \quad. \]

of the numbers of mesh points in each direction.

Inverse DFT

A periodic function \( F \) that is defined on the nodes of mesh can can be expanded in terms of its DFT \( \tilde{F} \) by applying an inverse discrete Fourier transform

\[ F(m_0, \ldots, m_{D-1}) = \sum_{n_{0}=0}^{M_{0}-1} \cdots \sum_{n_{D-1}=0}^{M_{D-1}-1} \tilde{F}(n_{0},\ldots,n_{D-1}) \exp \left ( 2\pi i \sum_{\alpha = 0}^{D-1} n_{\alpha} m_{\alpha}/M_{\alpha} \right ) \]

Here, as elsewhere in this discussion \( F(m_0, \ldots, m_{D-1}) \) denotes the value of the function \( F \) at the mesh point with integer mesh indices \( m_{0}, \ldots, m_{D-1} \). Note that the sums over the indices \( n_{0}, \ldots, n_{D-1} \) in an inverse discrete DFT extend only over a finite number of values \( n_{\alpha} = 0, \ldots, M_{\alpha}-1 \) for each index, giving a finite rather than an infinite sum.

Relationship to source code: The D-dimensional forward and inverse DFT, defined using the conventions given above, are implemented for the field container classes used by pscf_pc programs by the member functions forwardTransform and inverseTransform defined in the class Pscf::Pspc::FFT<D>, respectively. The class Pscf::Pspc::FFT<D> is an instantiation of the class template Pscf::Pspc::FFT , in which D=1, 2, or 3 is a template parameter that represents the dimensionality of space. These member functions are simple wrappers around functions provided by the open source FFTW package. Equivalent operations are implemented for the field container classes used by the pscf_pg programs by corresponding member functions of the class template Pscf::Pspg::FFT that are wrappers around functions provided by the cuFFT CUDA library, which implements FFTs on a GPU.

Wavevector Aliasing

The discrete Fourier transform \( \tilde{F} (n_0, \ldots, n_{D-1}) \) is a function of \( D \) integer indices, and is thus itself defined on a lattice in reciprocal space. It is straightforward to show that the DFT \( \tilde{F} \) defined above is itself a periodic function in reciprocal space, for which \( \tilde{F}(n_0, \ldots, n_{D-1}) \) is invariant under a shift of the index \( n_{\alpha} \) associated with direction \( \alpha \) by any multiple of the number of mesh points \( M_{\alpha} \) along direction \( \alpha \). That is, for any DFT \( \tilde{F} \), we have

\[ \tilde{F}(n_{0} + k_{0}M_{0}, \ldots, n_{D-1} + k_{D-1}M_{D-1} = \tilde{F}(n_{0}, \ldots, n_{D-1}) \]

for any choice of integers \( k_{0}, \ldots, k_{D=1} \).

Proof : To prove this statement, we simply note that the effect of changing \( n_{\alpha}\) to \( n_{\alpha} + k_{\alpha} M_{\alpha} \) within the sum that defines \( \tilde{f} \) is to shift the argument of the exponential by a corresponding integer multiple \( m_{\alpha}k_{\alpha}2\pi i \) of \( 2\pi i \), and that changing the argument of an exponential by an integer multiple of \( 2\pi i \) leaves the value of the exponential unchanged.

This periodicity of the value of any discrete Fourier transforms is sometimes referred to as "aliasing" of wavevectors. Any list of D Miller integer indices \( n_{0}, \ldots, n_{D-1} \) can be used to construct an associated reciprocal lattice wavevector \( {\bf G}(n_0, \ldots, n_{D-1}) \) in which the Miller indices are components of the vector in a basis of reciprocal lattice basis vectors, as above. Two lists of integer Miller indices are referred to as "aliases" of one another if these lists contain values for each index \( n_{i} \) that are either equal or that differ by some nonzero integer multiple of the number \( M_{i} \) of mesh points along direction \( i \), for each \( i = 0, \ldots , D-1 \) . Reciprocal lattice wavevectors for which the Miller indices are aliases are also said to be aliases of one another.

Reciprocal lattice vectors that are aliases of one another yield values for the complex exponential \( \exp( i {\bf G}\cdots {\bf r} ) \) that are equal for any position vector \( {\bf r} \) that is a node of the relevant mesh, but that are generally unequal if evaluated for values of \( {\bf r} \) that do not lie on this mesh. Wavevectors that are aliases of one another are thus treated as equivalent within forward and inverse DFT operations, in which values of the underlying field are only defined on the nodes of a mesh.

Wavevector Indexing Schemes

The phenomena of aliasing of wavevectors for fields that are defined on a mesh means that there are infinite number of ways of defining Miller indices so as to generate aliases of a particular wavevector, by adding or subtracting integer multiples of the number of mesh points to each index. Two different conventions are used in PSCF assign a unique list of indices to each wavevector within a DFT. We refer to the resulting conventions for lists of indices here as (1) discrete Fourier transform or DFT indices and (2) Minimal indices.

Standard or DFT Indices

The DFT indices of a wavevector are defined such that the index \( n_{i} \) associated with reduced coordinate \( i \) is required to be in the range \( 0, \ldots, N_{i} - 1 \), where \( N_{i} \) is the number of mesh points in direction \( i \). This is the indexing scheme used in the above definitions of the forward and inverse DFT operations.

Minimal or Brillouin Zone (BZ) Indices

The minimal indices of a wavevector are indices that are chosen so as to give a corresponding Cartesian wavevector \( {\bf G}(n_{0}, \ldots, n_{D-1}) \) as defined above as an expansion in reciprocal lattice basis vectors for which the magnitude \( |{\bf G}| \) is minimal. The indices obtained using this convention are referred to in parts of the PSCF source code as "Brillouin zone" indices because the underlying geometrical construction is analogous to that used to define the first Brillouin zone for crystal wavevectors in band theories of electronic states in crystalline solids.

If two or more equivalent choices of indices for a wavevector yield wavevector aliases for which \( |{\bf G}| \) is equal to its minimum possible value, then the minimal indices are defined to be the choice which for which the indices are "maximum" when lists of indices are compared while giving priority to indices that appear earlier in a list \( (n_{0}, \ldots, n_{D-1}) \), thus ordering first by \( n_{0} \), then ordering choices with equal values of \( n_{0} \) in order of increasing \( n_{1} \), etc.

Relationship to source code: The standard (or DFT) and minimal (or Brillouin zone) indexes of each wavevector used in a DFT are assigned to each wavevector within the Pscf::Basis<D> class.


Appendix: Self-Consistent Field Theory (Prev)         User Guide (Up)         Appendix: Space Group Symmetry (Next)