The Ewald method

We begin with the Fourier Transform of the charge distribution p which, for a crystal with reciprocal-space vectors G, is given by 

In the final expression for the Madelung energy, we must take care to exclude the contribution to the energy of the true Coulomb potential due to ion I evaluated at R/, to avoid infinities. This term arises from the l/ r - R/ part of the second sum in Eq. (F.17), which must be dropped. However, for convergence purposes, and in order to be consistent with the extension of the first sum in Eq. (F.17) over all values of I, we must include separately the term which has r = R/ = R/ in the argument of the error function in the second sum. In essence, we do not need to worry about infinities once we have smoothed the charge distribution to gaussians, so that we can extend the sums over all values of the indices I and J, but we must do this consistently in the two sums of Eq. (F.17). From the 

Putting this result together with the previous expressions, we obtain for the Madelung energy 

As a final step, we replace the summations over the positions of ions in the entire crystal by summations over Bravais lattice vectors R, R and positions of ions t/, tj within the primitive unit cell (PUC)  

Since for all Bravais lattice vectors R and all reciprocal lattice vectors G we have exp( iG R) = 1, the expression for the Madelung energy takes the form 

---

Fig. 6.21 The construction of a system of periodic cells in the Ewald method. (Figure adapted from Allen M P and D ] Tildesley 1987. Computer Simulation of Liquids, Oxford, Oxford University Press.)...

The fourth term on the right-hand side of eq. (11.3) is the electrostatic interaction (Coulomb s law) between pairs of charged atoms i and j, separated by distance r j. Since electrostatic interactions fall off slowly with r (only as r-1) they are referred to as long-range and, for an infinite system such as a periodic solid, special techniques, such as the Ewald method, are required to sum up all the electrostatic interactions (cf. Section 7.1) (see e.g. Leach, Jensen (Further reading)). The... 

Because of the long-range nature of the dipolar interaction, care must be taken in the evaluation of the dipolar field. For hnite systems the sums in Eq. (3.12) are performed over all particles in the system. Eor systems with periodic boundary conditions the Ewald method [57-59], can be used to correctly calculate the conditionally convergent sum involved. However, in most work [12,13] the simpler Lorentz-cavity method is used instead. 

The effeet of eharged niekel ions on the vibration spectrum of ZnO crystals was studied by ealeulating the S5mimetrized LVDS in perfect and defective crystals using a cluster of 1000 ions (region 1). For the calculation of the Coulombic component of the diagonal elements of the force constant matrixes of the ion-ion interactions by the Ewald method about 4000 ions were considered for region 2. 

The Ewald method splits the summation into a summation in direct space and a summation in reciprocal space ... 

In recent years, a number of models have been introduced which permit the inclusion of long-range electrostatic interactions in molecular dynamics simulation. For simulations of proteins and enzymes in a crystalline state, the Ewald summation is considered to be the correct treatment for long range electrostatic interactions (Ewald 1921 Allen and Tildesley 1989). Variations of the Ewald method for periodic systems include the particle-mesh Ewald method (York et al. 1993). To treat non-periodic systems, such as an enzyme in solution other methods are required. Kuwajima et al. (Kuwajima and Warshel 1988) have presented a model which extends the Ewald method to non-periodic systems. Other methods for treating explicitly long-range interactions... 

However, not all systems require periodic boundary conditions in all spatial directions.For membranes, for example, only two dimensions are periodic, while the third one is finite. In that case, the Ewald method is computationally highly inefficient and wouid not aiiow to treat more than a few hundred charged particles. We present two alternative approaches, the MMM2D and ELC methods, which allow for computational efficiency similar to the bulk case. It is also simple to adapt the MMM2D method for systems with only one periodic dimension. 

The Ewald method can also be formulated for partially periodic boundary conditions, i.e. systems where only one or two of the three spatial dimensions are periodic. These geometries are applied for example in simulations of membranes or nano-pores, where some dimensions are supposed to have finite extend. In this case, however, the method scales like 0 N ), which only allows the treatment of moderately large systems with a few hundred charges. 

In 3D periodic boundary conditions, the Ewald method can be further accelerated to a computation time scaling of 0 N log N) using Fast Fourier Transformations (FFT) by smearing the charges onto a regular mesh. Various variants of these mesh based methods exist, such as P M, PME or SPME, but P M is known to be the computationally optimal variant [8]. In addition, there are variants of the P M method for e.g. orthorhombic boundary conditions or dipolar systems. Here we will concentrate on the classical case of a cubic simulation box, and give a brief introduction into both the Ewald method and P M. 

Similar to the Ewald method, there are also error estimates for P M [14], allowing to tune the method for optimal computational efficiency. Assuming that the charges are homogeneously distributed in the simulation box and that the errors Xij originating from the pairwise interactions of the charges i and j are uncorrelated, one obtains for the force error AEi of the ith particle... 

The limit E exists and has as its value that of the Ewald method minus the dipole term [16]. Starting from this convergence factor approach, Strebel and Sperb constructed a method of computational order 0(A / ) or, with a more clever algorithm, (P(AllogAl), MMM [15]. Unlike the particle mesh Ewald methods, no mesh is introduced, so that no interpolation errors occur, and it is comparatively easy to find error estimates for the method. Consequently, this method allows much higher precisions compared to P M or (S)PME. 

Similar to the Ewald method in full periodic boundary conditions, MMM2D is only fast enough to be used for a couple of hundred charges. For larger systems, one needs a more efficient method, which in this case can be obtained by a combination of a standard method for full 3d periodicity, e.g. P M, and the far formula of MMM2D. 

In summing this term over the lattice, in the calculations described subsequently, the Ewald method (Born and Huang, 1954) had to be invoked. 

The MD calculation to study supercritical fluid extraction from ceramics was performed with the XDORTO program developed by Kawamura[12]. The Verlet algorithm was used for the calculation of atomic motions, while the Ewald method was applied to the calculation of electrostatic interactions. Temperature was controlled by means of scaling the atom velocities under 3-dimensional periodic boundary condition. The calculations were made for 40000 time steps with the time increment of 2.5 10 seconds. The two body central force interaction potential, as shown in equation (3) was used for all the calculations ... 

For the derivation of the Ewald method it is important to consider the charge density pi r) corresponding to the electrostatic potential (vi) given in Eq. (6.2). The link between these quantities is provided by Poisson s equation [242]... 

The coefficients Qjj> ( ) for given k only depend on the lattice structure (their explicit form, which is not needed here, can be found in the book by Born and Huang (16), eqn (30.31) a detailed description of the Ewald method of separation of the field is also given there). These coefficients for k —> 0 tend to values which do not depend on the direction s = k/k. It follows from formulas (2.40) and (2.42) that the quantity LQ/j(k) and, in consequence, the eigenvalues of the matrix .np (see (2.19)]) are, in general, nonanalytic functions of k. 

---