Charged Systems

The possibility to use fast Fourier transforms to calculate the electrostatic energy is one of the reasons for the high performance of plane wave calculations. However, plane wave based calculations imply periodic boundary conditions. This is appropriate for crystal calculations but very unnatural for molecule or slab calculations. For neutral systems this problem is circumvented by use of the supercell method. Namely, the molecule is periodically repeated but the distance between each molecule and its periodic images is so large that their interaction is negligible. This procedure is somewhat wasteful but can lead to satisfactory results.

Handling charged molecular systems is, however, considerably more difficult, due to the long range Coulomb forces. A charged periodic system has infinite energy and the interaction between images cannot really be completely eliminated. In order to circumvent this problem several solutions have been proposed. The simplest fix-up is to add to the system a neutralizing background charge. This is achieved trivially as the ${\bf G} = {\bf0}$ term in the electrostatic energy is already eliminated. This leads to finite energies but does not eliminate the interaction between the images and makes the calculation of absolute energies difficult. Other solutions involve performing a set of different calculations on the system such that extrapolation to the limit of infinitely separated images is possible. This procedure is lengthy and one cannot use it easily in molecular dynamics applications. It has been shown, that it is possible to estimate the correction to the total energy for the removal of the image charges [67]. Still it seems not easy to incorporate this scheme into the frameworks of molecular dynamics. Another method [68,69,70] works with the separation of the long and short range parts of the Coulomb forces. In this method the low-order multipole moments of the charge distribution are separated out and handled analytically. This method was used in the context of coupling ab initio and classical molecular dynamics [71].

The long-range forces in the electrostatic energy are contained in the first term. This term can be written

