ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipole/multipole_2/multipole2.tex
(Generate patch)

Comparing trunk/multipole/multipole_2/multipole2.tex (file contents):
Revision 4207 by gezelter, Thu Aug 7 21:13:00 2014 UTC vs.
Revision 4208 by gezelter, Fri Aug 8 15:11:49 2014 UTC

# Line 73 | Line 73 | preprint,
73    energy differences between configurations, molecular forces, and
74    torques were used to analyze how well the real-space models perform
75    relative to the more computationally expensive Ewald treatment.  We
76 <  have also investigated the energy conservation properties of the new
77 <  methods in molecular dynamics simulations. The SP method shows
78 <  excellent agreement with configurational energy differences, forces,
79 <  and torques, and would be suitable for use in Monte Carlo
80 <  calculations.  Of the two new shifted-force methods, the GSF
81 <  approach shows the best agreement with Ewald-derived energies,
82 <  forces, and torques and also exhibits energy conservation properties
83 <  that make it an excellent choice for efficient computation of
84 <  electrostatic interactions in molecular dynamics simulations.
76 >  have also investigated the energy conservation, structural, and
77 >  dynamical properties of the new methods in molecular dynamics
78 >  simulations. The SP method shows excellent agreement with
79 >  configurational energy differences, forces, and torques, and would
80 >  be suitable for use in Monte Carlo calculations.  Of the two new
81 >  shifted-force methods, the GSF approach shows the best agreement
82 >  with Ewald-derived energies, forces, and torques and also exhibits
83 >  energy conservation properties that make it an excellent choice for
84 >  efficient computation of electrostatic interactions in molecular
85 >  dynamics simulations.  Both SP and GSF are able to reproduce
86 >  structural and dyanamical properties in the liquid models with
87 >  excellent fidelity.
88   \end{abstract}
89  
90   %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
# Line 185 | Line 188 | point charges.\cite{Fukuda:2013sf}
188      orientational averaging helps to reduce the effective range of the
189      interactions in the fluid.  Placement of reversed image multipoles
190      on the surface of the cutoff sphere recovers the effective
191 <    higher-order multipole behavior.}
189 <  \label{fig:schematic}
191 >    higher-order multipole behavior. \label{fig:schematic}}
192   \end{figure}
193  
194   One can make a similar effective range argument for crystals of point
# Line 232 | Line 234 | interacting molecules traverse each other's cutoff bou
234   and reduce the effective range of the multipolar interactions as
235   interacting molecules traverse each other's cutoff boundaries.
236  
235 – % Because of this reason, although the nature of electrostatic
236 – % interaction short ranged, the hard cutoff sphere creates very large
237 – % fluctuation in the electrostatic energy for the perfect crystal. In
238 – % addition, the charge neutralized potential proposed by Wolf et
239 – % al. converged to correct Madelung constant but still holds oscillation
240 – % in the energy about correct Madelung energy.\cite{Wolf:1999dn}.  This
241 – % oscillation in the energy around its fully converged value can be due
242 – % to the non-neutralized value of the higher order moments within the
243 – % cutoff sphere.
244 –
237   Forces and torques acting on atomic sites are fundamental in driving
238   dynamics in molecular simulations, and the damped shifted force (DSF)
239   energy kernel provides consistent energies and forces on charged atoms
# Line 304 | Line 296 | reference method, a full multipolar Ewald treatment.\c
296   reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98}
297  
298  
307 – %\subsection{Conservation of total energy }
308 – %To conserve the total energy in MD simulations, the energy, force, and torque on a central molecule due to another molecule should smoothly tends to zero as second molecule approaches to cutoff radius. In addition, the force should be derivable from the energy and vice versa. If only the first condition holds but not the second, the total energy does not conserve.\cite{Fennell:2006lq}. The hard cutoff method does not ensure the smooth transition of the energy, force, and torque at the cutoff radius.\cite{Wolf:1999dn} By placing image charge on the surface of the cutoff sphere, the smooth transition of the energy can be ensured but the force and torque remains discontinuous. Therefore the purposed methods should have smooth transition of the energy, force and torque to ensure the total energy conservation and the expression should be close to idea of placing image multipole on the surface of the cutoff sphere.
309 –
299   \section{\label{sec:method}Review of Methods}
300   Any real-space electrostatic method that is suitable for MD
301   simulations should have the electrostatic energy, forces and torques
302   between two sites go smoothly to zero as the distance between the
303 < sites, $r_{\bf ab}$ approaches the cutoff radius, $r_c$.  Requiring
303 > sites, $r_{ab}$ approaches the cutoff radius, $r_c$.  Requiring
304   this continuity at the cutoff is essential for energy conservation in
305   MD simulations.  The mathematical details of the shifted potential
306   (SP), gradient-shifted-force (GSF) and Taylor shifted-force (TSF)
# Line 330 | Line 319 | $a$, etc.
319   where the multipole operator for site $a$, $M_{a}$, is
320   expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object
321   $a$, etc.
333 –
334 – % Interactions between multipoles can be expressed as higher derivatives
335 – % of the bare Coulomb potential, so one way of ensuring that the forces
336 – % and torques vanish at the cutoff distance is to include a larger
337 – % number of terms in the truncated Taylor expansion, e.g.,
338 – % %
339 – % \begin{equation}
340 – % f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert  _{r_c}  .
341 – % \end{equation}
342 – % %
343 – % The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
344 – % Thus, for $f(r)=1/r$, we find
345 – % %
346 – % \begin{equation}
347 – % f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} .
348 – % \end{equation}
349 – % This function is an approximate electrostatic potential that has
350 – % vanishing second derivatives at the cutoff radius, making it suitable
351 – % for shifting the forces and torques of charge-dipole interactions.
322  
323   The TSF potential for any multipole-multipole interaction can be
324   written
# Line 363 | Line 333 | performed at cutoff radius ($r_c$) to obtain $f_n(r)$.
333   quadrupole-quadrupole.  To ensure smooth convergence of the energy,
334   force, and torques, a Taylor expansion with $n$ terms must be
335   performed at cutoff radius ($r_c$) to obtain $f_n(r)$.
366 –
367 – % To carry out the same procedure for a damped electrostatic kernel, we
368 – % replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$.
369 – % Many of the derivatives of the damped kernel are well known from
370 – % Smith's early work on multipoles for the Ewald
371 – % summation.\cite{Smith82,Smith98}
336  
373 – % Note that increasing the value of $n$ will add additional terms to the
374 – % electrostatic potential, e.g., $f_2(r)$ includes orders up to
375 – % $(r-r_c)^3/r_c^4$, and so on.  Successive derivatives of the $f_n(r)$
376 – % functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) =
377 – % f^{\prime\prime}_2(r)$, etc.  These higher derivatives are required
378 – % for computing multipole energies, forces, and torques, and smooth
379 – % cutoffs of these quantities can be guaranteed as long as the number of
380 – % terms in the Taylor series exceeds the derivative order required.
381 –
337   For multipole-multipole interactions, following this procedure results
338   in separate radial functions for each of the distinct orientational
339   contributions to the potential, and ensures that the forces and
# Line 527 | Line 482 | in the test cases are given in table~\ref{tab:pars}.
482   in the test cases are given in table~\ref{tab:pars}.
483  
484   \begin{table}
530 – \label{tab:pars}
485   \caption{The parameters used in the systems used to evaluate the new
486    real-space methods.  The most comprehensive test was a liquid
487    composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-}
488    ions).  This test excercises all orders of the multipolar
489 <  interactions developed in the first paper.}
489 >  interactions developed in the first paper.\label{tab:pars}}
490   \begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline
491               & \multicolumn{2}{c|}{LJ parameters} &
492               \multicolumn{5}{c|}{Electrostatic moments} & & & & \\
# Line 739 | Line 693 | model must allow for long simulation times with minima
693  
694   \section{\label{sec:result}RESULTS}
695   \subsection{Configurational energy differences}
742 – %The magnitude of the fluctuation in the total electrostatic energy per molecule for a dipolar crystal is very high as shown in (Fig … paper I).\cite{PaperI}As soon as, the net dipole moment within a cutoff radius is neutralized in the SP method, the magnitude of the fluctuation in the total electrostatic energy per molecule reduced significantly and rapidly converged to the correct energy constant (Refer figure … Paper I).\cite{PaperI} The GSF potential energy also converged to the correct energy constant for the cutoff radius rc = 6a for the undamped case. The potential energy from the TSF method converges towards the correct value for a very large cutoff radius. The speed of convergence for the all the cutoff methods can be increased by using damping function as shown in figure … Paper I\cite{PaperI}. For the quadrupolar crystals, the fluctuation in the total electrostatic energy for the hard cutoff method is small and short ranged as compared to the dipolar crystals (figure … in the paper I).\cite{PaperI}  Similar to the dipolar crystals, the net quadrupolar neutralization of the cutoff sphere in the SP method reduces oscillation rapidly and converge electrostatic energy to the correct energy constant.
743 – %The oscillation in the the electrostatic energy for the hard cutoff method is even true for the dipolar liquids as shown in Figure ~\ref{fig:rcutConvergence_dipolarLiquid}. As we placed image on the surface of the cutoff sphere in SP method, the oscillation in the energy is reduced. The fluctuation in the energy in liquid is much smaller as compared to the crystal (This result is similar to the results observed by \textit{Wolf et al.} in the case of ionic crystal and Mgo melt). The large magnitude in the fluctuation of the electrostatic energy in the crystal is because of large range of multipole ordering in the crystal. When the energy is evaluated by the direct truncation, it breaks up large number of multipolar ordering leaving behind net multipole moment. But in the case of liquid, there is only local ordering of the multipoles and their ordering disappears in the long range. Therefore, the direct truncation results a small oscillation in the electrostatic energy (which is smaller than deviation SP energy from the Ewald Figure ~\ref{fig:rcutConvergence_dipolarLiquid}) in the case of liquid. Although, the oscillation in the energy is very small for the case of liquid, this affects the change in potential energy ($\triangle E$), which is observed when $\triangle E$ evaluated from the SP method compared with Ewald as shown in figure 4a and 4b.
744 – %\begin{figure}[h!]
745 – %        \centering
746 – %        \includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf}
747 – %        \caption{The energy per molecule plotted against cutoff radius, rc for i) Hard ii) SP iii) GSF, and iv) TSF method. The hard cutoff method shows fluctuation in the electorstatic energy and it disappers in all other methods.  }
748 – %        \label{fig:rcutConvergence_dipolarLiquid}
749 – %    \end{figure}
750 – %In MC simulations, the electrostatic differences between the molecules are important parameter for sampling. We have compared $\triangle E$ from the different methods (Hard, SP, GSF, and TSF) with the Ewald using linear regression analysis. The correlation coefficient ($R^2$) of the regression line measures the deviation of the evaluated quantities from the mean slope. We know that Ewald method evaluates accurate value of the electrostatic energy. Hence, if the proposed methods can quantify the electrostatic energy as good as Ewald then the correlation coefficient is 1.The correlation coefficient is 0 for the completely random result for any physical quantity measured by the proposed method. The slope is a measure of the accuracy of the average of a physical quantity obtained from the proposed methods. If the slope is 1 then we can conclude that the average of the physical quantity measured by the method is as good as Ewald. The deviation of the slope from 1 state that the method used in quantifying physical quantities is statistically biased as compared to Ewald.
751 – %\begin{figure}
752 – %        \centering
753 – %        \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf}
754 – %        \label{fig:barGraph1}
755 – %        \end{figure}
756 – %        \begin{figure}
757 – %        \centering
758 – %       \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf}
759 – %        \caption{}
760 –      
761 – %        \label{fig:barGraph2}
762 – %      \end{figure}
763 – %The correlation coefficient ($R^2$) and slope of the linear
764 – %regression plots for the energy differences for all six different
765 – %molecular systems is shown in figure 4a and 4b.The plot shows that
766 – %the correlation coefficient improves for the SP cutoff method as
767 – %compared to the undamped hard cutoff method in the case of SSDQC,
768 – %SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar
769 – %crystal and liquid, the correlation coefficient is almost unchanged
770 – %and close to 1.  The correlation coefficient is smallest (0.696276
771 – %for $r_c$ = 9 $A^\circ$) for the SSDQC liquid because of the presence of
772 – %charge-charge and charge-multipole interactions. Since the
773 – %charge-charge and charge-multipole interaction is long ranged, there
774 – %is huge deviation of correlation coefficient from 1. Similarly, the
775 – %quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with
776 – %compared to interactions in the other multipolar systems, thus the
777 – %correlation coefficient very close to 1 even for hard cutoff
778 – %method. The idea of placing image multipole on the surface of the
779 – %cutoff sphere improves the correlation coefficient and makes it close
780 – %to 1 for all types of multipolar systems. Similarly the slope is
781 – %hugely deviated from the correct value for the lower order
782 – %multipole-multipole interaction and slightly deviated for higher
783 – %order multipole – multipole interaction. The SP method improves both
784 – %correlation coefficient ($R^2$) and slope significantly in SSDQC and
785 – %dipolar systems.  The Slope is found to be deviated more in dipolar
786 – %crystal as compared to liquid which is associated with the large
787 – %fluctuation in the electrostatic energy in crystal. The GSF also
788 – %produced better values of correlation coefficient and slope with the
789 – %proper selection of the damping alpha (Interested reader can consult
790 – %accompanying supporting material). The TSF method gives good value of
791 – %correlation coefficient for the dipolar crystal, dipolar liquid,
792 – %SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the
793 – %regression slopes are significantly deviated.
696  
697   \begin{figure}
698    \centering
# Line 802 | Line 704 | model must allow for long simulation times with minima
704      from those obtained using the multipolar Ewald sum.  Different
705      values of the cutoff radius are indicated with different symbols
706      (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted
707 <    triangles).}
806 <  \label{fig:slopeCorr_energy}
707 >    triangles).\label{fig:slopeCorr_energy}}
708   \end{figure}
709  
710   The combined correlation coefficient and slope for all six systems is
# Line 865 | Line 766 | forces is desired.
766   systematic error in the forces is concerning if replication of Ewald
767   forces is desired.
768  
769 + It is important to note that the forces and torques from the SP and
770 + the Hard cutoffs are not identical. The SP method shifts each
771 + orientational contribution separately (e.g. the dipole-dipole dot
772 + product is shifted by a different function than the dipole-distance
773 + products), while the hard cutoff contains no orientation-dependent
774 + shifting.  The forces and torques for these methods therefore diverge
775 + for multipoles even though the forces for point charges are identical.
776 +
777   \begin{figure}
778    \centering
779    \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps}
# Line 874 | Line 783 | forces is desired.
783      line) indicate force magnitude values indistinguishable from those
784      obtained using the multipolar Ewald sum.  Different values of the
785      cutoff radius are indicated with different symbols (9\AA\ =
786 <    circles, 12\AA\ = squares, and 15\AA\ = inverted triangles). }
787 <  \label{fig:slopeCorr_force}
786 >    circles, 12\AA\ = squares, and 15\AA\ = inverted
787 >    triangles).\label{fig:slopeCorr_force}}
788   \end{figure}
789  
790  
# Line 888 | Line 797 | forces is desired.
797      line) indicate force magnitude values indistinguishable from those
798      obtained using the multipolar Ewald sum.  Different values of the
799      cutoff radius are indicated with different symbols (9\AA\ =
800 <    circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
801 <  \label{fig:slopeCorr_torque}
800 >    circles, 12\AA\ = squares, and 15\AA\ = inverted
801 >    triangles).\label{fig:slopeCorr_torque}}
802   \end{figure}
803  
804   The torques (Fig. \ref{fig:slopeCorr_torque}) appear to be
# Line 946 | Line 855 | systematically improved by varying $\alpha$ and $r_c$.
855      indicates direction of the force or torque vectors are
856      indistinguishable from those obtained from the Ewald sum. Here
857      different symbols represent different values of the cutoff radius
858 <    (9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)}
950 <  \label{fig:slopeCorr_circularVariance}
858 >    (9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}}
859   \end{figure}
860  
861   \subsection{Energy conservation\label{sec:conservation}}
# Line 982 | Line 890 | than the multipolar Ewald sum, even when utilizing a r
890  
891   We note that for all tested values of the cutoff radius, the new
892   real-space methods can provide better energy conservation behavior
893 < than the multipolar Ewald sum, even when utilizing a relatively large
894 < $k$-space cutoff values.
893 > than the multipolar Ewald sum, even when relatively large $k$-space
894 > cutoff values are utilized.
895  
896   \begin{figure}
897    \centering
898    \includegraphics[width=\textwidth]{newDrift_12.eps}
899 < \label{fig:energyDrift}        
900 < \caption{Analysis of the energy conservation of the real-space
901 <  methods. $\delta \mathrm{E}_1$ is the linear drift in energy over
902 <  time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is
903 <  the standard deviation of energy fluctuations around this drift (in
904 <  kcal / mol / particle).  Points that appear below the dashed grey
905 <  (Ewald) lines exhibit better energy conservation than commonly-used
906 <  parameters for Ewald-based electrostatics.  All simulations were of
999 <  a 2000-molecule simulation of SSDQ water with 48 ionic charges at
1000 <  300 K starting from the same initial configuration. All runs
1001 <  utilized the same real-space cutoff, $r_c = 12$\AA.}
899 > \caption{Analysis of the energy conservation of the real-space methods
900 >  for the SSDQ water/ion system. $\delta \mathrm{E}_1$ is the linear
901 >  drift in energy over time (in kcal/mol/particle/ns) and $\delta
902 >  \mathrm{E}_0$ is the standard deviation of energy fluctuations
903 >  around this drift (in kcal/mol/particle).  Points that appear in the
904 >  green region at the bottom exhibit better energy conservation than
905 >  would be obtained using common parameters for Ewald-based
906 >  electrostatics.\label{fig:energyDrift}}
907   \end{figure}
908  
909   \subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}}
# Line 1020 | Line 925 | Fig. \ref{fig:gofr}.
925   \begin{figure}
926    \centering
927    \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
1023 – \label{fig:gofr}        
928   \caption{The pair distribution functions, $g(r)$, for the SSDQ
929    water/ion system obtained using the different real-space methods are
930    essentially identical with the result from the Ewald
931 <  treatment.}
931 >  treatment.\label{fig:gofr}}
932   \end{figure}
933  
934   There is a very slight overstructuring of the first solvation shell
# Line 1054 | Line 958 | $r_c = 12$\AA appears to be an excellent compromise fo
958   These ``optimal'' values of the damping coefficient are slightly
959   larger than what were observed for DSF electrostatics for purely
960   point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for
961 < $r_c = 12$\AA appears to be an excellent compromise for mixed charge
962 < multipole systems.
961 > $r_c = 12$\AA\ appears to be an excellent compromise for mixed
962 > charge/multipolar systems.
963  
964   To test the fidelity of the electrostatic methods at reproducing
965   dynamics in a multipolar liquid, it is also useful to look at
# Line 1082 | Line 986 | $\hat{\mathbf{A}}_\gamma$ is the space-frame unit vect
986   \label{eq:OrientCorr}
987   \end{equation}
988   where $P_l$ is the Legendre polynomial of order $l$ and
989 < $\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis
990 < $\gamma$ on a molecule..  Th body-fixed reference frame used for our
991 < models has the $z$-axis running along the dipoles, and for the SSDQ
992 < water model, the $y$-axis connects the two implied hydrogen atom
993 < positions.  From the orientation autocorrelation functions, we can
994 < obtain time constants for rotational relaxation either by fitting an
995 < exponential function or by integrating the entire correlation
996 < function.  In a good water model, these decay times would be
997 < comparable to water orientational relaxation times from nuclear
998 < magnetic resonance (NMR). The relaxation constant obtained from
999 < $C_2^y(t)$ is normally of experimental interest because it describes
1000 < the relaxation of the principle axis connecting the hydrogen
1001 < atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
1002 < of the dipole-dipole relaxation from a proton NMR signal and should
1003 < provide an estimate of the NMR relaxation time constant.\cite{Impey82}
989 > $\hat{\mathbf{A}}_\gamma$ is the unit vector for body axis $\gamma$.
990 > The reference frame used for our sample dipolar systems has the
991 > $z$-axis running along the dipoles, and for the SSDQ water model, the
992 > $y$-axis connects the two implied hydrogen atom positions.  From the
993 > orientation autocorrelation functions, we can obtain time constants
994 > for rotational relaxation either by fitting an exponential function or
995 > by integrating the entire correlation function.  In a good water
996 > model, these decay times would be comparable to water orientational
997 > relaxation times from nuclear magnetic resonance (NMR). The relaxation
998 > constant obtained from $C_2^y(t)$ is normally of experimental interest
999 > because it describes the relaxation of the principle axis connecting
1000 > the hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the
1001 > intermolecular portion of the dipole-dipole relaxation from a proton
1002 > NMR signal and should provide an estimate of the NMR relaxation time
1003 > constant.\cite{Impey82}
1004  
1005   Results for the diffusion constants and orientational relaxation times
1006   are shown in figure \ref{tab:Props}. From this data, it is apparent
# Line 1107 | Line 1011 | between the real-space and the Ewald results.
1011   between the real-space and the Ewald results.
1012  
1013   \begin{table}
1110 – \label{tab:Props}
1014   \caption{Comparison of the structural and dynamic properties for the
1015 <  soft dipolar liquid test for all of the real-space methods.}
1015 >  soft dipolar liquid test for all of the real-space methods.\label{tab:Props}}
1016   \begin{tabular}{l|c|cccc|cccc|cccc}
1017           & Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\
1018   $\alpha$ (\AA$^{-1}$) & &      
# Line 1162 | Line 1065 | We also developed two natural extensions of the damped
1065   distance that prevents its use in molecular dynamics.
1066  
1067   We also developed two natural extensions of the damped shifted-force
1068 < (DSF) model originally proposed by Fennel and
1069 < Gezelter.\cite{Fennell:2006lq} The GSF and TSF approaches provide
1070 < smooth truncation of energies, forces, and torques at the real-space
1071 < cutoff, and both converge to DSF electrostatics for point-charge
1072 < interactions.  The TSF model is based on a high-order truncated Taylor
1073 < expansion which can be relatively perturbative inside the cutoff
1074 < sphere.  The GSF model takes the gradient from an images of the
1075 < interacting multipole that has been projected onto the cutoff sphere
1076 < to derive shifted force and torque expressions, and is a significantly
1077 < more gentle approach.
1068 > (DSF) model originally proposed by Zahn {\it et al.} and extended by
1069 > Fennel and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF
1070 > approaches provide smooth truncation of energies, forces, and torques
1071 > at the real-space cutoff, and both converge to DSF electrostatics for
1072 > point-charge interactions.  The TSF model is based on a high-order
1073 > truncated Taylor expansion which can be relatively perturbative inside
1074 > the cutoff sphere.  The GSF model takes the gradient from an images of
1075 > the interacting multipole that has been projected onto the cutoff
1076 > sphere to derive shifted force and torque expressions, and is a
1077 > significantly more gentle approach.
1078  
1079 < Of the two newly-developed shifted force models, the GSF method
1080 < produced quantitative agreement with Ewald energy, force, and torques.
1081 < It also performs well in conserving energy in MD simulations.  The
1082 < Taylor-shifted (TSF) model provides smooth dynamics, but these take
1083 < place on a potential energy surface that is significantly perturbed
1084 < from Ewald-based electrostatics.  
1079 > The GSF method produced quantitative agreement with Ewald energy,
1080 > force, and torques.  It also performs well in conserving energy in MD
1081 > simulations.  The Taylor-shifted (TSF) model provides smooth dynamics,
1082 > but these take place on a potential energy surface that is
1083 > significantly perturbed from Ewald-based electrostatics.  Because it
1084 > performs relatively poorly compared with GSF, it may seem odd that
1085 > that the TSF model was included in this work.  However, the functional
1086 > forms derived for the SP and GSF methods depend on the separation of
1087 > orientational contributions that were made visible by the Taylor
1088 > series of the electrostatic kernel at the cutoff radius. The TSF
1089 > method also has the unique property that a large number of derivatives
1090 > can be made to vanish at the cutoff radius.  This property has proven
1091 > useful in past treatments of the corrections to the fluctuation
1092 > formula for dielectric constants.\cite{Izvekov:2008wo}
1093  
1094 < % The direct truncation of any electrostatic potential energy without
1095 < % multipole neutralization creates large fluctuations in molecular
1096 < % simulations.  This fluctuation in the energy is very large for the case
1097 < % of crystal because of long range of multipole ordering (Refer paper
1098 < % I).\cite{PaperI} This is also significant in the case of the liquid
1099 < % because of the local multipole ordering in the molecules. If the net
1100 < % multipole within cutoff radius neutralized within cutoff sphere by
1101 < % placing image multiples on the surface of the sphere, this fluctuation
1102 < % in the energy reduced significantly. Also, the multipole
1103 < % neutralization in the generalized SP method showed very good agreement
1104 < % with the Ewald as compared to direct truncation for the evaluation of
1194 < % the $\triangle E$ between the configurations.  In MD simulations, the
1195 < % energy conservation is very important. The conservation of the total
1196 < % energy can be ensured by i) enforcing the smooth truncation of the
1197 < % energy, force and torque in the cutoff radius and ii) making the
1198 < % energy, force and torque consistent with each other. The GSF and TSF
1199 < % methods ensure the consistency and smooth truncation of the energy,
1200 < % force and torque at the cutoff radius, as a result show very good
1201 < % total energy conservation. But the TSF method does not show good
1202 < % agreement in the absolute value of the electrostatic energy, force and
1203 < % torque with the Ewald.  The GSF method has mimicked Ewald’s force,
1204 < % energy and torque accurately and also conserved energy.
1094 > Reproduction of both structural and dynamical features in the liquid
1095 > systems is remarkably good for both the SP and GSF models.  Pair
1096 > distribution functions are essentially equivalent to the same
1097 > functions produced using Ewald-based electrostatics, and with moderate
1098 > damping, a structural feature that directly probes the electrostatic
1099 > interaction (e.g. the mean electrostatic potential energy) can also be
1100 > made quantitative.  Dynamical features are sensitive probes of the
1101 > forces and torques produced by these methods, and even though the
1102 > smooth behavior of forces is produced by perturbing the overall
1103 > potential, the diffusion constants and orientational correlation times
1104 > are quite close to the Ewald-based results.
1105  
1106   The only cases we have found where the new GSF and SP real-space
1107   methods can be problematic are those which retain a bulk dipole moment
# Line 1212 | Line 1112 | Based on the results of this work, the GSF method is a
1112   replaced by the bare electrostatic kernel, and the energies return to
1113   the expected converged values.
1114  
1115 < Based on the results of this work, the GSF method is a suitable and
1116 < efficient replacement for the Ewald sum for evaluating electrostatic
1117 < interactions in MD simulations.  Both methods retain excellent
1118 < fidelity to the Ewald energies, forces and torques.  Additionally, the
1119 < energy drift and fluctuations from the GSF electrostatics are better
1120 < than a multipolar Ewald sum for finite-sized reciprocal spaces.
1121 < Because they use real-space cutoffs with moderate cutoff radii, the
1122 < GSF and SP models approach $\mathcal{O}(N)$ scaling as the system size
1123 < increases.  Additionally, they can be made extremely efficient using
1224 < spline interpolations of the radial functions.  They require no
1225 < Fourier transforms or $k$-space sums, and guarantee the smooth
1226 < handling of energies, forces, and torques as multipoles cross the
1227 < real-space cutoff boundary.
1115 > Based on the results of this work, we can conclude that the GSF method
1116 > is a suitable and efficient replacement for the Ewald sum for
1117 > evaluating electrostatic interactions in modern MD simulations, and
1118 > the SP meethod would be an excellent choice for Monte Carlo
1119 > simulations where smooth forces and energy conservation are not
1120 > important.  Both the SP and GSF methods retain excellent fidelity to
1121 > the Ewald energies, forces and torques.  Additionally, the energy
1122 > drift and fluctuations from the GSF electrostatics are significantly
1123 > better than a multipolar Ewald sum for finite-sized reciprocal spaces.
1124  
1125 + As in all purely pairwise cutoff methods, the SP, GSF and TSF methods
1126 + are expected to scale approximately {\it linearly} with system size,
1127 + and are easily parallelizable.  This should result in substantial
1128 + reductions in the computational cost of performing large simulations.
1129 + With the proper use of pre-computation and spline interpolation of the
1130 + radial functions, the real-space methods are essentially the same cost
1131 + as a simple real-space cutoff.  They require no Fourier transforms or
1132 + $k$-space sums, and guarantee the smooth handling of energies, forces,
1133 + and torques as multipoles cross the real-space cutoff boundary.
1134 +
1135 + We are not suggesting that there is any flaw with the Ewald sum; in
1136 + fact, it is the standard by which the SP, GSF, and TSF methods have
1137 + been judged in this work.  However, these results provide evidence
1138 + that in the typical simulations performed today, the Ewald summation
1139 + may no longer be required to obtain the level of accuracy most
1140 + researchers have come to expect.
1141 +
1142   \begin{acknowledgments}
1143    JDG acknowledges helpful discussions with Christopher
1144    Fennell. Support for this project was provided by the National

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines