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

Comparing trunk/chrisDissertation/Electrostatics.tex (file contents):
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC vs.
Revision 3045 by chrisfen, Fri Oct 13 21:03:14 2006 UTC

# Line 1 | Line 1
1 < \chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES}
1 > \chapter[ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES]{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION TECHNIQUES}
2  
3   In molecular simulations, proper accumulation of electrostatic
4 < interactions is essential and is one of the most
4 > interactions is essential and one of the most
5   computationally-demanding tasks.  The common molecular mechanics force
6   fields represent atomic sites with full or partial charges protected
7   by repulsive Lennard-Jones interactions.  This means that nearly every
# Line 32 | Line 32 | Wolf {\it et al.},\cite{Wolf99} which we further exten
32  
33   In this chapter, we focus on a new set of pairwise methods devised by
34   Wolf {\it et al.},\cite{Wolf99} which we further extend.  These
35 < methods along with a few other mixed methods (i.e. reaction field) are
36 < compared with the smooth particle mesh Ewald
35 > methods, along with a few other mixed methods (i.e. reaction field),
36 > are compared with the smooth particle mesh Ewald
37   sum,\cite{Onsager36,Essmann99} which is our reference method for
38   handling long-range electrostatic interactions. The new methods for
39   handling electrostatics have the potential to scale linearly with
40 < increasing system size since they involve only a simple modification
40 > increasing system size, since they involve only a simple modification
41   to the direct pairwise sum.  They also lack the added periodicity of
42   the Ewald sum, so they can be used for systems which are non-periodic
43   or which have one- or two-dimensional periodicity.  Below, these
44 < methods are evaluated using a variety of model systems to
45 < establish their usability in molecular simulations.
44 > methods are evaluated using a variety of model systems to establish
45 > their usability in molecular simulations.
46  
47   \section{The Ewald Sum}
48  
# Line 92 | Line 92 | $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the
92   where $\alpha$ is the damping or convergence parameter with units of
93   \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
94   $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
95 < constant of the surrounding medium. The final two terms of
96 < equation (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
97 < for interacting with a surrounding dielectric.\cite{Allen87} This
98 < dipolar term was neglected in early applications in molecular
99 < simulations,\cite{Brush66,Woodcock71} until it was introduced by de
100 < Leeuw {\it et al.} to address situations where the unit cell has a
101 < dipole moment which is magnified through replication of the periodic
102 < images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
103 < system is said to be using conducting (or ``tin-foil'') boundary
104 < conditions, $\epsilon_{\rm S} = \infty$.
95 > constant of the surrounding medium. The final two terms of equation
96 > (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term for
97 > interacting with a surrounding dielectric.\cite{Allen87} This dipolar
98 > term was neglected in early applications of this technique in
99 > molecular simulations,\cite{Brush66,Woodcock71} until it was
100 > introduced by de Leeuw {\it et al.} to address situations where the
101 > unit cell has a dipole moment which is magnified through replication
102 > of the periodic images.\cite{deLeeuw80,Smith81} If this term is taken
103 > to be zero, the system is said to be using conducting (or
104 > ``tin-foil'') boundary conditions, $\epsilon_{\rm S} = \infty$.
105  
106   \begin{figure}
107   \includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
# Line 155 | Line 155 | Several studies have recognized that the inherent peri
155   bringing them more in line with the cost of the full 3-D summation.
156  
157   Several studies have recognized that the inherent periodicity in the
158 < Ewald sum can also have an effect on three-dimensional
159 < systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
160 < Solvated proteins are essentially kept at high concentration due to
161 < the periodicity of the electrostatic summation method.  In these
162 < systems, the more compact folded states of a protein can be
163 < artificially stabilized by the periodic replicas introduced by the
164 < Ewald summation.\cite{Weber00} Thus, care must be taken when
158 > Ewald sum can have an effect not just on reduced dimensionality
159 > system, but on three-dimensional systems as
160 > well.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
161 > As an example, solvated proteins are essentially kept at high
162 > concentration due to the periodicity of the electrostatic summation
163 > method.  In these systems, the more compact folded states of a protein
164 > can be artificially stabilized by the periodic replicas introduced by
165 > the Ewald summation.\cite{Weber00} Thus, care must be taken when
166   considering the use of the Ewald summation where the assumed
167   periodicity would introduce spurious effects.
168  
# Line 174 | Line 175 | short-ranged in condensed phase systems and that neutr
175   periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.}
176   observed that the electrostatic interaction is effectively
177   short-ranged in condensed phase systems and that neutralization of the
178 < charge contained within the cutoff radius is crucial for potential
178 > charges contained within the cutoff radius is crucial for potential
179   stability. They devised a pairwise summation method that ensures
180   charge neutrality and gives results similar to those obtained with the
181   Ewald summation.  The resulting shifted Coulomb potential includes
182   image-charges subtracted out through placement on the cutoff sphere
183   and a distance-dependent damping function (identical to that seen in
184 < the real-space portion of the Ewald sum) to aid convergence
184 > the real-space portion of the Ewald sum) to aid convergence:
185   \begin{equation}
186   V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}
187   - \lim_{r_{ij}\rightarrow R_\textrm{c}}
# Line 240 | Line 241 | and showed that this potential does fairly well at cap
241   \label{eq:ZahnPot}
242   \end{equation}
243   and showed that this potential does fairly well at capturing the
244 < structural and dynamic properties of water compared the same
244 > structural and dynamic properties of water compared with the same
245   properties obtained using the Ewald sum.
246  
247   \section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation}
# Line 311 | Line 312 | v(r) = \frac{q_i q_j}{r},
312   v(r) = \frac{q_i q_j}{r},
313   \label{eq:Coulomb}
314   \end{equation}
315 < then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
316 < al.}'s undamped prescription:
315 > then the {\sc sp} form will give Wolf {\it et al.}'s undamped
316 > prescription:
317   \begin{equation}
318   V_\textrm{SP}(r) =
319   q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
# Line 334 | Line 335 | simulations.
335   forces at the cutoff radius which results in energy drift during MD
336   simulations.
337  
338 < The shifted force ({\sc sf}) form using the normal Coulomb potential
338 < will give,
338 > The {\sc sf} form using the normal Coulomb potential will give,
339   \begin{equation}
340   V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}
341   + \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right]
# Line 349 | Line 349 | This formulation has the benefits that there are no di
349   \label{eq:SFForces}
350   \end{equation}
351   This formulation has the benefits that there are no discontinuities at
352 < the cutoff radius, while the neutralizing image charges are present in
352 > the cutoff radius and the neutralizing image charges are present in
353   both the energy and force expressions.  It would be simple to add the
354   self-neutralizing term back when computing the total energy of the
355   system, thereby maintaining the agreement with the Madelung energies.
# Line 361 | Line 361 | insufficient for accurate determination of the energy
361   Wolf \textit{et al.} originally discussed the energetics of the
362   shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was
363   insufficient for accurate determination of the energy with reasonable
364 < cutoff distances.  The calculated Madelung energies fluctuated around
365 < the expected value as the cutoff radius was increased, but the
364 > cutoff distances.  The calculated Madelung energies fluctuated wildly
365 > around the expected value, but as the cutoff radius was increased, the
366   oscillations converged toward the correct value.\cite{Wolf99} A
367 < damping function was incorporated to accelerate the convergence; and
367 > damping function was incorporated to accelerate this convergence; and
368   though alternative forms for the damping function could be
369   used,\cite{Jones56,Heyes81} the complimentary error function was
370   chosen to mirror the effective screening used in the Ewald summation.
# Line 374 | Line 374 | v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
374   v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
375   \label{eq:dampCoulomb}
376   \end{equation}
377 < the shifted potential (Eq. (\ref{eq:SPPot})) becomes
377 > the {\sc sp} potential function (Eq. (\ref{eq:SPPot})) becomes
378   \begin{equation}
379   V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}
380   - \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right)
# Line 389 | Line 389 | + \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r
389   \quad r\leqslant R_\textrm{c}.
390   \label{eq:DSPForces}
391   \end{equation}
392 < Again, this damped shifted potential suffers from a
393 < force-discontinuity at the cutoff radius, and the image charges play
394 < no role in the forces.  To remedy these concerns, one may derive a
395 < {\sc sf} variant by including the derivative term in
396 < equation (\ref{eq:shiftingForm}),
392 > Again, this damped shifted potential suffers from a discontinuity in
393 > the forces at the cutoff radius, and the image charges play no role in
394 > the forces.  To remedy these concerns, one may derive a {\sc sf}
395 > variant by including the derivative term present in
396 > equation~(\ref{eq:shiftingForm}),
397   \begin{equation}
398   \begin{split}
399   V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
# Line 423 | Line 423 | If the damping parameter $(\alpha)$ is set to zero, th
423   \end{split}
424   \end{equation}
425   If the damping parameter $(\alpha)$ is set to zero, the undamped case,
426 < equations (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
427 < recovered from equations (\ref{eq:DSPPot} through \ref{eq:DSFForces}).
426 > equations (\ref{eq:SPPot}) through (\ref{eq:SFForces}) are correctly
427 > recovered from equations (\ref{eq:DSPPot}) through (\ref{eq:DSFForces}).
428  
429 < This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
429 > This new {\sc sf} potential is similar to equation (\ref{eq:ZahnPot})
430   derived by Zahn \textit{et al.}; however, there are two important
431   differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation
432   (\ref{eq:shiftingForm}) is equal to equation (\ref{eq:dampCoulomb})
# Line 436 | Line 436 | would be expected to have sudden jumps as particle dis
436   portion is different.  The missing $v_\textrm{c}$ term would not
437   affect molecular dynamics simulations (although the computed energy
438   would be expected to have sudden jumps as particle distances crossed
439 < $R_c$).  The sign problem is a potential source of errors, however.
440 < In fact, it introduces a discontinuity in the forces at the cutoff,
441 < because the force function is shifted in the wrong direction and
442 < doesn't cross zero at $R_\textrm{c}$.
439 > $R_c$); however, the sign problem is a potential source of errors.  In
440 > fact, equation~(\ref{eq:ZahnPot}) introduces a discontinuity in the
441 > forces at the cutoff, because the force function is shifted in the
442 > wrong direction and does not cross zero at $R_\textrm{c}$.
443  
444   Equations (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
445   electrostatic summation method in which the potential and forces are
446   continuous at the cutoff radius and which incorporates the damping
447   function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
448 < this paper, we will evaluate exactly how good these methods ({\sc sp},
449 < {\sc sf}, damping) are at reproducing the correct electrostatic
448 > this chapter, we will evaluate exactly how good these methods ({\sc
449 > sp}, {\sc sf}, damping) are at reproducing the correct electrostatic
450   summation performed by the Ewald sum.
451  
452  
# Line 455 | Line 455 | classical molecular mechanics simulations: Monte Carlo
455   As mentioned in the introduction, there are two primary techniques
456   utilized to obtain information about the system of interest in
457   classical molecular mechanics simulations: Monte Carlo (MC) and
458 < Molecular Dynamics (MD).  Both of these techniques utilize pairwise
458 > molecular dynamics (MD).  Both of these techniques utilize pairwise
459   summations of interactions between particle sites, but they use these
460   summations in different ways.
461  
# Line 476 | Line 476 | cumulative, one should expect greater deviation at lon
476   electrostatic summation techniques, the dynamics in the short term
477   will be indistinguishable.  Because error in MD calculations is
478   cumulative, one should expect greater deviation at longer times,
479 < although methods which have large differences in the force and torque
479 > and methods which have large differences in the force and torque
480   vectors will diverge from each other more rapidly.
481  
482   \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
# Line 503 | Line 503 | Sample correlation plots for two alternate methods are
503   correlation (slope) and correlation coefficient for these regressions
504   indicate perfect agreement between the alternative method and {\sc spme}.
505   Sample correlation plots for two alternate methods are shown in
506 < Fig. \ref{fig:linearFit}.
506 > figure \ref{fig:linearFit}.
507  
508   Each of the seven system types (detailed in section \ref{sec:RepSims})
509   were represented using 500 independent configurations.  Thus, each of
510   the alternative (non-Ewald) electrostatic summation methods was
511   evaluated using an accumulated 873,250 configurational energy
512 < differences.
513 <
514 < Results and discussion for the individual analysis of each of the
515 < system types appear in appendix \ref{app:IndividualResults}, while the
516 < cumulative results over all the investigated systems appear below in
517 < sections \ref{sec:EnergyResults}.
512 > differences. Results for and discussions regarding the individual
513 > analysis of each of the system types appear in appendix
514 > \ref{app:IndividualResults}, while the cumulative results over all the
515 > investigated systems appear below in section~\ref{sec:EnergyResults}.
516  
517   \subsection{Molecular Dynamics and the Force and Torque
518   Vectors}\label{sec:MDMethods} We evaluated the pairwise methods
# Line 527 | Line 525 | forces (and torques) on each molecule in each configur
525   comparing $\Delta E$ values.  Instead of a single energy difference
526   between two system configurations, we compared the magnitudes of the
527   forces (and torques) on each molecule in each configuration.  For a
528 < system of 1000 water molecules and 40 ions, there are 1040 force
529 < vectors and 1000 torque vectors.  With 500 configurations, this
530 < results in 520,000 force and 500,000 torque vector comparisons.
531 < Additionally, data from seven different system types was aggregated
532 < before the comparison was made.
528 > system of 1000 water molecules and 40 ions, there are 1040 force and
529 > 1000 torque vectors.  With 500 configurations, this results in 520,000
530 > force and 500,000 torque vector comparisons.  Additionally, data from
531 > seven different system types was aggregated before comparisons were
532 > made.
533  
534   The {\it directionality} of the force and torque vectors was
535   investigated through measurement of the angle ($\theta$) formed
# Line 546 | Line 544 | between two different electrostatic summation methods,
544   unit sphere.  Since this distribution is a measure of angular error
545   between two different electrostatic summation methods, there is no
546   {\it a priori} reason for the profile to adhere to any specific
547 < shape. Thus, gaussian fits were used to measure the width of the
547 > shape. Thus, Gaussian fits were used to measure the width of the
548   resulting distributions. The variance ($\sigma^2$) was extracted from
549   each of these fits and was used to compare distribution widths.
550   Values of $\sigma^2$ near zero indicate vector directions
# Line 609 | Line 607 | crystals), so the systems studied were:
607   \item a high ionic strength solution of NaCl in water (1.1 M), and
608   \item a 6~\AA\  radius sphere of Argon in water.
609   \end{enumerate}
612 –
610   By utilizing the pairwise techniques (outlined in section
611   \ref{sec:ESMethods}) in systems composed entirely of neutral groups,
612   charged particles, and mixtures of the two, we hope to discern under
# Line 733 | Line 730 | limitations, primarily that it was developed for use i
730  
731   The reaction field results illustrates some of that method's
732   limitations, primarily that it was developed for use in homogeneous
733 < systems; although it does provide results that are an improvement over
734 < those from an unmodified cutoff.
733 > systems. It does, however, provide results that are an improvement
734 > over those from an unmodified cutoff.
735  
736   \section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults}
737  
# Line 757 | Line 754 | inverted triangles).}
754   inverted triangles).}
755   \label{fig:frcMag}
756   \end{figure}
757 < Again, it is striking how well the Shifted Potential and Shifted Force
758 < methods are doing at reproducing the {\sc spme} forces.  The undamped and
759 < weakly-damped {\sc sf} method gives the best agreement with Ewald.
760 < This is perhaps expected because this method explicitly incorporates a
761 < smooth transition in the forces at the cutoff radius as well as the
757 > Again, it is striking how well the {\sc sp} and {\sc sf} methods
758 > reproduce the {\sc spme} forces.  The undamped and weakly-damped {\sc
759 > sf} method gives the best agreement with Ewald.  This is perhaps
760 > expected because this method explicitly incorporates a smooth
761 > transition in the forces at the cutoff radius as well as the
762   neutralizing image charges.
763  
764   Figure \ref{fig:frcMag}, for the most part, parallels the results seen
# Line 772 | Line 769 | method is generating usable forces.  Further increases
769  
770   With moderate damping and a large enough cutoff radius, the {\sc sp}
771   method is generating usable forces.  Further increases in damping,
772 < while beneficial for simulations with a cutoff radius of 9~\AA\ , is
772 > while beneficial for simulations with a cutoff radius of 9~\AA\ , are
773   detrimental to simulations with larger cutoff radii.
774  
775   The reaction field results are surprisingly good, considering the poor
776   quality of the fits for the $\Delta E$ results.  There is still a
777 < considerable degree of scatter in the data, but the forces correlate
778 < well with the Ewald forces in general.  We note that the reaction
779 < field calculations do not include the pure NaCl systems, so these
780 < results are partly biased towards conditions in which the method
781 < performs more favorably.
777 > considerable degree of scatter in the data, but in general, the forces
778 > correlate well with the Ewald forces.  We note that the pure NaCl
779 > systems were not included in the system set used in the reaction field
780 > calculations, so these results are partly biased towards conditions in
781 > which the method performs more favorably.
782  
783   \begin{figure}
784   \centering
# Line 797 | Line 794 | rigid molecules (i.e. the systems containing water).  
794   \end{figure}
795   Molecular torques were only available from the systems which contained
796   rigid molecules (i.e. the systems containing water).  The data in
797 < fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
797 > figure \ref{fig:trqMag} is taken from this smaller sampling pool.
798  
799 < Torques appear to be much more sensitive to charges at a longer
800 < distance.   The striking feature in comparing the new electrostatic
801 < methods with {\sc spme} is how much the agreement improves with increasing
802 < cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
803 < appears to reproduce the {\sc spme} torques most accurately.  
799 > Torques appear to be much more sensitive to charge interactions at
800 > longer distances.  The most noticeable feature in comparing the new
801 > electrostatic methods with {\sc spme} is how much the agreement
802 > improves with increasing cutoff radius.  Again, the weakly damped and
803 > undamped {\sc sf} method appears to reproduce the {\sc spme} torques
804 > most accurately.
805  
806   Water molecules are dipolar, and the reaction field method reproduces
807   the effect of the surrounding polarized medium on each of the
# Line 818 | Line 816 | angles formed between these vectors and the same vecto
816   will also be vital in calculating dynamical quantities accurately.
817   Force and torque directionalities were investigated by measuring the
818   angles formed between these vectors and the same vectors calculated
819 < using {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared
819 > using {\sc spme}.  The results (figure \ref{fig:frcTrqAng}) are compared
820   through the variance ($\sigma^2$) of the Gaussian fits of the angle
821   error distributions of the combined set over all system types.
822  
# Line 837 | Line 835 | total accumulated system data are tabulated in figure
835   \end{figure}
836   Both the force and torque $\sigma^2$ results from the analysis of the
837   total accumulated system data are tabulated in figure
838 < \ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
839 < sp}) method would be essentially unusable for molecular dynamics
840 < unless the damping function is added.  The Shifted Force ({\sc sf})
841 < method, however, is generating force and torque vectors which are
842 < within a few degrees of the Ewald results even with weak (or no)
845 < damping.
838 > \ref{fig:frcTrqAng}. Here it is clear that the {\sc sp} method would
839 > be essentially unusable for molecular dynamics unless the damping
840 > function is added.  The {\sc sf} method, however, is generating force
841 > and torque vectors which are within a few degrees of the Ewald results
842 > even with weak (or no) damping.
843  
844   All of the sets (aside from the over-damped case) show the improvement
845   afforded by choosing a larger cutoff radius.  Increasing the cutoff
# Line 924 | Line 921 | However, at larger values of $\alpha$, it is possible
921   The complimentary error function inserted into the potential weakens
922   the electrostatic interaction as the value of $\alpha$ is increased.
923   However, at larger values of $\alpha$, it is possible to over-damp the
924 < electrostatic interaction and to remove it completely.  Kast
924 > electrostatic interaction and remove it completely.  Kast
925   \textit{et al.}  developed a method for choosing appropriate $\alpha$
926   values for these types of electrostatic summation methods by fitting
927   to $g(r)$ data, and their methods indicate optimal values of 0.34,
928   0.25, and 0.16~\AA$^{-1}$ for cutoff values of 9, 12, and 15~\AA\
929   respectively.\cite{Kast03} These appear to be reasonable choices to
930 < obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on
930 > obtain proper MC behavior (figure \ref{fig:delE}); however, based on
931   these findings, choices this high would introduce error in the
932 < molecular torques, particularly for the shorter cutoffs.  Based on our
933 < observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial,
934 < but damping may be unnecessary when using the {\sc sf} method.
932 > molecular torques, particularly for the shorter cutoffs.  Based on the
933 > above observations, empirical damping up to 0.2~\AA$^{-1}$ is
934 > beneficial, but damping may be unnecessary when using the {\sc sf}
935 > method.
936  
937  
938   \section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics}
# Line 943 | Line 941 | that a method similar (but not identical with) the dam
941   using equations (\ref{eq:ZahnPot}) and
942   (\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated
943   that a method similar (but not identical with) the damped {\sc sf}
944 < method resulted in properties very similar to those obtained when
944 > method resulted in properties very close to those obtained when
945   using the Ewald summation.  The properties they studied (pair
946   distribution functions, diffusion constants, and velocity and
947   orientational correlation functions) may not be particularly sensitive
948   to the long-range and collective behavior that governs the
949   low-frequency behavior in crystalline systems.  Additionally, the
950 < ionic crystals are the worst case scenario for the pairwise methods
950 > ionic crystals are a worst case scenario for the pairwise methods
951   because they lack the reciprocal space contribution contained in the
952   Ewald summation.
953  
954 < We are using two separate measures to probe the effects of these
954 > We used two separate measures to probe the effects of these
955   alternative electrostatic methods on the dynamics in crystalline
956 < materials.  For short- and intermediate-time dynamics, we are
957 < computing the velocity autocorrelation function, and for long-time
958 < and large length-scale collective motions, we are looking at the
959 < low-frequency portion of the power spectrum.
956 > materials.  For short- and intermediate-time dynamics, we computed the
957 > velocity autocorrelation function, and for long-time and large
958 > length-scale collective motions, we looked at the low-frequency
959 > portion of the power spectrum.
960  
961   \begin{figure}
962   \centering
# Line 972 | Line 970 | are stiffer than the moderately damped and {\sc spme}
970   are stiffer than the moderately damped and {\sc spme} methods.}
971   \label{fig:vCorrPlot}
972   \end{figure}
973 < The short-time decay of the velocity autocorrelation function through
973 > The short-time decay of the velocity autocorrelation functions through
974   the first collision are nearly identical in figure
975   \ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
976   how the methods differ.  The undamped {\sc sf} method has deeper
977 < troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than
977 > troughs (see inset in figure \ref{fig:vCorrPlot}) and higher peaks than
978   any of the other methods.  As the damping parameter ($\alpha$) is
979   increased, these peaks are smoothed out, and the {\sc sf} method
980   approaches the {\sc spme} results.  With $\alpha$ values of 0.2~\AA$^{-1}$,
# Line 1004 | Line 1002 | the best-performing alternative methods are shown in
1002   collective long-time motion, we computed power spectra from long-time
1003   traces of the velocity autocorrelation function. The power spectra for
1004   the best-performing alternative methods are shown in
1005 < fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1005 > figure \ref{fig:methodPS}.  Apodization of the correlation functions via
1006   a cubic switching function between 40 and 50~ps was used to reduce the
1007   ringing resulting from data truncation.  This procedure had no
1008   noticeable effect on peak location or magnitude.
# Line 1016 | Line 1014 | damped {\sc sf}.  When using moderate damping ($\alpha
1014   figure \ref{fig:methodPS}), at frequencies below 100~cm$^{-1}$, the
1015   correlated motions are blue-shifted when using undamped or weakly
1016   damped {\sc sf}.  When using moderate damping ($\alpha =
1017 < 0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly
1018 < identical correlated motion to the Ewald method (which has a
1017 > 0.2$~\AA$^{-1}$), both the {\sc sf} and {\sc sp} methods produce
1018 > correlated motions nearly identical to the Ewald method (which has a
1019   convergence parameter of 0.3119~\AA$^{-1}$).  This weakening of the
1020   electrostatic interaction with increased damping explains why the
1021   long-ranged correlated motions are at lower frequencies for the
# Line 1037 | Line 1035 | spectra for a single method ({\sc sf}) with a range of
1035   \end{figure}
1036   To isolate the role of the damping constant, we have computed the
1037   spectra for a single method ({\sc sf}) with a range of damping
1038 < constants and compared this with the {\sc spme} spectrum.
1039 < Fig. \ref{fig:dampInc} shows more clearly that increasing the
1040 < electrostatic damping red-shifts the lowest frequency phonon modes.
1041 < However, even without any electrostatic damping, the {\sc sf} method
1042 < has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
1043 < Without the {\sc sf} modifications, an undamped (pure cutoff) method
1044 < would predict the lowest frequency peak near 325~cm$^{-1}$.  {\it
1045 < Most} of the collective behavior in the crystal is accurately captured
1046 < using the {\sc sf} method.  Quantitative agreement with Ewald can be
1047 < obtained using moderate damping in addition to the shifting at the
1048 < cutoff distance.
1038 > constants and compared this with the {\sc spme} spectrum.  Figure
1039 > \ref{fig:dampInc} shows more clearly that increasing the electrostatic
1040 > damping red-shifts the lowest frequency phonon modes.  However, even
1041 > without any electrostatic damping, the {\sc sf} method has at most a
1042 > 10 cm$^{-1}$ error in the lowest frequency phonon mode.  Without the
1043 > {\sc sf} modifications, an undamped (pure cutoff) method would predict
1044 > the lowest frequency peak near 325~cm$^{-1}$, an error significantly
1045 > larger than that of the undamped {\sc sf} technique.  This indicates
1046 > that {\it most} of the collective behavior in the crystal is
1047 > accurately captured using the {\sc sf} method.  Quantitative agreement
1048 > with Ewald can be obtained using moderate damping in addition to the
1049 > shifting at the cutoff distance.
1050  
1051   \section{An Application: TIP5P-E Water}\label{sec:t5peApplied}
1052  
# Line 1068 | Line 1067 | Without this correction, the pressure term on the cent
1067   long-range electrostatic
1068   correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87}
1069   Without this correction, the pressure term on the central particle
1070 < from the surroundings is missing. Because they expand to compensate
1071 < for this added pressure term when this correction is included, systems
1072 < composed of these particles tend to under-predict the density of water
1073 < under standard conditions. When using any form of long-range
1074 < electrostatic correction, it has become common practice to develop or
1075 < utilize a reparametrized water model that corrects for this
1070 > from the surroundings is missing. When this correction is included,
1071 > systems of these particles expand to compensate for this added
1072 > pressure term and under-predict the density of water under standard
1073 > conditions. When using any form of long-range electrostatic
1074 > correction, it has become common practice to develop or utilize a
1075 > reparametrized water model that corrects for this
1076   effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows
1077 < this practice and was optimized specifically for use with the Ewald
1077 > this practice and was optimized for use with the Ewald
1078   summation.\cite{Rick04} In his publication, Rick preserved the
1079   geometry and point charge magnitudes in TIP5P and focused on altering
1080 < the Lennard-Jones parameters to correct the density at
1081 < 298K.\cite{Rick04} With the density corrected, he compared common
1082 < water properties for TIP5P-E using the Ewald sum with TIP5P using a
1084 < 9~\AA\ cutoff.
1080 > the Lennard-Jones parameters to correct the density at 298~K. With the
1081 > density corrected, he compared common water properties for TIP5P-E
1082 > using the Ewald sum with TIP5P using a 9~\AA\ cutoff.
1083  
1084   In the following sections, we compared these same water properties
1085   calculated from TIP5P-E using the Ewald sum with TIP5P-E using the
# Line 1114 | Line 1112 | forces.  Since the {\sc sp} method does not modify the
1112   brackets of equation
1113   \ref{eq:MolecularPressure}) is directly dependent on the interatomic
1114   forces.  Since the {\sc sp} method does not modify the forces (see
1115 < section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
1115 > section \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
1116   will be identical to that obtained without an electrostatic
1117   correction.  The {\sc sf} method does alter the virial component and,
1118   by way of the modified pressures, should provide densities more in
# Line 1128 | Line 1126 | method for accumulating statistics, these sequences we
1126   temperatures. The average densities were calculated from the later
1127   three-fourths of each trajectory. Similar to Mahoney and Jorgensen's
1128   method for accumulating statistics, these sequences were spliced into
1129 < 200 segments to calculate the average density and standard deviation
1130 < at each temperature.\cite{Mahoney00}
1129 > 200 segments, each providing an average density. These 200 density
1130 > values were used to calculate the average and standard deviation of
1131 > the density at each temperature.\cite{Mahoney00}
1132  
1133   \begin{figure}
1134   \includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf}
# Line 1149 | Line 1148 | technique are close to, though typically lower than, t
1148   TIP5P-E using differing electrostatic corrections overlaid on the
1149   experimental values.\cite{CRC80} The densities when using the {\sc sf}
1150   technique are close to, though typically lower than, those calculated
1151 < while using the Ewald summation. These slightly reduced densities
1152 < indicate that the pressure component from the image charges at
1153 < R$_\textrm{c}$ is larger than that exerted by the reciprocal-space
1154 < portion of the Ewald summation. Bringing the image charges closer to
1155 < the central particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than
1156 < the preferred 12~\AA\ R$_\textrm{c}$) increases the strength of their
1157 < interactions, resulting in a further reduction of the densities.
1151 > using the Ewald summation. These slightly reduced densities indicate
1152 > that the pressure component from the image charges at R$_\textrm{c}$
1153 > is larger than that exerted by the reciprocal-space portion of the
1154 > Ewald summation. Bringing the image charges closer to the central
1155 > particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the
1156 > preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image
1157 > charge interactions on the central particle and results in a further
1158 > reduction of the densities.
1159  
1160   Because the strength of the image charge interactions has a noticeable
1161   effect on the density, we would expect the use of electrostatic
# Line 1211 | Line 1211 | check whether the choice of using the Ewald summation
1211   non-polarizable models.\cite{Sorenson00} This excellent agreement with
1212   experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To
1213   check whether the choice of using the Ewald summation or the {\sc sf}
1214 < technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K
1215 < and 1atm were determined for the systems compared in the previous
1216 < section.
1214 > technique alters the liquid structure, the $g_\textrm{OO}(r)$s at
1215 > 298~K and 1~atm were determined for the systems compared in the
1216 > previous section.
1217  
1218   \begin{figure}
1219   \includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf}
# Line 1226 | Line 1226 | sf} technique with a various parameters are overlaid o
1226   \end{figure}
1227   The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc
1228   sf} technique with a various parameters are overlaid on the
1229 < $g_\textrm{OO}(r)$ while using the Ewald summation in figure
1230 < \ref{fig:t5peGofRs}. The differences in density do not appear to have
1231 < any effect on the liquid structure as the $g_\textrm{OO}(r)$s are
1232 < indistinguishable. These results indicate that the $g_\textrm{OO}(r)$
1233 < is insensitive to the choice of electrostatic correction.
1229 > $g_\textrm{OO}(r)$ while using the Ewald summation in
1230 > figure~\ref{fig:t5peGofRs}. The differences in density do not appear
1231 > to have any effect on the liquid structure as the $g_\textrm{OO}(r)$s
1232 > are indistinguishable. These results indicate that the
1233 > $g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic
1234 > correction.
1235  
1236   \subsection{Thermodynamic Properties}\label{sec:t5peThermo}
1237  
# Line 1245 | Line 1246 | The $\Delta H_\textrm{vap}$ is the enthalpy change req
1246   good set for comparisons involving the {\sc sf} technique.
1247  
1248   The $\Delta H_\textrm{vap}$ is the enthalpy change required to
1249 < transform one mol of substance from the liquid phase to the gas
1249 > transform one mole of substance from the liquid phase to the gas
1250   phase.\cite{Berry00} In molecular simulations, this quantity can be
1251   determined via
1252   \begin{equation}
# Line 1360 | Line 1361 | property trends with temperature seen when using the E
1361  
1362   As observed for the density in section \ref{sec:t5peDensity}, the
1363   property trends with temperature seen when using the Ewald summation
1364 < are reproduced with the {\sc sf} technique. One noticable difference
1364 > are reproduced with the {\sc sf} technique. One noticeable difference
1365   between the properties calculated using the two methods are the lower
1366   $\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
1367   expected due to the direct weakening of the electrostatic interaction
# Line 1371 | Line 1372 | inflated $C_p$ values at all temperatures.
1372   steeper than the experimental trend, indirectly resulting in the
1373   inflated $C_p$ values at all temperatures.
1374  
1375 < Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$
1376 < values all overlap within error. As indicated for the $\Delta
1377 < H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph,
1378 < the deviations between experiment and simulation in this region are
1379 < not the fault of the electrostatic summation methods but are due to
1380 < the TIP5P class model itself. Like most rigid, non-polarizable,
1381 < point-charge water models, the density decreases with temperature at a
1382 < much faster rate than experiment (see figure
1383 < \ref{fig:t5peDensities}). The reduced density leads to the inflated
1375 > Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values
1376 > all overlap within error. As indicated for the $\Delta H_\textrm{vap}$
1377 > and $C_p$ results discussed in the previous paragraph, the deviations
1378 > between experiment and simulation in this region are not the fault of
1379 > the electrostatic summation methods but are due to the geometry and
1380 > parameters of the TIP5P class of water models. Like most rigid,
1381 > non-polarizable, point-charge water models, the density decreases with
1382 > temperature at a much faster rate than experiment (see figure
1383 > \ref{fig:t5peDensities}). This reduced density leads to the inflated
1384   compressibility and expansivity values at higher temperatures seen
1385   here in figure \ref{fig:t5peThermo}. Incorporation of polarizability
1386 < and many-body effects are required in order for simulation to overcome
1387 < these differences with experiment.\cite{Laasonen93,Donchev06}
1386 > and many-body effects are required in order for water models to
1387 > overcome differences between simulation-based and experimentally
1388 > determined densities at these higher
1389 > temperatures.\cite{Laasonen93,Donchev06}
1390  
1391   At temperatures below the freezing point for experimental water, the
1392   differences between {\sc sf} and the Ewald summation results are more
# Line 1392 | Line 1395 | temperatures below 250~K in the {\sc sf} simulations,
1395   particularly in the case of {\sc sf} without damping. This points to
1396   the onset of a more frustrated or glassy behavior for TIP5P-E at
1397   temperatures below 250~K in the {\sc sf} simulations, indicating that
1398 < disorder in the reciprical-space term of the Ewald summation might act
1398 > disorder in the reciprocal-space term of the Ewald summation might act
1399   to loosen up the local structure more than the image-charges in {\sc
1400 < sf}. Because the systems are locked in different regions of
1401 < phase-space, comparisons between properties at these temperatures are
1402 < not exactly fair. This observation is explored in more detail in
1403 < section \ref{sec:t5peDynamics}.
1400 > sf}. The damped {\sc sf} actually makes a better comparison with
1401 > experiment in this region, particularly for the $\alpha_p$ values. The
1402 > local interactions in the undamped {\sc sf} technique appear to be too
1403 > strong since the property change is much more dramatic than the damped
1404 > forms, while the Ewald summation appears to weight the
1405 > reciprocal-space interactions at the expense the local interactions,
1406 > disagreeing with the experimental results. This observation is
1407 > explored in more detail in section \ref{sec:t5peDynamics}.
1408  
1409   The final thermodynamic property displayed in figure
1410   \ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
# Line 1407 | Line 1414 | simulations. Lack of a damping function results in die
1414   conditions.\cite{Neumann80,Neumann83} This is readily apparent in the
1415   converged $\epsilon$ values accumulated for the {\sc sf}
1416   simulations. Lack of a damping function results in dielectric
1417 < constants significantly smaller than that obtained using the Ewald
1417 > constants significantly smaller than those obtained using the Ewald
1418   sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the
1419   agreement considerably. It should be noted that the choice of the
1420   ``Ewald coefficient'' value also has a significant effect on the
# Line 1419 | Line 1426 | further explored as optimal damping coefficients for d
1426   sf}; however, the choice of cutoff radius also plays an important
1427   role. In section \ref{sec:dampingDielectric}, this connection is
1428   further explored as optimal damping coefficients for different choices
1429 < of $R_\textrm{c}$ are determined for {\sc sf} for capturing the
1430 < dielectric behavior.
1429 > of $R_\textrm{c}$ are determined for {\sc sf} in order to best capture
1430 > the dielectric behavior.
1431  
1432   \subsection{Dynamic Properties}\label{sec:t5peDynamics}
1433  
1434   To look at the dynamic properties of TIP5P-E when using the {\sc sf}
1435 < method, 200~ps $NVE$ simulations were performed for each temperature at
1436 < the average density reported by the $NPT$ simulations. The
1437 < self-diffusion constants ($D$) were calculated with the Einstein
1438 < relation using the mean square displacement (MSD),
1435 > method, 200~ps $NVE$ simulations were performed for each temperature
1436 > at the average density reported by the $NPT$ simulations. The
1437 > self-diffusion constants ($D$) were calculated using the mean square
1438 > displacement (MSD) form of the Einstein relation,
1439   \begin{equation}
1440   D = \lim_{t\rightarrow\infty}
1441      \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
# Line 1444 | Line 1451 | regions:
1451   \item linear diffusive regime, and
1452   \item a region with poor statistics.
1453   \end{enumerate}
1454 < The slope from the linear region (region 2) is used to calculate $D$.
1454 > The slope from the linear regime (region 2) is used to calculate $D$.
1455   \begin{figure}
1456   \centering
1457   \includegraphics[width=3.5in]{./figures/ExampleMSD.pdf}
# Line 1463 | Line 1470 | labeled frame axes.}
1470   labeled frame axes.}
1471   \label{fig:waterFrame}
1472   \end{figure}
1473 < In addition to translational diffusion, reorientational time constants
1473 > In addition to translational diffusion, orientational relaxation times
1474   were calculated for comparisons with the Ewald simulations and with
1475 < experiments. These values were determined from 25~ps $NVE$ trajectories
1476 < through calculation of the orientational time correlation function,
1475 > experiments. These values were determined from 25~ps $NVE$
1476 > trajectories through calculation of the orientational time correlation
1477 > function,
1478   \begin{equation}
1479   C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t)
1480                  \cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle,
# Line 1533 | Line 1541 | of the choice of damping constant. Their are several p
1541   relaxes faster than experiment with the Ewald sum while tracking
1542   experiment fairly well when using the {\sc sf} technique, independent
1543   of the choice of damping constant. Their are several possible reasons
1544 < for this deviation between techniques. The Ewald results were taken
1545 < shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick
1546 < calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of
1547 < 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$,
1548 < placing the result more in line with that obtained using the Ewald
1549 < sum. These results support this explanation; however, recomputing the
1550 < results to meet a poorer statistical standard is
1544 > for this deviation between techniques. The Ewald results were
1545 > calculated using shorter (10ps) trajectories than the {\sc sf} results
1546 > (25ps). A quick calculation from a 10~ps trajectory with {\sc sf} with
1547 > an $\alpha$ of 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in
1548 > $\tau_2^y$, placing the result more in line with that obtained using
1549 > the Ewald sum. This example supports this explanation; however,
1550 > recomputing the results to meet a poorer statistical standard is
1551   counter-productive. Assuming the Ewald results are not the product of
1552   poor statistics, differences in techniques to integrate the
1553   orientational motion could also play a role. {\sc shake} is the most
1554   commonly used technique for approximating rigid-body orientational
1555 < motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and
1555 > motion,\cite{Ryckaert77} whereas in {\sc oopse}, we maintain and
1556   integrate the entire rotation matrix using the {\sc dlm}
1557   method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
1558   technique, if the convergence tolerances are raised for increased
# Line 1572 | Line 1580 | monopole (and use the monopole potential of equation (
1580   multipoles. In a mixed system of monopoles and multipoles, the
1581   undamped {\sc sf} potential needs only to shift the force terms of the
1582   monopole (and use the monopole potential of equation (\ref{eq:SFPot}))
1583 < and smoothly cutoff the multipole interactions with a switching
1583 > and smoothly truncate the multipole interactions with a switching
1584   function. The switching function is required in order to conserve
1585 < energy, because a discontinuity will exist at $R_\textrm{c}$ in the
1586 < absence of shifting terms.
1585 > energy, because a discontinuity will exist in both the potential and
1586 > forces at $R_\textrm{c}$ in the absence of shifting terms.
1587  
1588   If we consider damping the {\sc sf} potential (Eq. (\ref{eq:DSFPot})),
1589   then we need to incorporate the complimentary error function term into
# Line 1666 | Line 1674 | term. Continuing with higher rank tensors, we can obta
1674   \end{equation}
1675   Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional
1676   term. Continuing with higher rank tensors, we can obtain the damping
1677 < functions for higher multipoles as well as the forces. Each subsequent
1677 > functions for higher multipole potentials and forces. Each subsequent
1678   damping function includes one additional term, and we can simplify the
1679   procedure for obtaining these terms by writing out the following
1680   generating function,
# Line 1700 | Line 1708 | V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathb
1708                  c_1(r_{ij}),
1709   \label{eq:dampDipoleDipole}
1710   \end{equation}
1711 < $c_2(r_{ij})$ and $c_1(r_{ij})$ respectively dampen these two
1712 < parts. The forces for the damped dipole-dipole interaction,
1711 > $c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts
1712 > respectively. The forces for the damped dipole-dipole interaction,
1713   \begin{equation}
1714   \begin{split}
1715   F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
# Line 1731 | Line 1739 | going to be quite sensitive to the choice of damping p
1739   constant is calculated from the long-time fluctuations of the system's
1740   accumulated dipole moment (Eq. (\ref{eq:staticDielectric})), so it is
1741   going to be quite sensitive to the choice of damping parameter. We
1742 < would like to choose an optimal damping constant for any particular
1743 < cutoff radius choice that would properly capture the dielectric
1744 < behavior of the liquid.
1742 > would like to choose optimal damping constants such that any arbitrary
1743 > choice of cutoff radius will properly capture the dielectric behavior
1744 > of the liquid.
1745  
1746   In order to find these optimal values, we mapped out the static
1747   dielectric constant as a function of both the damping parameter and
# Line 1744 | Line 1752 | reaction field modified variant of the soft sticky dip
1752   four-point transferable intermolecular potential (TIP4P) for water
1753   targeted for use with the Ewald summation.\cite{Horn04} SSD/RF is the
1754   reaction field modified variant of the soft sticky dipole (SSD) model
1755 < for water\cite{Fennell04} This model is discussed in more detail in
1755 > for water.\cite{Fennell04} This model is discussed in more detail in
1756   the next chapter. One thing to note about it, electrostatic
1757   interactions are handled via dipole-dipole interactions rather than
1758   charge-charge interactions like the other three models. Damping of the
# Line 1777 | Line 1785 | $\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\
1785   with {\sc sf} these parameters give a dielectric constant of
1786   90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where
1787   $\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and
1788 < 0.35~\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal
1789 < to 63$\pm$1.\cite{Horn04} We did not perform calculations with these
1790 < exact parameters, but interpolating between surrounding values gives a
1791 < $\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric
1792 < constant on $\alpha$ and $R_\textrm{c}$ with the {\sc sf} technique,
1793 < it might be interesting to investigate the dielectric dependence of
1794 < the real-space Ewald parameters.
1788 > 0.35~\AA$^{-1}$, and these parameters resulted in a dielectric
1789 > constant equal to 63$\pm$1.\cite{Horn04} We did not perform
1790 > calculations with these exact parameters, but interpolating between
1791 > surrounding values gives a dielectric constant of 61$\pm$1. Since the
1792 > dielectric constant is dependent on $\alpha$ and $R_\textrm{c}$ with
1793 > the {\sc sf} technique, it might be interesting to investigate the
1794 > dielectric dependence of the real-space Ewald parameters.
1795  
1796   Although it is tempting to choose damping parameters equivalent to
1797   these Ewald examples, the results discussed in sections
# Line 1817 | Line 1825 | employing lattice summation techniques.  The cumulativ
1825   (\ref{eq:DSFForces}), shows a remarkable ability to reproduce the
1826   energetic and dynamic characteristics exhibited by simulations
1827   employing lattice summation techniques.  The cumulative energy
1828 < difference results showed the undamped {\sc sf} and moderately damped
1829 < {\sc sp} methods produced results nearly identical to the Ewald
1828 > difference results showed that the undamped {\sc sf} and moderately
1829 > damped {\sc sp} methods produce results nearly identical to the Ewald
1830   summation.  Similarly for the dynamic features, the undamped or
1831   moderately damped {\sc sf} and moderately damped {\sc sp} methods
1832   produce force and torque vector magnitude and directions very similar
# Line 1831 | Line 1839 | easily parallelizable.  This should result in substant
1839   As in all purely-pairwise cutoff methods, these methods are expected
1840   to scale approximately {\it linearly} with system size, and they are
1841   easily parallelizable.  This should result in substantial reductions
1842 < in the computational cost of performing large simulations.
1842 > in the computational cost associated with large-scale simulations.
1843  
1844   Aside from the computational cost benefit, these techniques have
1845   applicability in situations where the use of the Ewald sum can prove
# Line 1850 | Line 1858 | method with an electrostatic damping of 0.2~\AA$^{-1}$
1858   systems containing point charges, most structural features will be
1859   accurately captured using the undamped {\sc sf} method or the {\sc sp}
1860   method with an electrostatic damping of 0.2~\AA$^{-1}$.  These methods
1861 < would also be appropriate for molecular dynamics simulations where the
1862 < data of interest is either structural or short-time dynamical
1861 > would also be appropriate in molecular dynamics simulations where the
1862 > data of interest are either structural or short-time dynamical
1863   quantities.  For long-time dynamics and collective motions, the safest
1864   pairwise method we have evaluated is the {\sc sf} method with an
1865   electrostatic damping between 0.2 and 0.25~\AA$^{-1}$. It is also
# Line 1860 | Line 1868 | $R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease
1868   $R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc
1869   sf} method should use an $\alpha$ of 0.2175~\AA$^{-1}$ for an
1870   $R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease by
1871 < 0.025~\AA$^{-1}$ for every 1~\AA\ increase in cutoff radius.
1871 > 0.025~\AA$^{-1}$ for every 1~\AA\ increase in the cutoff radius.
1872  
1873   We are not suggesting that there is any flaw with the Ewald sum; in
1874 < fact, it is the standard by which these simple pairwise sums have been
1875 < judged.  However, these results do suggest that in the typical
1874 > fact, it is the standard by which these simple pairwise methods have
1875 > been judged.  However, these results do suggest that in the typical
1876   simulations performed today, the Ewald summation may no longer be
1877   required to obtain the level of accuracy most researchers have come to
1878   expect.

Diff Legend

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