--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/19 19:34:53 2638 +++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/20 15:43:13 2641 @@ -77,25 +77,56 @@ accurately incorporate their effect, and since the com leading to an effect excluded from the pair interactions within a unit box. In large systems, excessively large cutoffs need to be used to accurately incorporate their effect, and since the computational cost -increases proportionally with the cutoff sphere, it quickly becomes an -impractical task to perform these calculations. +increases proportionally with the cutoff sphere, it quickly becomes +very time-consuming to perform these calculations. +There have been many efforts to address this issue of both proper and +practical handling of electrostatic interactions, and these have +resulted in the availability of a variety of +techniques.\cite{Roux99,Sagui99,Tobias01} These are typically +classified as implicit methods (i.e., continuum dielectrics, static +dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e., +Ewald summations, interaction shifting or +truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., +reaction field type methods, fast multipole +methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are +often preferred because they incorporate dynamic solvent molecules in +the system of interest, but these methods are sometimes difficult to +utilize because of their high computational cost.\cite{Roux99} In +addition to this cost, there has been some question of the inherent +periodicity of the explicit Ewald summation artificially influencing +systems dynamics.\cite{Tobias01} + +In this paper, we focus on the common mixed and explicit methods of +reaction filed and smooth particle mesh +Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods +devised by Wolf {\it et al.} which we further extend.\cite{Wolf99} +These new methods for handling electrostatics are quite +computationally efficient, since they involve only a simple +modification to the direct pairwise sum, and they lack the added +periodicity of the Ewald sum. Below, these methods are evaluated using +a variety of model systems and comparison methodologies to establish +their usability in molecular simulations. + \subsection{The Ewald Sum} -The complete accumulation electrostatic interactions in a system with periodic boundary conditions (PBC) requires the consideration of the effect of all charges within a simulation box, as well as those in the periodic replicas, +The complete accumulation electrostatic interactions in a system with +periodic boundary conditions (PBC) requires the consideration of the +effect of all charges within a simulation box, as well as those in the +periodic replicas, \begin{equation} V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime \left[ \sum_{i=1}^N\sum_{j=1}^N \phi\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) \right], \label{eq:PBCSum} \end{equation} -where the sum over $\mathbf{n}$ is a sum over all periodic box replicas -with integer coordinates $\mathbf{n} = (l,m,n)$, and the prime indicates -$i = j$ are neglected for $\mathbf{n} = 0$.\cite{deLeeuw80} Within the -sum, $N$ is the number of electrostatic particles, $\mathbf{r}_{ij}$ is -$\mathbf{r}_j - \mathbf{r}_i$, $L$ is the cell length, $\bm{\Omega}_{i,j}$ are -the Euler angles for $i$ and $j$, and $\phi$ is Poisson's equation -($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge -interactions). In the case of monopole electrostatics, -eq. (\ref{eq:PBCSum}) is conditionally convergent and is discontiuous -for non-neutral systems. +where the sum over $\mathbf{n}$ is a sum over all periodic box +replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the +prime indicates $i = j$ are neglected for $\mathbf{n} = +0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic +particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is +the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and +$j$, and $\phi$ is Poisson's equation ($\phi(\mathbf{r}_{ij}) = q_i +q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge interactions). In the +case of monopole electrostatics, eq. (\ref{eq:PBCSum}) is +conditionally convergent and is discontinuous for non-neutral systems. This electrostatic summation problem was originally studied by Ewald for the case of an infinite crystal.\cite{Ewald21}. The approach he @@ -145,7 +176,7 @@ real-space or reciprocal space portion of the summatio direct and reciprocal-space portions of the summation. The choice of the magnitude of this value allows one to select whether the real-space or reciprocal space portion of the summation is an -$\mathscr{O}(N^2)$ calcualtion (with the other being +$\mathscr{O}(N^2)$ calculation (with the other being $\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$ and thoughtful algorithm development, this cost can be brought down to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to @@ -181,7 +212,7 @@ considering the use of the Ewald summation where the i artificially stabilized by the periodic replicas introduced by the Ewald summation.\cite{Weber00} Thus, care ought to be taken when considering the use of the Ewald summation where the intrinsic -perodicity may negatively affect the system dynamics. +periodicity may negatively affect the system dynamics. \subsection{The Wolf and Zahn Methods} @@ -199,7 +230,7 @@ the real-space portion of the Ewald sum) to aid conver and a distance-dependent damping function (identical to that seen in the real-space portion of the Ewald sum) to aid convergence \begin{equation} -V_{\textrm{Wolf}}(r_{ij})= \frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. +V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. \label{eq:WolfPot} \end{equation} Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted @@ -542,7 +573,7 @@ between those computed from the particular method and investigated through measurement of the angle ($\theta$) formed between those computed from the particular method and those from SPME, \begin{equation} -\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, +\theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right), \end{equation} where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force vector computed using method $M$. @@ -630,17 +661,17 @@ snapshots were taken at regular intervals from higher Generation of the system configurations was dependent on the system type. For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 -SPC/E water molecule trajectories and each equilibrated individually. -The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- -ions and were selected and equilibrated in the same fashion as the -water systems. For the low and high ionic strength NaCl solutions, 4 -and 40 ions were first solvated in a 1000 water molecule boxes -respectively. Ion and water positions were then randomly swapped, and -the resulting configurations were again equilibrated individually. -Finally, for the Argon/Water "charge void" systems, the identities of -all the SPC/E waters within 6 \AA\ of the center of the equilibrated -water configurations were converted to argon -(Fig. \ref{fig:argonSlice}). +SPC/E water molecule trajectories and each equilibrated +individually.\cite{Berendsen87} The solid and liquid NaCl systems +consisted of 500 Na+ and 500 Cl- ions and were selected and +equilibrated in the same fashion as the water systems. For the low +and high ionic strength NaCl solutions, 4 and 40 ions were first +solvated in a 1000 water molecule boxes respectively. Ion and water +positions were then randomly swapped, and the resulting configurations +were again equilibrated individually. Finally, for the Argon/Water +"charge void" systems, the identities of all the SPC/E waters within 6 +\AA\ of the center of the equilibrated water configurations were +converted to argon (Fig. \ref{fig:argonSlice}). \begin{figure} \centering @@ -695,19 +726,20 @@ realistic results using an unmodified cutoff. This is In this figure, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff. This is not all that -surprising since this results in large energy fluctuations as atoms -move in and out of the cutoff radius. These fluctuations can be -alleviated to some degree by using group based cutoffs with a -switching function.\cite{Steinbach94} The Group Switch Cutoff row -doesn't show a significant improvement in this plot because the salt -and salt solution systems contain non-neutral groups, see the +surprising since this results in large energy fluctuations as atoms or +molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79} +These fluctuations can be alleviated to some degree by using group +based cutoffs with a switching +function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff +row doesn't show a significant improvement in this plot because the +salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral. Correcting the resulting charged cutoff sphere is one of the purposes of the damped Coulomb summation proposed by Wolf \textit{et al.},\cite{Wolf99} and this correction indeed improves the results as -seen in the Shifted-Potental rows. While the undamped case of this +seen in the {\sc sp} rows. While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME. Inclusion of potential damping improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows @@ -920,7 +952,7 @@ the point charges for the pairwise summation methods; increased, these peaks are smoothed out, and approach the SPME curve. The damping acts as a distance dependent Gaussian screening of the point charges for the pairwise summation methods; thus, the -collisions are more elastic in the undamped {\sc sf} potental, and the +collisions are more elastic in the undamped {\sc sf} potential, and the stiffness of the potential is diminished as the electrostatic interactions are softened by the damping function. With $\alpha$ values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are