--- trunk/ssdePaper/nptSSD.tex	2004/02/04 22:42:52	1024
+++ trunk/ssdePaper/nptSSD.tex	2004/02/05 18:42:59	1027
@@ -215,13 +215,13 @@ field acting on dipole $i$ is
 field acting on dipole $i$ is
 \begin{equation}
 \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
-\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij})\  ,
+\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij}),
 \label{rfequation}
 \end{equation}
 where $\mathcal{R}$ is the cavity defined by the cutoff radius
 ($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
 system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole
-moment vector of particle $j$ and $f(r_{ij})$ is a cubic switching
+moment vector of particle $j$, and $f(r_{ij})$ is a cubic switching
 function.\cite{AllenTildesley} The reaction field contribution to the
 total energy by particle $i$ is given by $-\frac{1}{2}{\bf
 \mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf
@@ -251,14 +251,14 @@ symplectic splitting method proposed by Dullweber {\it
 need to be considered.
 
 Integration of the equations of motion was carried out using the
-symplectic splitting method proposed by Dullweber {\it et
-al.}\cite{Dullweber1997} Our reason for selecting this integrator
-centers on poor energy conservation of rigid body dynamics using
-traditional quaternion integration.\cite{Evans77,Evans77b} In typical
-microcanonical ensemble simulations, the energy drift when using
-quaternions was substantially greater than when using the symplectic
-splitting method (fig. \ref{timestep}).  This steady drift in the
-total energy has also been observed by Kol {\it et al.}\cite{Laird97}
+symplectic splitting method proposed by Dullweber, Leimkuhler, and
+McLachlan (DLM).\cite{Dullweber1997} Our reason for selecting this
+integrator centers on poor energy conservation of rigid body dynamics
+using traditional quaternion integration.\cite{Evans77,Evans77b} In
+typical microcanonical ensemble simulations, the energy drift when
+using quaternions was substantially greater than when using the DLM
+method (fig. \ref{timestep}).  This steady drift in the total energy
+has also been observed by Kol {\it et al.}\cite{Laird97}
 
 The key difference in the integration method proposed by Dullweber
 \emph{et al.} is that the entire rotation matrix is propagated from
@@ -267,19 +267,19 @@ The symplectic splitting method allows for Verlet styl
 rotation matrix into quaternions for storage purposes makes trajectory
 data quite compact.
 
-The symplectic splitting method allows for Verlet style integration of
-both translational and orientational motion of rigid bodies. In this
+The DML method allows for Verlet style integration of both
+translational and orientational motion of rigid bodies. In this
 integration method, the orientational propagation involves a sequence
 of matrix evaluations to update the rotation
 matrix.\cite{Dullweber1997} These matrix rotations are more costly
 than the simpler arithmetic quaternion propagation. With the same time
 step, a 1000 SSD particle simulation shows an average 7\% increase in
-computation time using the symplectic step method in place of
-quaternions. The additional expense per step is justified when one
-considers the ability to use time steps that are nearly twice as large
-under symplectic splitting than would be usable under quaternion
-dynamics.  The energy conservation of the two methods using a number
-of different time steps is illustrated in figure
+computation time using the DML method in place of quaternions. The
+additional expense per step is justified when one considers the
+ability to use time steps that are nearly twice as large under DML
+than would be usable under quaternion dynamics.  The energy
+conservation of the two methods using a number of different time steps
+is illustrated in figure
 \ref{timestep}.
 
 \begin{figure}
@@ -287,33 +287,33 @@ the symplectic step method proposed by Dullweber \emph
 \epsfxsize=6in
 \epsfbox{timeStep.epsi}
 \caption{Energy conservation using both quaternion based integration and
-the symplectic step method proposed by Dullweber \emph{et al.} with
-increasing time step. The larger time step plots are shifted from the
-true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.}
+the symplectic splitting method proposed by Dullweber \emph{et al.}
+with increasing time step. The larger time step plots are shifted from
+the true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.}
 \label{timestep}
 \end{center}
 \end{figure}
 
 In figure \ref{timestep}, the resulting energy drift at various time
-steps for both the symplectic step and quaternion integration schemes
-is compared.  All of the 1000 SSD particle simulations started with
-the same configuration, and the only difference was the method used to
-handle orientational motion. At time steps of 0.1 and 0.5 fs, both
-methods for propagating the orientational degrees of freedom conserve
-energy fairly well, with the quaternion method showing a slight energy
-drift over time in the 0.5 fs time step simulation. At time steps of 1
-and 2 fs, the energy conservation benefits of the symplectic step
-method are clearly demonstrated. Thus, while maintaining the same
-degree of energy conservation, one can take considerably longer time
-steps, leading to an overall reduction in computation time.
+steps for both the DML and quaternion integration schemes is compared.
+All of the 1000 SSD particle simulations started with the same
+configuration, and the only difference was the method used to handle
+orientational motion. At time steps of 0.1 and 0.5 fs, both methods
+for propagating the orientational degrees of freedom conserve energy
+fairly well, with the quaternion method showing a slight energy drift
+over time in the 0.5 fs time step simulation. At time steps of 1 and 2
+fs, the energy conservation benefits of the DML method are clearly
+demonstrated. Thus, while maintaining the same degree of energy
+conservation, one can take considerably longer time steps, leading to
+an overall reduction in computation time.
 
-Energy drift in the symplectic step simulations was unnoticeable for
-time steps up to 3 fs. A slight energy drift on the
-order of 0.012 kcal/mol per nanosecond was observed at a time step of
-4 fs, and as expected, this drift increases dramatically
-with increasing time step. To insure accuracy in our microcanonical
-simulations, time steps were set at 2 fs and kept at this value for
-constant pressure simulations as well.
+Energy drift in the simulations using DML integration was unnoticeable
+for time steps up to 3 fs. A slight energy drift on the order of 0.012
+kcal/mol per nanosecond was observed at a time step of 4 fs, and as
+expected, this drift increases dramatically with increasing time
+step. To insure accuracy in our microcanonical simulations, time steps
+were set at 2 fs and kept at this value for constant pressure
+simulations as well.
 
 Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
 were generated as starting points for all simulations. The $I_h$
@@ -430,7 +430,7 @@ was what lead Ichiye {\it et al.} to recently reparame
 simulations.\cite{Clancy94,Jorgensen98b} However, even without the
 reaction field, the density around 300 K is still significantly lower
 than experiment and comparable water models. This anomalous behavior
-was what lead Ichiye {\it et al.} to recently reparameterize
+was what lead Tan {\it et al.} to recently reparameterize
 SSD.\cite{Ichiye03} Throughout the remainder of the paper our
 reparamaterizations of SSD will be compared with the newer SSD1 model.
 
@@ -642,7 +642,7 @@ and broadened the first peak of SSD/E and SSD/RF.}
 \begin{figure}
 \begin{center} 
 \epsfxsize=6in
-\epsfbox{dualsticky.ps} 
+\epsfbox{dualsticky_bw.eps} 
 \caption{Isosurfaces of the sticky potential for SSD1 (left) and SSD/E \& 
 SSD/RF (right). Light areas correspond to the tetrahedral attractive
 component, and darker areas correspond to the dipolar repulsive
@@ -812,23 +812,24 @@ without an active reaction field at the densities calc
 the densities, it is important that the excellent diffusive behavior
 of SSD be maintained or improved. Figure \ref{ssdediffuse} compares
 the temperature dependence of the diffusion constant of SSD/E to SSD1
-without an active reaction field at the densities calculated from the
-NPT simulations at 1 atm. The diffusion constant for SSD/E is
-consistently higher than experiment, while SSD1 remains lower than
-experiment until relatively high temperatures (around 360 K). Both
-models follow the shape of the experimental curve well below 300 K but
-tend to diffuse too rapidly at higher temperatures, as seen in SSD1's
-crossing above 360 K.  This increasing diffusion relative to the
-experimental values is caused by the rapidly decreasing system density
-with increasing temperature.  Both SSD1 and SSD/E show this deviation
-in diffusive behavior, but this trend has different implications on
-the diffusive behavior of the models.  While SSD1 shows more
-experimentally accurate diffusive behavior in the high temperature
-regimes, SSD/E shows more accurate behavior in the supercooled and
-biologically relevant temperature ranges.  Thus, the changes made to
-improve the liquid structure may have had an adverse affect on the
-density maximum, but they improve the transport behavior of SSD/E
-relative to SSD1 under the most commonly simulated conditions.
+without an active reaction field at the densities calculated from
+their respective NPT simulations at 1 atm. The diffusion constant for
+SSD/E is consistently higher than experiment, while SSD1 remains lower
+than experiment until relatively high temperatures (around 360
+K). Both models follow the shape of the experimental curve well below
+300 K but tend to diffuse too rapidly at higher temperatures, as seen
+in SSD1's crossing above 360 K.  This increasing diffusion relative to
+the experimental values is caused by the rapidly decreasing system
+density with increasing temperature.  Both SSD1 and SSD/E show this
+deviation in particle mobility, but this trend has different
+implications on the diffusive behavior of the models.  While SSD1
+shows more experimentally accurate diffusive behavior in the high
+temperature regimes, SSD/E shows more accurate behavior in the
+supercooled and biologically relevant temperature ranges.  Thus, the
+changes made to improve the liquid structure may have had an adverse
+affect on the density maximum, but they improve the transport behavior
+of SSD/E relative to SSD1 under the most commonly simulated
+conditions.
 
 \begin{figure}
 \begin{center} 
@@ -862,7 +863,7 @@ reparameterization when using an altered long-range co
 
 \begin{table}
 \begin{center}
-\caption{Calculated and experimental properties of the single point waters and liquid water at 298 K and 1 atm. (a) Ref. [\citen{Mills73}]. (b) Calculated by integrating the data in ref. \citen{Head-Gordon00_1}. (c) Calculated by integrating the data in ref. \citen{Soper86}. (d) Ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.}
+\caption{Calculated and experimental properties of the single point waters and liquid water at 298 K and 1 atm. (a) Ref. [\citen{Mills73}]. (b) Calculated by integrating the data in ref. \citen{Head-Gordon00_1}. (c) Calculated by integrating the data in ref. \citen{Soper86}. (d) Calculated for 298 K from data in ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.}
 \begin{tabular}{ l  c  c  c  c  c } 
 \hline \\[-3mm]
 \ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ 
@@ -872,9 +873,9 @@ reparameterization when using an altered long-range co
 \ \ \ $C_p$ (cal/mol K) & 28.80 $\pm$0.11 & 25.45 $\pm$0.09 & 28.28 $\pm$0.06 & 23.83 $\pm$0.16 & 17.98 \\
 \ \ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78 $\pm$0.07 & 2.51 $\pm$0.18 & 2.00 $\pm$0.17 & 2.32 $\pm$0.06 & 2.299$^\text{a}$ \\
 \ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\
-\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.4$^\text{c}$ \\
-\ \ \ $\tau_1^\mu$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 4.76$^\text{d}$ \\
-\ \ \ $\tau_2^\mu$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2 $\pm$0.2 & 2.3$^\text{e}$ \\
+\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.5$^\text{c}$ \\
+\ \ \ $\tau_1$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 5.7$^\text{d}$ \\
+\ \ \ $\tau_2$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2 $\pm$0.2 & 2.3$^\text{e}$ \\
 \end{tabular}
 \label{liquidproperties}
 \end{center}
@@ -883,15 +884,17 @@ coordination number and hydrogen bonds per particle we
 Table \ref{liquidproperties} gives a synopsis of the liquid state
 properties of the water models compared in this study along with the
 experimental values for liquid water at ambient conditions. The
-coordination number and hydrogen bonds per particle were calculated by
-integrating the following relation:
+coordination number ($N_C$) and hydrogen bonds per particle ($N_H$)
+were calculated by integrating the following relations:
 \begin{equation}
-4\pi\rho\int_{0}^{a}r^2\text{g}(r)dr,
+N_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr,
 \end{equation}
-where $\rho$ is the number density of pair interactions, $a$ is the
-radial location of the minima following the first solvation shell
-peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for
-calculation of the coordination number or hydrogen bonds per particle
+\begin{equation}
+N_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr,
+\end{equation}
+where $\rho$ is the number density of the specified pair interactions,
+$a$ and $b$ are the radial locations of the minima following the first
+solvation shell peak in g$_\text{OO}(r)$ or g$_\text{OH}(r)$
 respectively. The number of hydrogen bonds stays relatively constant
 across all of the models, but the coordination numbers of SSD/E and
 SSD/RF show an improvement over SSD1. This improvement is primarily
@@ -916,30 +919,31 @@ regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculat
 is the unit vector of the particle dipole.\cite{Rahman71} From these
 correlation functions, the orientational relaxation time of the dipole
 vector can be calculated from an exponential fit in the long-time
-regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculation of these
+regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these
 time constants were averaged from five detailed NVE simulations
 performed at the STP density for each of the respective models. It
 should be noted that the commonly cited value for $\tau_2$ of 1.9 ps
 was determined from the NMR data in reference \citen{Krynicki66} at a
-temperature near 34$^\circ$C.\cite{Rahman73} Because of the strong
+temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong
 temperature dependence of $\tau_2$, it is necessary to recalculate it
 at 298 K to make proper comparisons. The value shown in Table
 \ref{liquidproperties} was calculated from the same NMR data in the 
-fashion described in reference \citen{Krynicki66}. Again, SSD/E and
-SSD/RF show improved behavior over SSD1, both with and without an
-active reaction field. Turning on the reaction field leads to much
-improved time constants for SSD1; however, these results also include
-a corresponding decrease in system density. Numbers published from the
-original SSD dynamics studies appear closer to the experimental
-values, and this difference can be attributed to the use of the Ewald
-sum technique versus a reaction field.\cite{Ichiye99}
+fashion described in reference \citen{Krynicki66}. Similarly, $\tau_1$
+was recomputed for 298 K from the data in ref \citen{Eisenberg69}.
+Again, SSD/E and SSD/RF show improved behavior over SSD1, both with
+and without an active reaction field. Turning on the reaction field
+leads to much improved time constants for SSD1; however, these results
+also include a corresponding decrease in system density. Numbers
+published from the original SSD dynamics studies are shorter than the
+values observed here, and this difference can be attributed to the use
+of the Ewald sum technique versus a reaction field.\cite{Ichiye99}
 
 \subsection{Additional Observations}
 
 \begin{figure}
 \begin{center}
 \epsfxsize=6in
-\epsfbox{povIce.ps} 
+\epsfbox{icei_bw.eps} 
 \caption{A water lattice built from the crystal structure assumed by 
 SSD/E when undergoing an extremely restricted temperature NPT
 simulation. This form of ice is referred to as ice-{\it i} to
@@ -1019,7 +1023,7 @@ super-cooled liquid regimes.
 significantly lower than the densities obtained from other water
 models (and experiment). Analysis of self-diffusion showed SSD to
 capture the transport properties of water well in both the liquid and
-super-cooled liquid regimes.
+supercooled liquid regimes.
 
 In order to correct the density behavior, the original SSD model was
 reparameterized for use both with and without a reaction field (SSD/RF
@@ -1035,8 +1039,8 @@ effectively simulations of super-cooled or metastable 
 The existence of a novel low-density ice structure that is preferred
 by the SSD family of water models is somewhat troubling, since liquid
 simulations on this family of water models at room temperature are
-effectively simulations of super-cooled or metastable liquids.  One
-way to de-stabilize this unphysical ice structure would be to make the
+effectively simulations of supercooled or metastable liquids.  One
+way to destabilize this unphysical ice structure would be to make the
 range of angles preferred by the attractive part of the sticky
 potential much narrower.  This would require extensive
 reparameterization to maintain the same level of agreement with the