Electrostatic Corrections for the Multiple Time Step Simulation

In flexible molecular systems of large size, the Ewald summation presents computational problems which are crucial to constructing efficient and stable multiple time step integrators [103,2]. We have seen that intra-molecular Coulomb interactions between bonded atoms or between atoms bonded to a common atom are excluded in most of the standard force fields for protein simulation. In any practical implementation of the Ewald method, the intra-molecular energy $V_{\rm intra}$ is automatically included in the reciprocal space summation and is subtracted in direct space (see Eqs. (4.23,4.21). In actual simulations the reciprocal space sum is computed with a finite accuracy whereas the intra-molecular term $V_{\rm intra}$, due to the excluded Coulomb interactions, is computed exactly. This clearly prevents the cancellation of the intra-molecular forces and energies. When the stretching and bending forces are integrated explicitly, the intra-molecular term due to the excluded contacts varies rapidly with time and so does the cancellation error. Consequently, instability may be observed when integrating the reciprocal lattice forces in reference systems with large time steps. The correction due to the truncation can be evaluated by approximating the reciprocal lattice sum for the excluded contacts in Eq. (4.21) to an integral in the 3-dimensional $k$ space and evaluating this integral from the cutoff $k_{cut} \equiv 2 \pi \vert{\bf m}\vert _{max}$ to infinity in polar coordinates. The neglected reciprocal lattice intra-molecular energy is then [104]

$$V_{\rm corr} = {\alpha \over \pi^{1/2}} {\rm erfc}(k_{cut}/2 \alp... ...q_{i}^{2} + \sum_{ij-\textrm{excl.}} q_{i}q_{j} \chi(r_{ij},k_{\rm cut},\alpha)$$ (4.44)

with

$$\chi(r,k_{\rm cut},\alpha) = {2 \over \pi } \int_{k_{\rm cut}}^{\infty} e^{-k^{2}/4\alpha^{2}} {\sin(kr) \over kr} dk .$$ (4.45)

The first constant term in (4.46) refers to the self energy, while the second accounts for the intra-molecular excluded interactions^4.4. This correction must be included in the same reference systems to which $V_{\rm qr}$ is assigned, e.g. $V_{l}$ in our potential separation (see Table 2).

In principle the correction in Eq. (4.46) applies only to standard Ewald and not to the reciprocal lattice energy computed via SPME. We can still, however, use the correction Eq. (4.46), if a spherical cutoff $k_{\rm cut}$ is applied to SPME. This can be done easily by setting $\exp(-\pi^2{\bf m}^2/\alpha^{2})/{\bf m}^{2}=0$ for $2 \pi {\bf m} > k_{\rm cut} \equiv f_{f} \pi N_{f}/L$ where $L$ is the side length of the cubic box and $N_{f}$ is the number of grid points in each direction. The factor $f_{f}$ must be chosen slightly less than unity. This simple device decreases the effective cutoff in reciprocal space while maintaining the same grid spacing, thus reducing the B-spline interpolation error (the error in the B-spline interpolation of the complex exponential is, indeed, maximum precisely at the tail of the reciprocal sums [35]). In Ref. [104] the effect of including or not such correction in electrostatic systems using multiple time step algorithms is studied and discussed thoroughly.

Figure: The correction potential $\chi (r,k_{c},\alpha )$ as a function of the distance for different values of the parameters $\alpha$ (left) and $k_{c}$ (right). The solid line on the top right corner is the bare Coulomb potential $1/r$

.

The potential $\chi(r,k_{\rm cut},\alpha)$ yields, in direct space, the neglected reciprocal energy due to the truncation of the reciprocal lattice sums, and must, in principle, be included for each atom pair distance in direct space. Thus, the corrected direct space potential is then

$$V_{\rm qd}' = {1 \over 2} \sum_{ij}^{N} q_i q_j \left [ \sum_n { {\rm erfc} \al... ...}_n \vert } + \chi(\vert r_{ij} + {\bf r}_{n}\vert,k_{\rm cut},\alpha) \right ]$$ (4.46)

which is then split as usual in short-medium-long range according to (4.39). The correction is certainly more crucial for the excluded intramolecular contacts because $V_{\rm corr}$ is essentially a short-ranged potential which is non negligible only for intramolecular short distances. For systems with hydrogen bonds, however, the correction is also important for intermolecular interactions.

In Fig. 4.4 the correction potential is compared to the Coulomb potential (solid line in the top right corner) for different value of the reciprocal space cutoff $k_{c}$ and of the convergence parameter $\alpha$. For practical values of $\alpha$ and $k_{c}$, the potential is short ranged and small compared to the bare $1/r$ Coulomb interaction. In the asymptotic limit $V_{corr}$ goes to zero as ${\rm sin}\left ( a r \right )/r^{2}$ where $a$ is a constant. This oscillatory long range behavior of the correction potential $V_{\rm corr}$ is somewhat nasty: In Fig. 4.5 we show the integral

$$I (r,k_{c},\alpha) = \int_{0}^{r} \chi(x,k_{c},\alpha) x^{2} dx$$ (4.47)

as a function of the distance. If this integral converges then the $\chi(r,k)$ is absolutely convergent in 3D.

Figure: The integral $I(r)$ of Eq. (4.49) as a function of the distance for different values of the parameters $\alpha$ (left) and $k_{c}$ (right)

We see that the period of the oscillations in $I(r)$ increases with $k_{c}$ while $\alpha$ affects only the amplitude. The total energy is hence again conditionally convergent, since the limit $lim_{r\rightarrow \infty} I(r)$ does not exist. However, unlike for the $1/r$ bare potential, the energy integral remains in this case bounded. Due to this, a cutoff on the small potential $V_{\rm corr}$ is certainly far less dangerous that a cutoff on the bare $1/r$ term. In order to verify this, we have calculated some properties of liquid water using the SPC model[105] from a 200 ps MD simulation in the $NPT$ ensemble at temperature of 300 K and pressure of 0.1 MPa with i) a very accurate Ewald sum (column EWALD in Table 4.2), ii) with inaccurate Ewald but corrected in direct space using Eq. (4.48) (CORRECTED) and iii) with simple cutoff truncation of the bare Coulomb potential and no Ewald (CUTOFF). Results are reported in Table 4.2

Table: Properties of liquid water computed from a 200 ps simulation at 300 K and 0.1 Mpa on a sample of 343 molecules in PBC with accurate Ewald ( $\alpha = 0.35 ~\AA^{-1}$; $k_{c} = 2.8 ~\AA^{-1}$ ) and no correction Eq. (4.46) (column EWALD), with inaccurate Ewald ( $\alpha = 0.35 ~\AA^{-1}$; $k_{c} = 0.9 ~\AA^{-1}$) but including the correction Eq. (4.46) and with no Ewald and cutoff at 10.0 Å. $R_{0-0}$ is the distance corresponding to the first peak in the Oxygen-Oxygen pair distribution function.

	EWALD	CORRECTED	CUTOFF
Coulomb energy (KJ/mole)	-55.2 $\pm 0.1$	-55.1 $\pm 0.1$	-56.4 $\pm 0.1$
Potential energy (KJ/mole)	-46.2 $\pm 0.1$	-46.1 $\pm 0.1$	-47.3 $\pm 0.1$
Heat Capacity (KJ/mole/K)	74 $\pm 24.5$	94 $\pm 22.0$	87 $\pm 23.2$
Volume (cm$^{3}$)	18.2 $\pm 0.1$	18.3 $\pm 0.1$	18.1 $\pm 0.1$
Volume Fluctuation (Å$^{3}$)	136.9$\pm 3.5$	147.0 $\pm 3.5$	138.7$\pm 3.5$
$R_{0-0}$ (Å)	2.81 $\pm 0.01$	2.81 $\pm 0.01$	2.81 $\pm 0.01$
Dielectric constant	59 $\pm 25.8$	47 $\pm 27.3$	3 $\pm 2$

We notice that almost all the computed properties of water are essentially independent, within statistical error, of the truncation method. The dielectric properties, on the contrary, appear very sensitive to the method for dealing with long range tails: Accurate and inaccurate Ewald (corrected in direct space through 4.48) yields, within statistical error, comparable results whereas the dielectric constant predicted by the spherical cutoff method is more than order of magnitude smaller. We should remark that method ii) (CORRECTED) is almost twice as efficient as the ``exact'' method i).