$${1 \over 2} \int \! \! \int \! d{\bf r} d{\bf r'} { n_{\rm... ...\over 2} \int \! d{\bf r} V_{\rm H}({\bf r}) n_{\rm tot}({\bf r}) \enspace ,$$ (124)

where the electrostatic potential $V_{\rm H}({\bf r})$ is the solution of Poisson's equation. We will discuss two approaches to solve Poisson's equation subject to the boundary conditions $V_{\rm H}({\bf r}) \rightarrow 0$ for ${\bf r} \rightarrow \infty$ . Both of them rely on fast Fourier transforms, thus keeping the same framework as for the periodic case.

Figure 3: Schematic view of the Hockney method to decouple images in the electrostatic energy. The first line shows the artificially replicated system consisting of the original computational cell and an empty duplicated cell. The Green's function of this new periodic system is shown in the lower part of the figure.

The first method is due to Hockney [72] and was first applied to density functional plane wave calculations in Ref. [73]. In the following outline, for the sake of simplicity, a one-dimensional case is presented. The charge density is assumed to be non-zero only within an interval $L$ and sampled on $N$ equidistant points. These points are denoted by $x_{p}$ . The potential can then be written

$$V_{\rm H}(x_{p}) = {L \over N} \sum_{p'=-\infty}^{\infty} \! \! G(x_{p} - x_{p'}) n({x_{p'}})$$ (125)

$$= {L \over N} \sum_{p'=0}^{N} G(x_{p} - x_{p'}) n({x_{p'}})$$ (126)

for $p = 0,1,2,\ldots N$ , where $G(x_{p} - x_{p'})$ is the corresponding Green's function. In Hockney's algorithm this equation is replaced by the cyclic convolution

$$\tilde V_{\rm H}(x_{p}) = {L \over N} \sum_{p'=0}^{2N+1} \tilde G(x_{p} - x_{p'}) \tilde n({x_{p'}})$$ (127)

where $p = 0,1,2,\ldots 2N+1$ , and

$$\tilde n({x_{p}}) = \left\{ { \begin{array}{lll} n({x_{p}}) & 0 & \leq p \leq N \\ 0 & N & \leq p \leq 2 N +1 \end{array} } \right.$$ (128)

$$\tilde G(x_{p}) = G(x_{p}) -(N+1) \leq p \leq N$$ (129)

$$\tilde n({x_{p}}) = \tilde n(x_{p} + L)$$ (130)

$$\tilde G(x_{p}) = \tilde G(x_{p} + L)$$ (131)

The solution $\tilde V_{\rm H}(x_{p})$ can be obtained by a series of fast Fourier transforms and has the desired property

$$\tilde V_{\rm H}(x_{p}) = V_{\rm H}(x_{p}) 0 \leq p \leq N \enspace .$$ (132)

To remove the singularity of the Green's function at $x = 0$ , $G(x)$ is modified for small $x$ and the error is corrected by using the identity

$$G(x) = {1 \over x} \left[{x \over r_{\rm c}}\right] + {1 \over x} \left[{x \over r_{\rm c}}\right] \enspace ,$$ (133)

where $r_{\rm c}$ is chosen such, that the short-ranged part can be accurately described by a plane wave expansion with the density cutoff. In an optimized implementation Hockney's method requires the double amount of memory and two additional fast Fourier transforms on the box of double size. Hockney's method can be generalized to systems with periodicity in one (wires) and two (slabs) dimensions. It was pointed out [74] that Hockney's method gives the exact solution to Poisson's equation for isolated systems if the boundary condition (zero density at the edges of the box) are fulfilled.

A different, fully reciprocal space based method, that can be seen as an approximation to Hockney's method, was recently proposed [75]. The final expression for the Hartree energy is also based on the splitting of the Green's function in Eq. (133)

$$E_{\rm ES} = 2 \pi \Omega \sum_{\bf G} V_{\rm H}^{\rm MT}({\bf G}) n^\star_{\rm tot}({\bf G}) + E_{\rm ovrl} - E_{\rm self} \enspace .$$ (134)

The potential function is calculated from two parts,

$$V_{\rm H}^{\rm MT}({\bf G}) = \bar V_{\rm H}({\bf G}) + \tilde V_{\rm H}({\bf G}) \enspace ,$$ (135)

where $\tilde V_{\rm H}({\bf G})$ is the analytic part, calculated from a Fourier transform of erfc

$$\tilde V_{\rm H}({\bf G}) = {4 \pi \over G^2} \left( 1 - \exp \! \left[ - {G^2 r_{\rm c}^2 \over 4 } \right] \right) n({\bf G})$$ (136)

and $\bar V_{\rm H}({\bf G})$ is calculated from a discrete Fourier transform of the Green's function on an appropriate grid. The calculation of the Green's function can be done at the beginning of the calculation and has not to be repeated again. It is reported [75] that a cutoff of ten to twenty percent higher than the one employed for the charge density gives converged results. The same technique can also be applied for systems that are periodic in one and two dimensions.

If the boundary conditions are appropriately chosen, the discrete Fourier transforms for the calculation of $\bar V_{\rm H}({\bf G})$ can be performed analytically [76]. This is possible for the limiting case where $r_{\rm c} = 0$ and the boundary conditions are on a sphere of radius $R$ for the cluster. For a one-dimensional system we choose a torus of radius $R$ and for the two-dimensional system a slab of thickness $Z$ . The electrostatic potential for these systems are listed in Table 2, where $G_{\rm xy} = \left[ g_{\rm x}^2 + g_{\rm y}^2 \right]^{1/2}$ and $J_{n}$ and $K_{n}$ are the Bessel functions of the first and second kind of integer order $n$ .

Table 2: Fourier space formulas for the Hartree energy, see text for definitions.

Dim.	periodic	$(G^2/4\pi) V_{\rm H}({\bf G})$	$V_{\rm H}({\bf 0})$
0	-	$\left( 1 - \cos \left[ R G \right] \right) n({\bf G})$	$2 \pi R^2 n(0)$
1	z	$\left( 1 + R \left( G_{\rm xy} J_1 (R G_{\rm xy}) K_0(R g_{\rm z}) \right. \right.$
		$\left. \left. - g_{\rm z} J_0 (R G_{\rm xy}) K_1(R g_{\rm z}) \right) \right) n({\bf G})$	0
2	x, y	$\left( 1 - (-1)^{g_{\rm z}} \exp \left[ - G Z / 2 \right] \right) n({\bf G})$	0
3	x, y, z	n(G)	0

Hockney's method requires a computational box such that the charge density is negligible at the edges. This is equivalent to the supercell approach [77]. Practical experience tells that a minimum distance of about 3 Å of all atoms to the edges of the box is sufficient for most systems. The Green's function is then applied to the charge density in a box double this size. The Green's function has to be calculated only once at the beginning of the calculation. The other methods presented in this chapter require a computational box of double the size of the Hockney method as they are applying the artificially periodic Green's function within the computational box. This can only be equivalent to the exact Hockney method if the box is enlarged to double the size. In plane wave calculations computational costs grow linearly with the volume of the box. Therefore Hockney's method will prevail over the others in accuracy, speed, and memory requirements in the limit of large systems. The direct Fourier space methods have advantages through their easy implementation and for small systems, if not full accuracy is required, i.e. if they are used with smaller computational boxes. In addition, they can be of great use in calculations with classical potentials.