The smooth particle mesh Ewald method

Before we discuss the non bonded multiple time step separation it is useful to describe in some details one of the most advanced techniques to handle long range forces. Indeed, this type of non bonded forces are the most cumbersome to handle and deserve closer scrutiny.

In the recent literature, a variety of techniques are available to handle the problem of long range interactions in computer simulations of charged particles at different level of approximation [30,31,12]. In this section, we shall focus on the Ewald summation method for the treatment of long range interactions in periodic systems [32,33,101]. The Ewald method gives the exact result for the electrostatic energy of a periodic system consisting of an infinitely replicated neutral box of charged particles. The method is the natural choice in MD simulations of complex molecular system with PBC.

The Ewald potential [33] is given by

$$V_{\rm qd}' = {1 \over 2} \sum_{ij}^{N} q_i q_j \sum_n { 1 \over ... ...erfc} \left ( \alpha \left \vert {\bf r}_{ij} + {\bf r}_n \right \vert \right )$$ (4.20)

$$V_{\rm qr} = \left [ {1 \over 2 \pi V} \sum_{{\bf m} \ne 0}^{\inf... ... {\bf m} \right ) - {\alpha \over \pi^{1/2}} \sum_i q_i^2 \right ] - V_{intra}.$$ (4.21)

with

$$S \left ( {\bf m} \right ) = \sum_{i}^{N} q_i \exp \left ( 2 \pi i {\bf m} \cdot {\bf r}_i \right )$$ (4.22)

$$V_{\rm intra} = \sum_{ij-\textrm{excl.}} q_{i}q_{j} { {\rm erf} \left ( \alpha r_{ij} \right ) \over r_{ij}},$$ (4.23)

where, ${\bf r_{i}}$ is the vector position of the atomic charge $q_{i}$, ${\bf r}_{ij} = {\bf r}_{i} - {\bf r}_{j}$, ${\bf r}_{n}$ is a vector of the direct lattice, ${\rm erfc}(x) = \pi^{-1/2}\int_{x}^{\infty} e^{-t^{2}} dt$ is the complementary error function, ${\rm erf}(x) = 1-erfc(x)$, $V$ the unit cell volume, ${\bf m}$ a reciprocal lattice vector and $\alpha$ is the Ewald convergence parameter. In the direct lattice part, Eq. (4.20), the prime indicates that intramolecular excluded contacts^4.1are omitted. In addition, in Eq. (4.21) the term $V_{\rm intra}$ subtracts, in direct space, the intra-molecular energy between bonded pairs, which is automatically included in the right hand side of that equation. Consequently, the summation on $i$ and $j$ in Eq. (4.23) goes over all the excluded intra-molecular contacts. We must point out that in the Ewald potential given above, we have implicitly assumed the so-called ``tin-foil'' boundary conditions: the Ewald sphere is immersed in a perfectly conducting medium and hence the dipole term on the surface of the Ewald sphere is zero [33]. For increasingly large systems the computational cost of standard Ewald summation, which scales with $N^2$, becomes too large for practical applications. Alternative algorithms which scale with a smaller power of $N$ than standard Ewald have been proposed in the past. Among the fastest algorithms designed for periodic systems is the particle mesh Ewald algorithm (PME)[34,35], inspired by the particle mesh method of Hockney and Eastwood [36]. Here, a multidimensional piece-wise interpolation approach is used to compute the reciprocal lattice energy, $V_{\rm qr}$, of Eq. 4.21, while the direct part, $V_{\rm qd}$, is computed straightforwardly. The low computational cost of the PME method allows the choice of large values of the Ewald convergence parameter $\alpha$, as compared to those used in conventional Ewald. Correspondingly, shorter cutoffs in the direct space Ewald sum $V_{\rm qd}$ may be adopted. If ${\bf u_{j}}$ is the scaled fractional coordinate of the $i$-th particle, the charge weighted structure factor, $S \left ( {\bf m} \right )$ in Eq. (4.22), can be rewritten as:

$$S \left ( {\bf m} \right ) = \sum_{j=1}^{N} q_{j} \exp \left [ 2 ... ...{1}} + {m_{2}u_{2j} \over K_{2}} + {m_{3}u_{3j} \over K_{3}} \right ) \right ].$$ (4.24)

Where, $N$ is the number of particles, $K_{1}, K_{2}, K_{3}$ and $m_1$, $m_2$, $m_3$ are integers. The $\alpha$ component of the scaled fractional coordinate for the $i$-th atom can be written as:^4.2

$$u_{i\alpha} = K_{\alpha} {\bf k}_{\alpha} \cdot {\bf r_i},$$ (4.25)

where ${\bf k}_{\alpha}$, $\alpha = 1,2,3$ are the reciprocal lattice basic vectors.

