Appendix A: Delta functions

In this paper, whenever we refer to `delta functions' we will assume a periodic and bandwidth limited version of the Dirac delta functions. These delta functions are three-dimensional versions of the `impulse functions' that are common in signal processing applications of FFTs [53]. In electronic structure, similar functions have been used as `mesh delta functions' in the `exact finite difference method' of Hoshi et al. [54] and in recent studies of their possible application when we consider the limit of an infinite simulation cell [55].

In our derivations we will assume that we have a simulation cell of any symmetry, which in general is a parallelepiped defined by its primitive lattice vectors $\mathbf{A}_1$, $\mathbf{A}_2$ and $\mathbf{A}_3$. In this simulation cell we define a regular grid with an odd number of points $N_1=2J_1+1$, $N_2=2J_2+1$, and $N_3=2J_3+1$ in every direction (the adaptation of our results to the case of even numbers of points is straightforward). Therefore point $\mathbf{r}_{KLM}$ of this regular grid is defined as

$$\mathbf{r}_{KLM}=\frac{K}{N_1}\mathbf{A}_1 + \frac{L}{N_2}\mathbf{A}_2 + \frac{M}{N_3}\mathbf{A}_3$$ (31)

with $K=0,1,\ldots(N_1-1)$, etc.

Bandwidth limited delta functions centred at points of the regular grid are defined as

$$D_{KLM}(\mathbf{r}) = D_{000}(\mathbf{r}-\mathbf{r}_{KLM})$$

$$= \frac{1}{N_1 N_2 N_3} \sum_{P= -J_1}^{J_1} \sum_{Q= -J_2}^{J_2} \... ...hbf{B}_1 + Q\mathbf{B}_2 + R\mathbf{B}_3 ) \cdot(\mathbf{r}-\mathbf{r}_{KLM}) }$$ (32)

where $\mathbf{B}_1$, $\mathbf{B}_2$ and $\mathbf{B}_3$ are the primitive reciprocal lattice vectors of the simulation cell. Plane-waves whose wavevector is a linear combination of these reciprocal lattice vectors have periodicity compatible with the simulation cell and therefore so do our delta functions, or any other function expanded in terms of these plane-waves. These periodic bandwidth limited delta functions are our basis set. A plot of a two-dimensional version of one of these delta functions is shown in Figure 7. It is obvious from (32) that the delta functions are real-valued everywhere in space. They are not normalised to unity but they are normalised to the grid point volume ($V$ is the volume of the simulation cell)

$$W=\frac{V}{N_1 N_2 N_3} \,\,.$$ (33)

Their value at grid points is equal to one when the grid point coincides with the centre of the function and zero for all other grid points

$$D_{KLM}(\mathbf{r}_{FGH})= \delta_{KF} \delta_{LG} \delta_{MH} \,\, .$$ (34)

Figure 7: A two-dimensional version of one of the functions that constitute our basis set. Here function $D_{00}(\mathbf{r})$ is shown which is identically equal to 1 at its centre (point $\mathbf{r}_{00}$) and equal to zero at the centres (shown as black dots in the picture) of all other functions in the basis set.

The delta functions act as Dirac delta functions with the added effect of filtering out any plane-wave components that are not part of them. For example, if $f(\mathbf{r})$ is a function periodic with the periodicity of the simulation cell but not bandwidth limited, it can be expressed in terms of its discrete Fourier transform (plane-wave) expansion

$$f(\mathbf{r})= \frac{1}{V} \sum_{S=-\infty}^{\infty} \sum_{... ...mathbf{B}_1 + T\mathbf{B}_2 + U\mathbf{B}_3 )\cdot \mathbf{r}}$$ (35)

where $V$ is the volume of the simulation cell.

It is straightforward to show that the projection of $f(\mathbf{r})$ onto $D_{KLM}(\mathbf{r})$ is

$$\begin{eqnarray*} & & \int_V D_{KLM}(\mathbf{r}) f(\mathbf{r}) d\mathbf{r} \\ &... ..._3 ) \cdot \mathbf{r}_{KLM} } \\ &=& W f_{D}(\mathbf{r}_{KLM}) \end{eqnarray*}$$

We define here $f_{D}(\mathbf{r})$ to be the bandwidth limited version of the function $f(\mathbf{r})$, limited to the same frequency components as $D_{KLM}(\mathbf{r})$.

As the NGWFs are linear combinations of the delta functions according to (7), the result of equation (36) is very important since it leads to the following relation

$$\int_V \phi_{\alpha}(\mathbf{r}) f(\mathbf{r}) d\mathbf{r} =... ...-1} \sum_{M=0}^{N_3-1} C_{KLM,\alpha} f_{D}(\mathbf{r}_{KLM})$$ (36)

which means that the integral in the lefthand side of the above equation is exactly equal to a discrete summation of values on the grid, provided we use the bandwidth limited version of $f(\mathbf{r})$.

As a corollary we observe that the delta functions are an orthogonal set since

$$\int_V D_{FGH}(\mathbf{r}) D_{KLM}(\mathbf{r}) d\mathbf{r} =... ...thbf{r}_{KLM}) = W \delta_{FK} \delta_{GL} \delta_{HM} \,\, .$$ (37)

We also need to define the fine grid delta functions $B_{XYZ}(\mathbf{r})$ (here the $XYZ$ are just grid point indices for the fine grid, they are not related to any Cartesian coordinates). These functions are the analogues of the delta functions we have just described that would be obtained if we doubled the minimum and maximum values that their wavevectors can take. Consequently, they have the same periodicity but they correspond to a grid with twice the number of points in every direction, i.e. $2N_1$, $2N_2$ and $2N_3$ points. They are defined by

$$B_{XYZ}(\mathbf{r}) = B_{000}(\mathbf{r}- \mathbf{r}_{XYZ})$$

$$= \frac{1}{8 N_1 N_2 N_3} \sum_{P=(-N_1+1)}^{N_1} \sum_{Q=(-N_2+1)}... ...thbf{B}_1 + Q\mathbf{B}_2 + R\mathbf{B}_3 ) \cdot(\mathbf{r}-\mathbf{r}_{XYZ})}$$ (38)

As expected, the fine grid delta functions also satisfy an equation similar to (36)

$$\int_V B_{XYZ}(\mathbf{r}) f(\mathbf{r}) d\mathbf{r} = \frac{W}{8}f_{B}(\mathbf{r}_{XYZ})$$ (39)

where $f_{B}(\mathbf{r})$ is again a bandwidth limited version of $f(\mathbf{r})$ but this time it is limited to contain any of the plane-waves that constitute $B_{XYZ}(\mathbf{r})$ rather than $D_{KLM}(\mathbf{r})$. It is easy to verify that any function that can be written as a sum of products of pairs of delta functions can also be written as a fine grid delta function expansion. We define and use the fine grid delta functions because of this `product rule' property.