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 3022 by chrisfen, Thu Sep 21 23:21:37 2006 UTC vs.
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES}
2  
3 < In molecular simulations, proper accumulation of the electrostatic
3 > In molecular simulations, proper accumulation of electrostatic
4   interactions is essential and is 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 101 | Line 101 | system is said to be using conducting (or ``tin-foil''
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}
104 > 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 481 | Line 481 | 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 496 | Line 505 | Fig. \ref{fig:linearFit}.
505   Sample correlation plots for two alternate methods are shown in
506   Fig. \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
# Line 699 | Line 698 | inverted triangles).}
698   inverted triangles).}
699   \label{fig:delE}
700   \end{figure}
702 –
701   The most striking feature of this plot is how well the Shifted Force
702   ({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
703   differences.  For the undamped {\sc sf} method, and the
# Line 721 | Line 719 | shows an excellent correlation and quality of fit with
719   For the {\sc sp} method, inclusion of electrostatic damping improves
720   the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$
721   shows an excellent correlation and quality of fit with the {\sc spme}
722 < results, particularly with a cutoff radius greater than 12~\AA\.  Use
722 > results, particularly with a cutoff radius greater than 12~\AA . Use
723   of a larger damping parameter is more helpful for the shortest cutoff
724   shown, but it has a detrimental effect on simulations with larger
725   cutoffs.
# Line 759 | Line 757 | inverted triangles).}
757   inverted triangles).}
758   \label{fig:frcMag}
759   \end{figure}
762 –
760   Again, it is striking how well the Shifted Potential and Shifted Force
761   methods are doing at reproducing the {\sc spme} forces.  The undamped and
762   weakly-damped {\sc sf} method gives the best agreement with Ewald.
# Line 798 | Line 795 | inverted triangles).}
795   inverted triangles).}
796   \label{fig:trqMag}
797   \end{figure}
801 –
798   Molecular torques were only available from the systems which contained
799   rigid molecules (i.e. the systems containing water).  The data in
800   fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
# Line 807 | Line 803 | cutoff radius.  Again, the weakly damped and undamped
803   distance.   The striking feature in comparing the new electrostatic
804   methods with {\sc spme} is how much the agreement improves with increasing
805   cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
806 < appears to be reproducing the {\sc spme} torques most accurately.  
806 > appears to reproduce the {\sc spme} torques most accurately.  
807  
808   Water molecules are dipolar, and the reaction field method reproduces
809   the effect of the surrounding polarized medium on each of the
# Line 816 | Line 812 | performs best of all of the methods on molecular torqu
812  
813   \section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}
814  
815 < It is clearly important that a new electrostatic method can reproduce
816 < the magnitudes of the force and torque vectors obtained via the Ewald
817 < sum. However, the {\it directionality} of these vectors will also be
818 < vital in calculating dynamical quantities accurately.  Force and
819 < torque directionalities were investigated by measuring the angles
820 < formed between these vectors and the same vectors calculated using
821 < {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
822 < variance ($\sigma^2$) of the Gaussian fits of the angle error
823 < distributions of the combined set over all system types.
815 > It is clearly important that a new electrostatic method should be able
816 > to reproduce the magnitudes of the force and torque vectors obtained
817 > via the Ewald sum. However, the {\it directionality} of these vectors
818 > will also be vital in calculating dynamical quantities accurately.
819 > Force and torque directionalities were investigated by measuring the
820 > angles formed between these vectors and the same vectors calculated
821 > using {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared
822 > through the variance ($\sigma^2$) of the Gaussian fits of the angle
823 > error distributions of the combined set over all system types.
824  
825   \begin{figure}
826   \centering
# Line 839 | Line 835 | and 15~\AA\ = inverted triangles).}
835   and 15~\AA\ = inverted triangles).}
836   \label{fig:frcTrqAng}
837   \end{figure}
842 –
838   Both the force and torque $\sigma^2$ results from the analysis of the
839   total accumulated system data are tabulated in figure
840   \ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
# Line 977 | Line 972 | are stiffer than the moderately damped and {\sc spme}
972   are stiffer than the moderately damped and {\sc spme} methods.}
973   \label{fig:vCorrPlot}
974   \end{figure}
980 –
975   The short-time decay of the velocity autocorrelation function through
976   the first collision are nearly identical in figure
977   \ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
# Line 996 | Line 990 | important.
990  
991   \section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}
992  
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 –
993   \begin{figure}
994   \centering
995   \includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
# Line 1015 | Line 1000 | functions of NaCl crystals at 1000~K while using {\sc
1000   100~cm$^{-1}$ to highlight where the spectra differ.}
1001   \label{fig:methodPS}
1002   \end{figure}
1003 + To evaluate how the differences between the methods affect the
1004 + collective long-time motion, we computed power spectra from long-time
1005 + traces of the velocity autocorrelation function. The power spectra for
1006 + the best-performing alternative methods are shown in
1007 + fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1008 + a cubic switching function between 40 and 50~ps was used to reduce the
1009 + ringing resulting from data truncation.  This procedure had no
1010 + noticeable effect on peak location or magnitude.
1011  
1012   While the high frequency regions of the power spectra for the
1013   alternative methods are quantitatively identical with Ewald spectrum,
1014   the low frequency region shows how the summation methods differ.
1015   Considering the low-frequency inset (expanded in the upper frame of
1016 < figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the
1016 > figure \ref{fig:methodPS}), at frequencies below 100~cm$^{-1}$, the
1017   correlated motions are blue-shifted when using undamped or weakly
1018   damped {\sc sf}.  When using moderate damping ($\alpha =
1019   0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly
# Line 1030 | Line 1023 | moderately damped methods than for undamped or weakly
1023   long-ranged correlated motions are at lower frequencies for the
1024   moderately damped methods than for undamped or weakly damped methods.
1025  
1026 + \begin{figure}
1027 + \centering
1028 + \includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
1029 + \caption{Effect of damping on the two lowest-frequency phonon modes in
1030 + the NaCl crystal at 1000~K.  The undamped shifted force ({\sc sf})
1031 + method is off by less than 10~cm$^{-1}$, and increasing the
1032 + electrostatic damping to 0.25~\AA$^{-1}$ gives quantitative agreement
1033 + with the power spectrum obtained using the Ewald sum.  Over-damping can
1034 + result in underestimates of frequencies of the long-wavelength
1035 + motions.}
1036 + \label{fig:dampInc}
1037 + \end{figure}
1038   To isolate the role of the damping constant, we have computed the
1039   spectra for a single method ({\sc sf}) with a range of damping
1040   constants and compared this with the {\sc spme} spectrum.
# Line 1044 | Line 1049 | cutoff distance.
1049   obtained using moderate damping in addition to the shifting at the
1050   cutoff distance.
1051  
1047 – \begin{figure}
1048 – \centering
1049 – \includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
1050 – \caption{Effect of damping on the two lowest-frequency phonon modes in
1051 – the NaCl crystal at 1000~K.  The undamped shifted force ({\sc sf})
1052 – method is off by less than 10~cm$^{-1}$, and increasing the
1053 – electrostatic damping to 0.25~\AA$^{-1}$ gives quantitative agreement
1054 – with the power spectrum obtained using the Ewald sum.  Over-damping can
1055 – result in underestimates of frequencies of the long-wavelength
1056 – motions.}
1057 – \label{fig:dampInc}
1058 – \end{figure}
1059 –
1052   \section{An Application: TIP5P-E Water}\label{sec:t5peApplied}
1053  
1054   The above sections focused on the energetics and dynamics of a variety
# Line 1114 | Line 1106 | p(t) =         \frac{1}{\textrm{d}V}\sum_\mu
1106          + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
1107   \label{eq:MolecularPressure}
1108   \end{equation}
1109 < where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of
1110 < molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass
1111 < ($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on
1112 < atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the
1113 < right term in the brackets of equation \ref{eq:MolecularPressure}) is
1114 < directly dependent on the interatomic forces.  Since the {\sc sp}
1115 < method does not modify the forces (see
1116 < section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will
1117 < be identical to that obtained without an electrostatic correction.
1118 < The {\sc sf} method does alter the virial component and, by way of the
1119 < modified pressures, should provide densities more in line with those
1120 < obtained using the Ewald summation.
1109 > where d is the dimensionality of the system, $V$ is the volume,
1110 > $\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
1111 > is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
1112 > $\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
1113 > $\mu$.\cite{Melchionna93} The virial term (the right term in the
1114 > brackets of equation
1115 > \ref{eq:MolecularPressure}) is directly dependent on the interatomic
1116 > forces.  Since the {\sc sp} method does not modify the forces (see
1117 > section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
1118 > will be identical to that obtained without an electrostatic
1119 > correction.  The {\sc sf} method does alter the virial component and,
1120 > by way of the modified pressures, should provide densities more in
1121 > line with those obtained using the Ewald summation.
1122  
1123   To compare densities, $NPT$ simulations were performed with the same
1124   temperatures as those selected by Rick in his Ewald summation
# Line 1147 | Line 1140 | greater force on the central particle. The error bars
1140   Ewald summation, leading to slightly lower densities. This effect is
1141   more visible with the 9~\AA\ cutoff, where the image charges exert a
1142   greater force on the central particle. The error bars for the {\sc sf}
1143 < methods show plus or minus the standard deviation of the density
1144 < measurement at each temperature.}
1143 > methods show the average one-sigma uncertainty of the density
1144 > measurement, and this uncertainty is the same for all the {\sc sf}
1145 > curves.}
1146   \label{fig:t5peDensities}
1147   \end{figure}
1154 –
1148   Figure \ref{fig:t5peDensities} shows the densities calculated for
1149   TIP5P-E using differing electrostatic corrections overlaid on the
1150   experimental values.\cite{CRC80} The densities when using the {\sc sf}
# 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 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 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.
1234  
1235   \subsection{Thermodynamic Properties}\label{sec:t5peThermo}
1236  
# Line 1368 | Line 1360 | property trends with temperature seen when using the E
1360  
1361   As observed for the density in section \ref{sec:t5peDensity}, the
1362   property trends with temperature seen when using the Ewald summation
1363 < are reproduced with the {\sc sf} technique. Differences include the
1364 < calculated values of $\Delta H_\textrm{vap}$ under-predicting the Ewald
1365 < values. This is to be expected due to the direct weakening of the
1366 < electrostatic interaction through forced neutralization in {\sc
1367 < sf}. This results in an increase of the intermolecular potential
1368 < producing lower values from equation (\ref{eq:DeltaHVap}). The slopes of
1369 < these values with temperature are similar to that seen using the Ewald
1370 < summation; however, they are both steeper than the experimental trend,
1371 < indirectly resulting in the inflated $C_p$ values at all temperatures.
1363 > are reproduced with the {\sc sf} technique. One noticable difference
1364 > between the properties calculated using the two methods are the lower
1365 > $\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
1366 > expected due to the direct weakening of the electrostatic interaction
1367 > through forced neutralization. This results in an increase of the
1368 > intermolecular potential producing lower values from equation
1369 > (\ref{eq:DeltaHVap}). The slopes of these values with temperature are
1370 > similar to that seen using the Ewald summation; however, they are both
1371 > steeper than the experimental trend, indirectly resulting in the
1372 > inflated $C_p$ values at all temperatures.
1373  
1374   Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$
1375   values all overlap within error. As indicated for the $\Delta
# Line 1398 | Line 1391 | the onset of a more frustrated or glassy behavior for
1391   indicate a more pronounced transition in the supercooled regime,
1392   particularly in the case of {\sc sf} without damping. This points to
1393   the onset of a more frustrated or glassy behavior for TIP5P-E at
1394 < temperatures below 250~K in these simulations. Because the systems are
1395 < locked in different regions of phase-space, comparisons between
1396 < properties at these temperatures are not exactly fair. This
1397 < observation is explored in more detail in section
1398 < \ref{sec:t5peDynamics}.
1394 > temperatures below 250~K in the {\sc sf} simulations, indicating that
1395 > disorder in the reciprical-space term of the Ewald summation might act
1396 > to loosen up the local structure more than the image-charges in {\sc
1397 > sf}. Because the systems are locked in different regions of
1398 > phase-space, comparisons between properties at these temperatures are
1399 > not exactly fair. This observation is explored in more detail in
1400 > section \ref{sec:t5peDynamics}.
1401  
1402   The final thermodynamic property displayed in figure
1403   \ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
# Line 1435 | Line 1430 | relation using the mean square displacement (MSD),
1430   self-diffusion constants ($D$) were calculated with the Einstein
1431   relation using the mean square displacement (MSD),
1432   \begin{equation}
1433 < D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1433 > D = \lim_{t\rightarrow\infty}
1434 >    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1435   \label{eq:MSD}
1436   \end{equation}
1437   where $t$ is time, and $\mathbf{r}_i$ is the position of particle
# Line 1446 | Line 1442 | regions:
1442   \begin{enumerate}[itemsep=0pt]
1443   \item parabolic short-time ballistic motion,
1444   \item linear diffusive regime, and
1445 < \item poor statistic region at long-time.
1445 > \item a region with poor statistics.
1446   \end{enumerate}
1447   The slope from the linear region (region 2) is used to calculate $D$.
1448   \begin{figure}
# Line 1515 | Line 1511 | easier comparisons in the more relevant temperature re
1511   easier comparisons in the more relevant temperature regime.}
1512   \label{fig:t5peDynamics}
1513   \end{figure}
1514 < Results for the diffusion constants and reorientational time constants
1514 > Results for the diffusion constants and orientational relaxation times
1515   are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
1516   apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
1517   the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
# Line 1524 | Line 1520 | diffuse a little faster than with the Ewald sum; howev
1520   insight into differences between the electrostatic summation
1521   techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
1522   diffuse a little faster than with the Ewald sum; however, use of light
1523 < to moderate damping results in indistinguishable $D$ values. Though not
1524 < apparent in this figure, {\sc sf} values at the lowest temperature are
1525 < approximately an order of magnitude lower than with Ewald. These
1523 > to moderate damping results in indistinguishable $D$ values. Though
1524 > not apparent in this figure, {\sc sf} values at the lowest temperature
1525 > are approximately an order of magnitude lower than with Ewald. These
1526   values support the observation from section \ref{sec:t5peThermo} that
1527   there appeared to be a change to a more glassy-like phase with the
1528   {\sc sf} technique at these lower temperatures.
# Line 1540 | Line 1536 | calculation from a 10~ps trajectory with {\sc sf} with
1536   for this deviation between techniques. The Ewald results were taken
1537   shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick
1538   calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of
1539 < 0.2~\AA$^-1$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, placing
1540 < the result more in line with that obtained using the Ewald sum. These
1541 < results support this explanation; however, recomputing the results to
1542 < meet a poorer statistical standard is counter-productive. Assuming the
1543 < Ewald results are not the product of poor statistics, differences in
1544 < techniques to integrate the orientational motion could also play a
1545 < role. {\sc shake} is the most commonly used technique for
1546 < approximating rigid-body orientational motion,\cite{Ryckaert77} where
1547 < as in {\sc oopse}, we maintain and integrate the entire rotation
1548 < matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake}
1549 < is an iterative constraint technique, if the convergence tolerances
1550 < are raised for increased performance, error will accumulate in the
1551 < orientational motion. Finally, the Ewald results were calculated using
1552 < the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf}
1539 > 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$,
1540 > placing the result more in line with that obtained using the Ewald
1541 > sum. These results support this explanation; however, recomputing the
1542 > results to meet a poorer statistical standard is
1543 > counter-productive. Assuming the Ewald results are not the product of
1544 > poor statistics, differences in techniques to integrate the
1545 > orientational motion could also play a role. {\sc shake} is the most
1546 > commonly used technique for approximating rigid-body orientational
1547 > motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and
1548 > integrate the entire rotation matrix using the {\sc dlm}
1549 > method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
1550 > technique, if the convergence tolerances are raised for increased
1551 > performance, error will accumulate in the orientational
1552 > motion. Finally, the Ewald results were calculated using the $NVT$
1553 > ensemble, while the $NVE$ ensemble was used for {\sc sf}
1554   calculations. The additional mode of motion due to the thermostat will
1555   alter the dynamics, resulting in differences between $NVT$ and $NVE$
1556   results. These differences are increasingly noticeable as the
# Line 1565 | Line 1562 | consider an extension of these techniques to include p
1562   neutralizing the cutoff sphere with charge-charge interaction shifting
1563   and by damping the electrostatic interactions. Now we would like to
1564   consider an extension of these techniques to include point multipole
1565 < interactions. How will the shifting and damping need to develop in
1565 > interactions. How will the shifting and damping need to be modified in
1566   order to accommodate point multipoles?
1567  
1568 < Of the two techniques, the least to vary is shifting. Shifting is
1568 > Of the two techniques, the easiest to adapt is shifting. Shifting is
1569   employed to neutralize the cutoff sphere; however, in a system
1570   composed purely of point multipoles, the cutoff sphere is already
1571   neutralized. This means that shifting is not necessary between point
# Line 1586 | Line 1583 | than considering only the interactions between single
1583   replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole
1584   expansion.\cite{Hirschfelder67} In the multipole expansion, rather
1585   than considering only the interactions between single point charges,
1586 < the electrostatic interactions is reformulated such that it describes
1586 > the electrostatic interaction is reformulated such that it describes
1587   the interaction between charge distributions about central sites of
1588   the respective sets of charges. This procedure is what leads to the
1589   familiar charge-dipole,

Diff Legend

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