$S \left ( {\bf m} \right )$ in Eq. (4.24) can be looked at as a discrete Fourier transform (FT) of a set of charges placed irregularly within the unit cell. Techniques have been devised in the past to approximate $S \left ( {\bf m} \right )$ with expressions involving Fourier transforms on a regular grid of points. Such approximations of the weighted structure factor are computationally advantageous because they can be evaluated by fast Fourier transforms (FFT). All these FFT-based approaches involve, in some sense, a smearing of the charges over nearby grid points to produce a regularly gridded charge distribution. The PME method accomplishes this task by interpolation. Thus, the complex exponential $\exp (2 \pi i m_{\alpha}u_{i\alpha}/K_{\alpha})$, computed at the position of the $i$-th charge in Eq. (4.24), are rewritten as a sum of interpolation coefficients multiplied by their values at the nearby grid points. In the smooth version of PME (SPME) [35], which uses cardinal B-splines in place of the Lagrangian coefficients adopted by PME, the sum is further multiplied by an appropriate factor, namely:

$$\exp \left (2 \pi i m_{\alpha}u_{i\alpha}/K_{\alpha} \right ) = b... ..._{k} M_{n}(u_{i\alpha}-k) \exp \left ( 2 \pi i m_{\alpha}k/K_{\alpha} \right ),$$ (4.26)

where $n$ is the order of the spline interpolation, $M_{n}(u_{i\alpha}-k)$ defines the coefficients of the cardinal B-spline interpolation at the scaled coordinate $u_{i\alpha}$. In Eq. (4.26) the sum over $k$, representing the grid points, is only over a finite range of integers, since the functions $M_{n}(u)$ are zero outside the interval $0 \le u \le n$. It must be stressed that the complex coefficients $b(m_{i})$ are independent of the charge coordinates ${\bf u_i}$ and need be computed only at the very beginning of a simulation. A detailed derivation of the $M_{n}\left ( u \right )$ functions and of the $b_{\alpha}$ coefficients is given in Ref. [35]. By inserting Eq. (4.26) into Eq. (4.24), $S \left ( {\bf m} \right )$ can be rewritten as:

$$S({\bf m}) = b_{1}(m_{1}) b_{2}(m_{2}) b_{3}(m_{3}) {\cal F} \left [ Q \right ] \left ( m_{1},m_{2},m_{3} \right ),$$ (4.27)

where ${\cal F} \left [ Q \right ] \left ( m_{1},m_{2},m_{3} \right )$ stands for the discrete FT at the grid point $m_{1},m_{2},m_{3}$ of the array $Q \left ( k_{1},k_{2},k_{3} \right )$ with $1 \le k_{i} \le K_{i}$, $i=1,2,3$. The gridded charge array, $Q \left ( k_{1},k_{2},k_{3} \right )$, is defined as:

$$Q \left ( k_{1},k_{2},k_{3} \right ) = \sum_{i=1,N} q_i M_n \left... ..._1 \right ) M_n \left ( u_{i2} - k_2 \right ) M_n \left ( u_{i3} - k_3 \right )$$ (4.28)

Inserting the approximated structure factor of Eq. (4.27) into Eq. (4.21) and using the fact that ${\cal F}\left [ Q \right ] \left ( -m_{1}, -m_{2}, -m_{3} \right ) = K_{1} K_{2} K_{3} {\cal F}^{-1} \left [ Q \right ] \left ( m_{1} , m_{2},m_{3} \right )$, the SPME reciprocal lattice energy can be then written as

$$V_{\rm qr} = {1 \over 2 } \sum_{m_{1}=1}^{K_{1}}\sum_{m_{2}=1}^{K_{2}}\sum_{m_{3}=1}^{K_{3}} B (m_{1},m_{2},m_{3}) C(m_{1},m_{2},m_{3}) \times$$

$$\times {\cal F} \left [ Q \right ] (m_{1},m_{2},m_{3}) {\cal F} \left [ Q \right ] (-m_{1},-m_{2},-m_{3})$$ (4.29)

$$= {1 \over 2 } \sum_{m_{1}=1}^{K_{1}}\sum_{m_{2}=1}^{K_{2}}\sum_{m_... ...ht ] (m_{1},m_{2},m_{3}) {\cal F} \left [ Q \right ] (m_{1},m_{2},m_{3}) \times$$

$$\times K_{1}K_{2}K_{3} {\cal F}^{-1} \left [ Q \right ] (m_{1},m_{2},m_{3}),$$ (4.30)

with

$$B \left ( m_1, m_2, m_3 \right ) = \vert b_{1}(m_{1})\vert^{2}\vert b_{2}(m_{2})\vert^{2}\vert b_{3}(m_{3})\vert^{2}$$ (4.31)

