Finite Differences

The most straightforward approach to the evaluation of the Laplacian operator applied to a function at every grid point is to approximate the second derivative by finite differences of increasing order of accuracy [7]. For example, the $\partial^2 \phi / \partial x^2$ part of the Laplacian on a grid of orthorhombic symmetry is

$$\frac{\partial^2 \phi}{\partial x^2} (x_i, y_j, z_k) \sime... ...{A/2} C_n^{(A)} \phi(x_i+nh_{x}, y_j, z_k) + O(h_{x}^{A}) \ ,$$ (5)

where $h_{x}$ is the grid spacing in the $x$-direction, $A$ is the order of accuracy and is an even integer, and the weights $C_n^{(A)}$ are even with respect to $n$, i.e. $C_n^{(A)}=C_{-n}^{(A)}$. This equation is exact when $\phi$ is a polynomial of degree less than or equal to $A$. The leading contribution to the error is of order $h_{x}^{A}$. The full Laplacian operator for a single grid point in three dimensions consists of a sum of $(3A+1)$ terms.

In principle, for well behaved functions, the second order form of equation (5) should converge to the exact Laplacian as $h \rightarrow 0$. Therefore to increase the accuracy of a calculation one would need to proceed to smaller grid spacings. However, in most cases of interest, this is computationally undesirable and instead, formulae of increasing order are used to improve the accuracy at an affordable cost [8]. Chelikowsky et al. [9], in their finite difference pseudopotential method, have tested the finite difference expression for up to $A=18$ on calculations of a variety of diatomic molecules and have suggested $A=12$ as the most appropriate for their purpose, as the higher orders did not provide any significant improvement.

Alternative discretisations of the Laplacian operator are possible, such as the Mehrstellen discretisation of Briggs et al. [5]. This is a fourth order discretisation that includes off-diagonal terms, but only nearest neighbours to the point of interest. It is more costly to compute than the standard fourth order formula of equation (5) and it is still not clear whether its fourth order is sufficient. One may also use FD methods on a grid with variable spatial resolution, such as that of Modine et al. [10] which is denser near the ionic positions. Such a scheme, however, has the added overhead of a transformation of the Laplacian from Cartesian to curvilinear coordinates. In this paper we use only the FD scheme of equation (5).

The FD approach has desirable properties, both in terms of computational scaling and parallelisation. The Laplacian in the FD representation is a near-local operator, becoming more delocalised with increasing order. Therefore, the cost of applying it to $N$ grid points is strictly linear (compared to $N\log N$ for Fourier transform methods). Also, as a result of its near-locality, ideal load balancing can be achieved in parallel implementations by partitioning the real space grid into subregions of equal size and distributing them amongst processing elements (PEs) while requiring little communication for applying the Laplacian at the bordering points of the subregions.

If $N_s$ represents the size of the system, then the number of support functions will be proportional to $N_s$ and so will the total number of grid points in the simulation cell, resulting in a total computational cost proportional to $N_s^2$ for the application of the Laplacian on all support functions. More favourable scaling can be achieved by predicting the region in space whithin which the values of a particular function will be of significant magnitude and operating only on this region [4,11]. Linear-scaling can be achieved by strictly restricting from the outset the support functions to spherical regions centred on atoms [12]. In this case, the cost is $qN_s$ with $q$ being the cost of applying the Laplacian on the points of a spherical region, which is constant with system size.

FD methods nevertheless have disadvantages that do not appear in the plane-wave formalism. Firstly, there is no a priori way of knowing whether a particular order of FD approximation will be sufficient to represent a particular support function accurately. In addition, while plane-wave methods can handle different symmetry groups trivially through the reciprocal lattice vectors of the simulation cell, real space implementations need to consider every symmetry separately and require considerable modifications to the code and higher computational cost. Briggs et al. [5] have demonstrated this difficulty by performing calculations with hexagonal grids while most common applications of real space methods in the literature are limited to grids of cubic or orthorhombic symmetry [4,6,9,12].

The computational cost for the calculation of the Laplacian of a single support function with the FD method scales as $(3A+1)(1+A/D)^3 N_{\mathrm{reg}}$ where $N_{\mathrm{reg}}$ is the number of grid points within the support region, and $D$ is the number of grid points along the support region diameter and is proportional to $N_{\mathrm{reg}}^{1/3}$. This estimate of cost includes all the nonzero values of the Laplacian, which in general occur not only at the grid points inside the support region but also at points outside, up to a distance of $A/2$ points from the region's boundary. It is important to include the contribution to the Laplacian from outside the support region in the sum of equation (4) in order to obtain the best possible accuracy for a given order $A$ and also to ensure the Hermiticity of the discretised representation of the Laplacian, $\hat{T}$, and hence of the kinetic energy matrix elements $T_{\alpha \beta}$.