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 3016 by chrisfen, Thu Sep 21 23:21:37 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 the electrostatic
4 < interactions is essential and is one of the most
3 > In molecular simulations, proper accumulation of electrostatic
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
8 < every pair interaction involves a calculation of charge-charge forces.
7 > by repulsive Lennard-Jones interactions.  This means that nearly every
8 > pair interaction involves a calculation of charge-charge forces.
9   Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
10   interactions quickly become the most expensive part of molecular
11   simulations.  Historically, the electrostatic pair interaction would
# 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$. Figure
105 < \ref{fig:ewaldTime} shows how the Ewald sum has been applied over
106 < time.  Initially, due to the small system sizes that could be
107 < simulated feasibly, the entire simulation box was replicated to
108 < convergence.  In more modern simulations, the systems have grown large
109 < enough that a real-space cutoff could potentially give convergent
110 < behavior.  Indeed, it has been observed that with the choice of a
111 < small $\alpha$, the reciprocal-space portion of the Ewald sum can be
112 < rapidly convergent and small relative to the real-space
113 < portion.\cite{Karasawa89,Kolafa92}
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 121 | Line 112 | convergence for the larger systems of charges that are
112   convergence for the larger systems of charges that are common today.}
113   \label{fig:ewaldTime}
114   \end{figure}
115 + Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied
116 + over time.  Initially, due to the small system sizes that could be
117 + simulated feasibly, the entire simulation box was replicated to
118 + convergence.  In more modern simulations, the systems have grown large
119 + enough that a real-space cutoff could potentially give convergent
120 + behavior.  Indeed, it has been observed that with the choice of a
121 + small $\alpha$, the reciprocal-space portion of the Ewald sum can be
122 + rapidly convergent and small relative to the real-space
123 + portion.\cite{Karasawa89,Kolafa92}
124  
125   The original Ewald summation is an $\mathcal{O}(N^2)$ algorithm.  The
126   convergence parameter $(\alpha)$ plays an important role in balancing
# 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}
483  
484 + \begin{figure}
485 + \centering
486 + \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
487 + \caption{Example least squares regressions of the configuration energy
488 + differences for SPC/E water systems. The upper plot shows a data set
489 + with a poor correlation coefficient ($R^2$), while the lower plot
490 + shows a data set with a good correlation coefficient.}
491 + \label{fig:linearFit}
492 + \end{figure}
493   The pairwise summation techniques (outlined in section
494   \ref{sec:ESMethods}) were evaluated for use in MC simulations by
495   studying the energy differences between conformations.  We took the
# Line 494 | 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  
499 – \begin{figure}
500 – \centering
501 – \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
502 – \caption{Example least squares regressions of the configuration energy
503 – differences for SPC/E water systems. The upper plot shows a data set
504 – with a poor correlation coefficient ($R^2$), while the lower plot
505 – shows a data set with a good correlation coefficient.}
506 – \label{fig:linearFit}
507 – \end{figure}
508 –
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
517 < cumulative results over all the investigated systems appear below in
518 < 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 528 | 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 547 | 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 610 | 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}
613 –
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 699 | Line 695 | inverted triangles).}
695   inverted triangles).}
696   \label{fig:delE}
697   \end{figure}
702 –
698   The most striking feature of this plot is how well the Shifted Force
699   ({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
700   differences.  For the undamped {\sc sf} method, and the
# Line 721 | Line 716 | shows an excellent correlation and quality of fit with
716   For the {\sc sp} method, inclusion of electrostatic damping improves
717   the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$
718   shows an excellent correlation and quality of fit with the {\sc spme}
719 < results, particularly with a cutoff radius greater than 12~\AA\.  Use
719 > results, particularly with a cutoff radius greater than 12~\AA . Use
720   of a larger damping parameter is more helpful for the shortest cutoff
721   shown, but it has a detrimental effect on simulations with larger
722   cutoffs.
# Line 735 | 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 759 | Line 754 | inverted triangles).}
754   inverted triangles).}
755   \label{fig:frcMag}
756   \end{figure}
757 <
758 < Again, it is striking how well the Shifted Potential and Shifted Force
759 < methods are doing at reproducing the {\sc spme} forces.  The undamped and
760 < weakly-damped {\sc sf} method gives the best agreement with Ewald.
761 < This is perhaps expected because this method explicitly incorporates a
767 < 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 775 | 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 798 | Line 792 | inverted triangles).}
792   inverted triangles).}
793   \label{fig:trqMag}
794   \end{figure}
801 –
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 be reproducing 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 816 | Line 810 | performs best of all of the methods on molecular torqu
810  
811   \section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}
812  
813 < It is clearly important that a new electrostatic method can reproduce
814 < the magnitudes of the force and torque vectors obtained via the Ewald
815 < sum. However, the {\it directionality} of these vectors will also be
816 < vital in calculating dynamical quantities accurately.  Force and
817 < torque directionalities were investigated by measuring the angles
818 < formed between these vectors and the same vectors calculated using
819 < {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
820 < variance ($\sigma^2$) of the Gaussian fits of the angle error
821 < distributions of the combined set over all system types.
813 > It is clearly important that a new electrostatic method should be able
814 > to reproduce the magnitudes of the force and torque vectors obtained
815 > via the Ewald sum. However, the {\it directionality} of these vectors
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 (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  
823   \begin{figure}
824   \centering
# Line 839 | Line 833 | and 15~\AA\ = inverted triangles).}
833   and 15~\AA\ = inverted triangles).}
834   \label{fig:frcTrqAng}
835   \end{figure}
842 –
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)
850 < 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 929 | 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 948 | 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 977 | 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 <
981 < 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 996 | Line 988 | important.
988  
989   \section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}
990  
999 – To evaluate how the differences between the methods affect the
1000 – collective long-time motion, we computed power spectra from long-time
1001 – traces of the velocity autocorrelation function. The power spectra for
1002 – the best-performing alternative methods are shown in
1003 – fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1004 – a cubic switching function between 40 and 50~ps was used to reduce the
1005 – ringing resulting from data truncation.  This procedure had no
1006 – noticeable effect on peak location or magnitude.
1007 –
991   \begin{figure}
992   \centering
993   \includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
# Line 1015 | Line 998 | functions of NaCl crystals at 1000~K while using {\sc
998   100~cm$^{-1}$ to highlight where the spectra differ.}
999   \label{fig:methodPS}
1000   \end{figure}
1001 + To evaluate how the differences between the methods affect the
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 + 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.
1009  
1010   While the high frequency regions of the power spectra for the
1011   alternative methods are quantitatively identical with Ewald spectrum,
1012   the low frequency region shows how the summation methods differ.
1013   Considering the low-frequency inset (expanded in the upper frame of
1014 < figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the
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
1022   moderately damped methods than for undamped or weakly damped methods.
1023  
1033 – To isolate the role of the damping constant, we have computed the
1034 – spectra for a single method ({\sc sf}) with a range of damping
1035 – constants and compared this with the {\sc spme} spectrum.
1036 – Fig. \ref{fig:dampInc} shows more clearly that increasing the
1037 – electrostatic damping red-shifts the lowest frequency phonon modes.
1038 – However, even without any electrostatic damping, the {\sc sf} method
1039 – has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
1040 – Without the {\sc sf} modifications, an undamped (pure cutoff) method
1041 – would predict the lowest frequency peak near 325~cm$^{-1}$.  {\it
1042 – Most} of the collective behavior in the crystal is accurately captured
1043 – using the {\sc sf} method.  Quantitative agreement with Ewald can be
1044 – obtained using moderate damping in addition to the shifting at the
1045 – cutoff distance.
1046 –
1024   \begin{figure}
1025   \centering
1026   \includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
# Line 1056 | Line 1033 | motions.}
1033   motions.}
1034   \label{fig:dampInc}
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.  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 1076 | 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
1092 < 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 1104 | p(t) =         \frac{1}{\textrm{d}V}\sum_\mu
1104          + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
1105   \label{eq:MolecularPressure}
1106   \end{equation}
1107 < where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of
1108 < molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass
1109 < ($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on
1110 < atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the
1111 < right term in the brackets of equation \ref{eq:MolecularPressure}) is
1112 < directly dependent on the interatomic forces.  Since the {\sc sp}
1113 < method does not modify the forces (see
1114 < section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will
1115 < be identical to that obtained without an electrostatic correction.
1116 < The {\sc sf} method does alter the virial component and, by way of the
1117 < modified pressures, should provide densities more in line with those
1118 < obtained using the Ewald summation.
1107 > where d is the dimensionality of the system, $V$ is the volume,
1108 > $\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
1109 > is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
1110 > $\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
1111 > $\mu$.\cite{Melchionna93} The virial term (the right term in 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}
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
1119 > line with those obtained using the Ewald summation.
1120  
1121   To compare densities, $NPT$ simulations were performed with the same
1122   temperatures as those selected by Rick in his Ewald summation
# Line 1135 | 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 1147 | Line 1139 | greater force on the central particle. The error bars
1139   Ewald summation, leading to slightly lower densities. This effect is
1140   more visible with the 9~\AA\ cutoff, where the image charges exert a
1141   greater force on the central particle. The error bars for the {\sc sf}
1142 < methods show plus or minus the standard deviation of the density
1143 < measurement at each temperature.}
1142 > methods show the average one-sigma uncertainty of the density
1143 > measurement, and this uncertainty is the same for all the {\sc sf}
1144 > curves.}
1145   \label{fig:t5peDensities}
1146   \end{figure}
1154 –
1147   Figure \ref{fig:t5peDensities} shows the densities calculated for
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 1176 | Line 1169 | As a final note, all of the above density calculations
1169   important role in the resulting densities.
1170  
1171   As a final note, all of the above density calculations were performed
1172 < with systems of 512 water molecules. Rick observed a system sized
1172 > with systems of 512 water molecules. Rick observed a system size
1173   dependence of the computed densities when using the Ewald summation,
1174   most likely due to his tying of the convergence parameter to the box
1175   dimensions.\cite{Rick04} For systems of 256 water molecules, the
# Line 1218 | 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 1231 | Line 1224 | identical.}
1224   identical.}
1225   \label{fig:t5peGofRs}
1226   \end{figure}
1234 –
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. The differences in
1230 < density do not appear to have any effect on the liquid structure as
1231 < the $g_\textrm{OO}(r)$s are indistinguishable. These results indicate
1232 < that the $g_\textrm{OO}(r)$ is insensitive to the choice of
1233 < 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 1253 | 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 1368 | 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. Differences include the
1365 < calculated values of $\Delta H_\textrm{vap}$ under-predicting the Ewald
1366 < values. This is to be expected due to the direct weakening of the
1367 < electrostatic interaction through forced neutralization in {\sc
1368 < sf}. This results in an increase of the intermolecular potential
1369 < producing lower values from equation (\ref{eq:DeltaHVap}). The slopes of
1370 < these values with temperature are similar to that seen using the Ewald
1371 < summation; however, they are both steeper than the experimental trend,
1372 < indirectly resulting in the inflated $C_p$ values at all temperatures.
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
1368 > through forced neutralization. This results in an increase of the
1369 > intermolecular potential producing lower values from equation
1370 > (\ref{eq:DeltaHVap}). The slopes of these values with temperature are
1371 > similar to that seen using the Ewald summation; however, they are both
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 1398 | Line 1394 | the onset of a more frustrated or glassy behavior for
1394   indicate a more pronounced transition in the supercooled regime,
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 these simulations. Because the systems are
1398 < locked in different regions of phase-space, comparisons between
1399 < properties at these temperatures are not exactly fair. This
1400 < observation is explored in more detail in section
1401 < \ref{sec:t5peDynamics}.
1397 > temperatures below 250~K in the {\sc sf} simulations, indicating that
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}. 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 1412 | 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 1424 | 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 = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1440 > D = \lim_{t\rightarrow\infty}
1441 >    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1442   \label{eq:MSD}
1443   \end{equation}
1444   where $t$ is time, and $\mathbf{r}_i$ is the position of particle
# Line 1446 | Line 1449 | regions:
1449   \begin{enumerate}[itemsep=0pt]
1450   \item parabolic short-time ballistic motion,
1451   \item linear diffusive regime, and
1452 < \item poor statistic region at long-time.
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 1467 | 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 1515 | Line 1519 | easier comparisons in the more relevant temperature re
1519   easier comparisons in the more relevant temperature regime.}
1520   \label{fig:t5peDynamics}
1521   \end{figure}
1522 < Results for the diffusion constants and reorientational time constants
1522 > Results for the diffusion constants and orientational relaxation times
1523   are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
1524   apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
1525   the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
# Line 1524 | Line 1528 | diffuse a little faster than with the Ewald sum; howev
1528   insight into differences between the electrostatic summation
1529   techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
1530   diffuse a little faster than with the Ewald sum; however, use of light
1531 < to moderate damping results in indistinguishable $D$ values. Though not
1532 < apparent in this figure, {\sc sf} values at the lowest temperature are
1533 < approximately an order of magnitude lower than with Ewald. These
1531 > to moderate damping results in indistinguishable $D$ values. Though
1532 > not apparent in this figure, {\sc sf} values at the lowest temperature
1533 > are approximately an order of magnitude lower than with Ewald. These
1534   values support the observation from section \ref{sec:t5peThermo} that
1535   there appeared to be a change to a more glassy-like phase with the
1536   {\sc sf} technique at these lower temperatures.
# Line 1537 | 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$, placing
1548 < the result more in line with that obtained using the Ewald sum. These
1549 < results support this explanation; however, recomputing the results to
1550 < meet a poorer statistical standard is counter-productive. Assuming the
1551 < Ewald results are not the product of poor statistics, differences in
1552 < techniques to integrate the orientational motion could also play a
1553 < role. {\sc shake} is the most commonly used technique for
1554 < approximating rigid-body orientational motion,\cite{Ryckaert77} where
1555 < as in {\sc oopse}, we maintain and integrate the entire rotation
1556 < matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake}
1557 < is an iterative constraint technique, if the convergence tolerances
1558 < are raised for increased performance, error will accumulate in the
1559 < orientational motion. Finally, the Ewald results were calculated using
1560 < the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf}
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} 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
1559 > performance, error will accumulate in the orientational
1560 > motion. Finally, the Ewald results were calculated using the $NVT$
1561 > ensemble, while the $NVE$ ensemble was used for {\sc sf}
1562   calculations. The additional mode of motion due to the thermostat will
1563   alter the dynamics, resulting in differences between $NVT$ and $NVE$
1564   results. These differences are increasingly noticeable as the
# Line 1565 | Line 1570 | consider an extension of these techniques to include p
1570   neutralizing the cutoff sphere with charge-charge interaction shifting
1571   and by damping the electrostatic interactions. Now we would like to
1572   consider an extension of these techniques to include point multipole
1573 < interactions. How will the shifting and damping need to develop in
1573 > interactions. How will the shifting and damping need to be modified in
1574   order to accommodate point multipoles?
1575  
1576 < Of the two techniques, the least to vary is shifting. Shifting is
1576 > Of the two techniques, the easiest to adapt is shifting. Shifting is
1577   employed to neutralize the cutoff sphere; however, in a system
1578   composed purely of point multipoles, the cutoff sphere is already
1579   neutralized. This means that shifting is not necessary between point
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 1586 | Line 1591 | than considering only the interactions between single
1591   replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole
1592   expansion.\cite{Hirschfelder67} In the multipole expansion, rather
1593   than considering only the interactions between single point charges,
1594 < the electrostatic interactions is reformulated such that it describes
1594 > the electrostatic interaction is reformulated such that it describes
1595   the interaction between charge distributions about central sites of
1596   the respective sets of charges. This procedure is what leads to the
1597   familiar charge-dipole,
# Line 1669 | 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 1703 | 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 1734 | 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 1747 | 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 1780 | 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 1820 | 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 1834 | 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 1853 | 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 1863 | 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