$$C \left ( m_1, m_2, m_3 \right ) = (1/\pi V) \exp (-\pi^{2} {\bf m}^{2}/\alpha^{2})/{\bf m}^{2}$$ (4.32)

$$\Theta_{rec} = {\cal F}\left [ BC \right ].$$ (4.33)

Using the convolution theorem for FFT the energy (4.30) can be rewritten as

$$V_{\rm qr} = {1 \over 2 } \sum_{m_{1}=1}^{K_{1}}\sum_{m_{2}=1}^{K_{2}}\sum_{m_... ... Q \right ] (m_{1},m_{2},m_{3}) {\cal F} \left [ Q \right ] (m_{1},m_{2},m_{3})$$ (4.34)

We now use the identity $\sum_{\bf m} F(A)({\bf m}) B({\bf m}) = \sum_{\bf m} A({\bf m}) F(B)({\bf m})$ to arrive at

$$V_{\rm qr} = {1 \over 2 } \sum_{m_{1}=1}^{K_{1}}\sum_{m_{2}=1}^{K_{2}}\sum_{m_{3}=1}^{K_{3}} (\Theta_{rec} \star Q) ( m_{1},m_{2},m_{3}) Q (m_{1},m_{2},m_{3})$$

We first notice that $\Theta_{rec}$ does not depend on the charge positions and that $M_n \left ( u_{i\alpha} - k \right )$ is differentiable for $n > 2$ (which is always the case in practical applications). Thus the force on each charge can be obtained by taking the derivative of Eq. (4.35), namely

$$F_{i\alpha}^{({\rm qr})} = - { \partial V_{\rm qr} \over \partial... ...m_{3}) \over \partial r_{i\alpha}} (\Theta_{rec} \star Q )( m_{1},m_{2},m_{3}).$$ (4.35)

In practice, the calculation is carried out according to the following scheme: i) At each simulation step one computes the grid scaled fractional coordinates $u_{i\alpha}$ and fills an array with $Q$ according to Eq. (4.28). At this stage, the derivative of the $M_{n}$ functions are also computed and stored in memory. ii) The array containing $Q$ is then overwritten by ${\cal F} \left [ Q \right ]$, i.e. $Q$'s 3-D Fourier transform. iii) Subsequently, the electrostatic energy is computed via Eq. (4.30). At the same time, the array containing ${\cal F} \left [ Q \right ]$ is overwritten by the product of itself with the array containing $BC$ (computed at the very beginning of the run). iv) The resulting array is then Fourier transformed to obtain the convolution $\Theta_{rec} \star Q$. v) Finally, the forces are computed via Eq. (4.36) using the previously stored derivatives of the $M_{n}$ functions to recast $\partial Q /\partial r_{i\alpha}$.

The memory requirements of the SPME method are limited. $2K_{1}K_{2}K_{3}$ double precision real numbers are needed for the grid charge array $Q$, while the calculation of the functions $M_{n}(u_{i\alpha}-j)$ and their derivatives requires only $6 \times n \times N$ double precision real numbers. The $K_{\alpha}$ integers determines the fineness of the grid along the $\alpha$-th lattice vector of the unit cell. The output accuracy of the energy and forces depends on the SPME parameters: The $\alpha$ convergence parameter, the grid spacing and the order $n$ of the B-spline interpolation. For a typical $\alpha \simeq 0.4$ Å $^{-1}$ a relative accuracy between $10^{-4} - 10^{-5}$ for the electrostatic energy are obtained when the grid spacing is around 1 Å along each axis, and the order $n$ of the $B$-spline interpolation is 4 or 5. A rigorous error analysis and a comparison with standard Ewald summation can be found in Refs. [35] and [102]. For further readings on the PME and SPME techniques we refer to the original papers [34,102,29,35].

Figure: CPU time versus number of particles for the SPME algorithm as measured on a 43P/160MH IBM workstation

The power of the SPME algorithm, compared to the straightforward implementation of the standard Ewald method, is indeed astonishing. In Fig. (4.3) we report CPU timing obtained on a low end 43P/160MH IBM workstation for the evaluation of the reciprocal lattice energy and forces via SPME for cyanobiphenil as a function of the number of atoms in the system. Public domain 3-D FFT routines were used. The algorithm is practically linear and for 12,000 particles SPME takes only 2 CPU seconds to perform the calculation. A standard Ewald simulation for a box $64\times 64 \times 64$ Å$^{3}$ (i.e. with a grid spacing in $k$ space of $k = 2 \pi /64 \simeq 0.01 ~{\rm\AA}^{-1}$) for the same sample and at the same level of accuracy would have taken several minutes.