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
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
|
|
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
and sampled on
equidistant points. These points are denoted by
.
The potential can then be written
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)
If the boundary conditions are appropriately chosen, the discrete Fourier transforms
for the calculation of
can be performed analytically [76].
This is possible for the limiting case where
and the boundary
conditions are on a sphere of radius
for the cluster. For a one-dimensional
system we choose a torus of radius
and for the two-dimensional system
a slab of thickness
. The electrostatic potential for these systems are
listed in Table 2,
where
and
and
are the
Bessel functions of the first and second kind of integer order
.
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.