--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/20 15:43:13 2641 +++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/20 17:32:33 2643 @@ -65,54 +65,54 @@ In molecular simulations, proper accumulation of the e \section{Introduction} In molecular simulations, proper accumulation of the electrostatic -interactions is considered one of the most essential and -computationally demanding tasks. The common molecular mechanics force -fields are founded on representation of the atomic sites centered on -full or partial charges shielded by Lennard-Jones type interactions. -This means that nearly every pair interaction involves an -charge-charge calculation. Coupled with $r^{-1}$ decay, the monopole -interactions quickly become a burden for molecular systems of all -sizes. For example, in small systems, the electrostatic pair -interaction may not have decayed appreciably within the box length -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 -very time-consuming to perform these calculations. +interactions is essential and is one of the most +computationally-demanding tasks. The common molecular mechanics force +fields represent atomic sites with full or partial charges protected +by Lennard-Jones (short range) interactions. This means that nearly +every pair interaction involves a calculation of charge-charge forces. +Coupled with relatively long-ranged $r^{-1}$ decay, the monopole +interactions quickly become the most expensive part of molecular +simulations. Historically, the electrostatic pair interaction would +not have decayed appreciably within the typical box lengths that could +be feasibly simulated. In the larger systems that are more typical of +modern simulations, large cutoffs should be used to incorporate +electrostatics correctly. -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 +There have been many efforts to address the proper and practical +handling of electrostatic interactions, and these have resulted in 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} +often preferred because they physically incorporate 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 the computational cost, there have been some questions +regarding possible artifacts caused by the inherent periodicity of the +explicit Ewald summation.\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 +In this paper, we focus on a new set of shifted methods devised by +Wolf {\it et al.},\cite{Wolf99} which we further extend. These +methods along with a few other mixed methods (i.e. reaction field) are +compared with the smooth particle mesh Ewald +sum,\cite{Onsager36,Essmann99} which is our reference method for +handling long-range electrostatic interactions. The new methods for +handling electrostatics have the potential to scale linearly with +increasing system size since they involve only a simple modification +to the direct pairwise sum. They also lack the added periodicity of +the Ewald sum, so they can be used for systems which are non-periodic +or which have one- or two-dimensional periodicity. Below, these +methods are evaluated using a variety of model systems 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, +effect of all charges within a (cubic) 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} @@ -123,12 +123,13 @@ the cell length, $\bm{\Omega}_{i,j}$ are the Euler ang 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. +$j$, and $\phi$ is the solution to 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 divergent for +non-neutral systems. -This electrostatic summation problem was originally studied by Ewald +The electrostatic summation problem was originally studied by Ewald for the case of an infinite crystal.\cite{Ewald21}. The approach he took was to convert this conditionally convergent sum into two absolutely convergent summations: a short-ranged real-space summation @@ -140,23 +141,27 @@ where $\alpha$ is a damping parameter, or separation c \label{eq:EwaldSum} \end{equation} where $\alpha$ is a damping parameter, or separation constant, with -units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and equal -$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric -constant of the encompassing medium. The final two terms of +units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are +equal to $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the +dielectric constant of the surrounding medium. The final two terms of eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term for interacting with a surrounding dielectric.\cite{Allen87} This dipolar term was neglected in early applications in molecular simulations,\cite{Brush66,Woodcock71} until it was introduced by de Leeuw {\it et al.} to address situations where the unit cell has a -dipole moment and this dipole moment gets magnified through -replication of the periodic images.\cite{deLeeuw80,Smith81} If this -term is taken to be zero, the system is using conducting boundary +dipole moment which is magnified through replication of the periodic +images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the +system is said to be using conducting (or ``tin-foil'') boundary conditions, $\epsilon_{\rm S} = \infty$. Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied over -time. Initially, due to the small size of systems, the entire -simulation box was replicated to convergence. Currently, we balance a -spherical real-space cutoff with the reciprocal sum and consider the -surrounding dielectric. +time. Initially, due to the small sizes of the systems that could be +feasibly simulated, the entire simulation box was replicated to +convergence. In more modern simulations, the simulation boxes have +grown large enough that a real-space cutoff could potentially give +convergent behavior. Indeed, it has often been observed that the +reciprocal-space portion of the Ewald sum can be vanishingly +small compared to the real-space portion.\cite{XXX} + \begin{figure} \centering \includegraphics[width = \linewidth]{./ewaldProgression.pdf} @@ -170,65 +175,63 @@ a surrounding dielectric term is included.} \label{fig:ewaldTime} \end{figure} -The Ewald summation in the straight-forward form is an -$\mathscr{O}(N^2)$ algorithm. The separation constant $(\alpha)$ -plays an important role in the computational cost balance between the -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)$ 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 -reduce the cost of the Ewald summation further is to set $\alpha$ such -that the real-space interactions decay rapidly, allowing for a short -spherical cutoff, and then optimize the reciprocal space summation. -These optimizations usually involve the utilization of the fast -Fourier transform (FFT),\cite{Hockney81} leading to the +The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The +separation constant $(\alpha)$ plays an important role in balancing +the computational cost between the direct and reciprocal-space +portions of the summation. The choice of this value allows one to +select whether the real-space or reciprocal space portion of the +summation is an $\mathscr{O}(N^2)$ calculation (with the other being +$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of +$\alpha$ and thoughtful algorithm development, this cost can be +reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route +taken to reduce the cost of the Ewald summation even further is to set +$\alpha$ such that the real-space interactions decay rapidly, allowing +for a short spherical cutoff. Then the reciprocal space summation is +optimized. These optimizations usually involve utilization of the +fast Fourier transform (FFT),\cite{Hockney81} leading to the particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these methods, the cost of the reciprocal-space portion of the Ewald -summation is from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$. +summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N +\log N)$. -These developments and optimizations have led the use of the Ewald -summation to become routine in simulations with periodic boundary -conditions. However, in certain systems the intrinsic three -dimensional periodicity can prove to be problematic, such as two -dimensional surfaces and membranes. The Ewald sum has been -reformulated to handle 2D -systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new -methods have been found to be computationally -expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the -full Ewald summation is a possible direction for enabling the handling -of 2D systems and the inclusion of the optimizations described -previously.\cite{Yeh99} +These developments and optimizations have made the use of the Ewald +summation routine in simulations with periodic boundary +conditions. However, in certain systems, such as vapor-liquid +interfaces and membranes, the intrinsic three-dimensional periodicity +can prove problematic. The Ewald sum has been reformulated to handle +2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the +new methods are computationally expensive.\cite{Spohr97,Yeh99} +Inclusion of a correction term in the Ewald summation is a possible +direction for handling 2D systems while still enabling the use of the +modern optimizations.\cite{Yeh99} Several studies have recognized that the inherent periodicity in the -Ewald sum can also have an effect on systems that have the same -dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} -Good examples are solvated proteins kept at high relative -concentration due to the periodicity of the electrostatics. In these +Ewald sum can also have an effect on three-dimensional +systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} +Solvated proteins are essentially kept at high concentration due to +the periodicity of the electrostatic summation method. In these systems, the more compact folded states of a protein can be 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 -periodicity may negatively affect the system dynamics. +Ewald summation.\cite{Weber00} Thus, care must be taken when +considering the use of the Ewald summation where the assumed +periodicity would introduce spurious effects in the system dynamics. - \subsection{The Wolf and Zahn Methods} In a recent paper by Wolf \textit{et al.}, a procedure was outlined for the accurate accumulation of electrostatic interactions in an -efficient pairwise fashion and lacks the inherent periodicity of the -Ewald summation.\cite{Wolf99} Wolf \textit{et al.} observed that the -electrostatic interaction is effectively short-ranged in condensed -phase systems and that neutralization of the charge contained within -the cutoff radius is crucial for potential stability. They devised a -pairwise summation method that ensures charge neutrality and gives -results similar to those obtained with the Ewald summation. The -resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes -image-charges subtracted out through placement on the cutoff sphere -and a distance-dependent damping function (identical to that seen in -the real-space portion of the Ewald sum) to aid convergence +efficient pairwise fashion. This procedure lacks the inherent +periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} +observed that the electrostatic interaction is effectively +short-ranged in condensed phase systems and that neutralization of the +charge contained within the cutoff radius is crucial for potential +stability. They devised a pairwise summation method that ensures +charge neutrality and gives results similar to those obtained with the +Ewald summation. The resulting shifted Coulomb potential +(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through +placement on the cutoff sphere 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_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} @@ -255,22 +258,20 @@ the potential are not commensurate. Attempts to use b force expressions for use in simulations involving water.\cite{Zahn02} In their work, they pointed out that the forces and derivative of the potential are not commensurate. Attempts to use both -Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead +eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead to poor energy conservation. They correctly observed that taking the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the derivatives gives forces for a different potential energy function -than the one shown in Eq. (\ref{eq:WolfPot}). +than the one shown in eq. (\ref{eq:WolfPot}). -Zahn \textit{et al.} proposed a modified form of this ``Wolf summation -method'' as a way to use this technique in Molecular Dynamics -simulations. Taking the integral of the forces shown in equation -\ref{eq:WolfForces}, they proposed a new damped Coulomb -potential, +Zahn \textit{et al.} introduced a modified form of this summation +method as a way to use the technique in Molecular Dynamics +simulations. They proposed a new damped Coulomb potential, \begin{equation} -V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. +V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}, \label{eq:ZahnPot} \end{equation} -They showed that this potential does fairly well at capturing the +and showed that this potential does fairly well at capturing the structural and dynamic properties of water compared the same properties obtained using the Ewald sum. @@ -301,7 +302,7 @@ shifted potential, \textit{et al.} and Zahn \textit{et al.} by considering the standard shifted potential, \begin{equation} -v_\textrm{SP}(r) = \begin{cases} +V_\textrm{SP}(r) = \begin{cases} v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} \end{cases}, @@ -309,7 +310,7 @@ and shifted force, \end{equation} and shifted force, \begin{equation} -v_\textrm{SF}(r) = \begin{cases} +V_\textrm{SF}(r) = \begin{cases} v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} \end{cases}, @@ -325,16 +326,16 @@ of the unshifted potential itself (when inside the cut The forces associated with the shifted potential are simply the forces of the unshifted potential itself (when inside the cutoff sphere), \begin{equation} -f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), +F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), \end{equation} and are zero outside. Inside the cutoff sphere, the forces associated with the shifted force form can be written, \begin{equation} -f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d +F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d v(r)}{dr} \right)_{r=R_\textrm{c}}. \end{equation} -If the potential ($v(r)$) is taken to be the normal Coulomb potential, +If the potential, $v(r)$, is taken to be the normal Coulomb potential, \begin{equation} v(r) = \frac{q_i q_j}{r}, \label{eq:Coulomb} @@ -342,14 +343,14 @@ al.}'s undamped prescription: then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et al.}'s undamped prescription: \begin{equation} -v_\textrm{SP}(r) = +V_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, \label{eq:SPPot} \end{equation} with associated forces, \begin{equation} -f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. +F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. \label{eq:SPForces} \end{equation} These forces are identical to the forces of the standard Coulomb @@ -364,18 +365,18 @@ will give, The shifted force ({\sc sf}) form using the normal Coulomb potential will give, \begin{equation} -v_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}. +V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}. \label{eq:SFPot} \end{equation} with associated forces, \begin{equation} -f_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}. +F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}. \label{eq:SFForces} \end{equation} This formulation has the benefits that there are no discontinuities at -the cutoff distance, while the neutralizing image charges are present -in both the energy and force expressions. It would be simple to add -the self-neutralizing term back when computing the total energy of the +the cutoff radius, while the neutralizing image charges are present in +both the energy and force expressions. It would be simple to add the +self-neutralizing term back when computing the total energy of the system, thereby maintaining the agreement with the Madelung energies. A side effect of this treatment is the alteration in the shape of the potential that comes from the derivative term. Thus, a degree of @@ -383,13 +384,13 @@ Wolf \textit{et al.} originally discussed the energeti to gain functionality in dynamics simulations. Wolf \textit{et al.} originally discussed the energetics of the -shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that -it was still insufficient for accurate determination of the energy -with reasonable cutoff distances. The calculated Madelung energies -fluctuate around the expected value with increasing cutoff radius, but -the oscillations converge toward the correct value.\cite{Wolf99} A +shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was +insufficient for accurate determination of the energy with reasonable +cutoff distances. The calculated Madelung energies fluctuated around +the expected value as the cutoff radius was increased, but the +oscillations converged toward the correct value.\cite{Wolf99} A damping function was incorporated to accelerate the convergence; and -though alternative functional forms could be +though alternative forms for the damping function could be used,\cite{Jones56,Heyes81} the complimentary error function was chosen to mirror the effective screening used in the Ewald summation. Incorporating this error function damping into the simple Coulomb @@ -398,37 +399,37 @@ v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, \label{eq:dampCoulomb} \end{equation} -the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using -eq. (\ref{eq:shiftingForm}), +the shifted potential (eq. (\ref{eq:SPPot})) becomes \begin{equation} -v_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, +V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, \label{eq:DSPPot} \end{equation} with associated forces, \begin{equation} -f_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}. +F_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}. \label{eq:DSPForces} \end{equation} -Again, this damped shifted potential suffers from a discontinuity and -a lack of the image charges in the forces. To remedy these concerns, -one may derive a {\sc sf} variant by including the derivative -term in eq. (\ref{eq:shiftingForm}), +Again, this damped shifted potential suffers from a +force-discontinuity at the cutoff radius, and the image charges play +no role in the forces. To remedy these concerns, one may derive a +{\sc sf} variant by including the derivative term in +eq. (\ref{eq:shiftingForm}), \begin{equation} \begin{split} -v_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}. +V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}. \label{eq:DSFPot} \end{split} \end{equation} The derivative of the above potential will lead to the following forces, \begin{equation} \begin{split} -f_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}. +F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}. \label{eq:DSFForces} \end{split} \end{equation} -If the damping parameter $(\alpha)$ is chosen to be zero, the undamped -case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered -from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}). +If the damping parameter $(\alpha)$ is set to zero, the undamped case, +eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly +recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important @@ -440,55 +441,48 @@ would be expected to have sudden jumps as particle dis portion is different. The missing $v_\textrm{c}$ term would not affect molecular dynamics simulations (although the computed energy would be expected to have sudden jumps as particle distances crossed -$R_c$). The sign problem would be a potential source of errors, -however. In fact, it introduces a discontinuity in the forces at the -cutoff, because the force function is shifted in the wrong direction -and doesn't cross zero at $R_\textrm{c}$. +$R_c$). The sign problem is a potential source of errors, however. +In fact, it introduces a discontinuity in the forces at the cutoff, +because the force function is shifted in the wrong direction and +doesn't cross zero at $R_\textrm{c}$. Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an -electrostatic summation method that is continuous in both the -potential and forces and which incorporates the damping function -proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of this -paper, we will evaluate exactly how good these methods ({\sc sp}, {\sc -sf}, damping) are at reproducing the correct electrostatic summation -performed by the Ewald sum. +electrostatic summation method in which the potential and forces are +continuous at the cutoff radius and which incorporates the damping +function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of +this paper, we will evaluate exactly how good these methods ({\sc sp}, +{\sc sf}, damping) are at reproducing the correct electrostatic +summation performed by the Ewald sum. \subsection{Other alternatives} -In addition to the methods described above, we will consider some -other techniques that commonly get used in molecular simulations. The +In addition to the methods described above, we considered some other +techniques that are commonly used in molecular simulations. The simplest of these is group-based cutoffs. Though of little use for -non-neutral molecules, collecting atoms into neutral groups takes +charged molecules, collecting atoms into neutral groups takes advantage of the observation that the electrostatic interactions decay faster than those for monopolar pairs.\cite{Steinbach94} When -considering these molecules as groups, an orientational aspect is -introduced to the interactions. Consequently, as these molecular -particles move through $R_\textrm{c}$, the energy will drift upward -due to the anisotropy of the net molecular dipole -interactions.\cite{Rahman71} To maintain good energy conservation, -both the potential and derivative need to be smoothly switched to zero -at $R_\textrm{c}$.\cite{Adams79} This is accomplished using a -switching function, -\begin{equation} -S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\ -\frac{(R_\textrm{c}+2r-3r_\textrm{sw})(R_\textrm{c}-r)^2}{(R_\textrm{c}-r_\textrm{sw})^3} &\quad r_\textrm{sw}R_\textrm{c} -\end{cases}, -\end{equation} -where the above form is for a cubic function. If a smooth second -derivative is desired, a fifth (or higher) order polynomial can be -used.\cite{Andrea83} +considering these molecules as neutral groups, the relative +orientations of the molecules control the strength of the interactions +at the cutoff radius. Consequently, as these molecular particles move +through $R_\textrm{c}$, the energy will drift upward due to the +anisotropy of the net molecular dipole interactions.\cite{Rahman71} To +maintain good energy conservation, both the potential and derivative +need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79} +This is accomplished using a standard switching function. If a smooth +second derivative is desired, a fifth (or higher) order polynomial can +be used.\cite{Andrea83} Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, -and to incorporate their effect, a method like Reaction Field ({\sc -rf}) can be used. The original theory for {\sc rf} was originally -developed by Onsager,\cite{Onsager36} and it was applied in -simulations for the study of water by Barker and Watts.\cite{Barker73} -In application, it is simply an extension of the group-based cutoff -method where the net dipole within the cutoff sphere polarizes an -external dielectric, which reacts back on the central dipole. The -same switching function considerations for group-based cutoffs need to -made for {\sc rf}, with the additional pre-specification of a -dielectric constant. +and to incorporate the effects of the surroundings, a method like +Reaction Field ({\sc rf}) can be used. The original theory for {\sc +rf} was originally developed by Onsager,\cite{Onsager36} and it was +applied in simulations for the study of water by Barker and +Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply +an extension of the group-based cutoff method where the net dipole +within the cutoff sphere polarizes an external dielectric, which +reacts back on the central dipole. The same switching function +considerations for group-based cutoffs need to made for {\sc rf}, with +the additional pre-specification of a dielectric constant. \section{Methods}