--- trunk/ssdePaper/nptSSD.tex	2004/02/05 22:54:57	1030
+++ trunk/ssdePaper/nptSSD.tex	2004/02/09 14:42:27	1040
@@ -1,4 +1,5 @@
 %\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
+%\documentclass[preprint,aps,endfloat]{revtex4}
 \documentclass[11pt]{article}
 \usepackage{endfloat}
 \usepackage{amsmath}
@@ -8,45 +9,45 @@
 \usepackage{tabularx}
 \usepackage{graphicx}
 \usepackage[ref]{overcite}
-%\usepackage{berkeley}
-%\usepackage{curves}
 \pagestyle{plain}
 \pagenumbering{arabic}
 \oddsidemargin 0.0cm \evensidemargin 0.0cm
 \topmargin -21pt \headsep 10pt
 \textheight 9.0in \textwidth 6.5in
 \brokenpenalty=10000
-\renewcommand{\baselinestretch}{1.2}
-\renewcommand\citemid{\ } % no comma in optional reference note
 
+%\renewcommand\citemid{\ } % no comma in optional reference note
+
 \begin{document}
 
 \title{On the structural and transport properties of the soft sticky
 dipole (SSD) and related single point water models}
 
-\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
-Department of Chemistry and Biochemistry\\ University of Notre Dame\\
+\author{Christopher J. Fennell and J. Daniel
+Gezelter\footnote{Corresponding author. \ Electronic mail:
+gezelter@nd.edu} \\ Department of Chemistry and Biochemistry\\ University of Notre Dame\\
 Notre Dame, Indiana 46556}
 
 \date{\today}
 
 \maketitle
+\doublespacing
 
 \begin{abstract}
 The density maximum and temperature dependence of the self-diffusion
 constant were investigated for the soft sticky dipole (SSD) water
-model and two related re-parameterizations of this single-point model.
+model and two related reparameterizations of this single-point model.
 A combination of microcanonical and isobaric-isothermal molecular
 dynamics simulations were used to calculate these properties, both
 with and without the use of reaction field to handle long-range
 electrostatics.  The isobaric-isothermal (NPT) simulations of the
 melting of both ice-$I_h$ and ice-$I_c$ showed a density maximum near
-260 K.  In most cases, the use of the reaction field resulted in
+260~K.  In most cases, the use of the reaction field resulted in
 calculated densities which were were significantly lower than
 experimental densities.  Analysis of self-diffusion constants shows
 that the original SSD model captures the transport properties of
-experimental water very well in both the normal and super-cooled
-liquid regimes.  We also present our re-parameterized versions of SSD
+experimental water very well in both the normal and supercooled
+liquid regimes.  We also present our reparameterized versions of SSD
 for use both with the reaction field or without any long-range
 electrostatic corrections.  These are called the SSD/RF and SSD/E
 models respectively.  These modified models were shown to maintain or
@@ -62,7 +63,6 @@ family.
 
 %\narrowtext 
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %	                       BODY OF TEXT
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -94,7 +94,7 @@ which has an interaction site that is both a point dip
 was developed by Ichiye \emph{et al.} as a modified form of the
 hard-sphere water model proposed by Bratko, Blum, and
 Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model
-which has an interaction site that is both a point dipole along with a
+which has an interaction site that is both a point dipole and a
 Lennard-Jones core.  However, since the normal aligned and
 anti-aligned geometries favored by point dipoles are poor mimics of
 local structure in liquid water, a short ranged ``sticky'' potential
@@ -164,12 +164,12 @@ simulations using this model, Ichiye {\it et al.} repo
 Since SSD is a single-point {\it dipolar} model, the force
 calculations are simplified significantly relative to the standard
 {\it charged} multi-point models. In the original Monte Carlo
-simulations using this model, Ichiye {\it et al.} reported that using
-SSD decreased computer time by a factor of 6-7 compared to other
+simulations using this model, Liu and Ichiye reported that using SSD
+decreased computer time by a factor of 6-7 compared to other
 models.\cite{Ichiye96} What is most impressive is that this savings
 did not come at the expense of accurate depiction of the liquid state
-properties.  Indeed, SSD maintains reasonable agreement with the
-Soper data for the structural features of liquid
+properties.  Indeed, SSD maintains reasonable agreement with the Soper
+data for the structural features of liquid
 water.\cite{Soper86,Ichiye96} Additionally, the dynamical properties
 exhibited by SSD agree with experiment better than those of more
 computationally expensive models (like TIP3P and
@@ -203,16 +203,16 @@ utilizing the Reaction Field.
 follows, we compare our reparamaterization of SSD with both the
 original SSD and SSD1 models with the goal of improving
 the bulk phase behavior of an SSD-derived model in simulations
-utilizing the Reaction Field.
+utilizing the reaction field.
 
 \section{Methods}
 
 Long-range dipole-dipole interactions were accounted for in this study
-by using either the reaction field method or by resorting to a simple
-cubic switching function at a cutoff radius.  The reaction field
-method was actually first used in Monte Carlo simulations of liquid
-water.\cite{Barker73} Under this method, the magnitude of the reaction
-field acting on dipole $i$ is
+by using either the reaction field technique or by resorting to a
+simple cubic switching function at a cutoff radius.  One of the early
+applications of a reaction field was actually in Monte Carlo
+simulations of liquid water.\cite{Barker73} Under this method, the
+magnitude of the reaction 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} s(r_{ij}),
@@ -226,12 +226,12 @@ field is known to alter the bulk orientational propert
 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
 \mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley}  Use of the reaction
-field is known to alter the bulk orientational properties, such as the
-dielectric relaxation time.  There is particular sensitivity of this
-property on changes in the length of the cutoff
-radius.\cite{Berendsen98} This variable behavior makes reaction field
-a less attractive method than the Ewald sum.  However, for very large
-systems, the computational benefit of reaction field is dramatic.
+field is known to alter the bulk orientational properties of simulated
+water, and there is particular sensitivity of these properties on
+changes in the length of the cutoff radius.\cite{Berendsen98} This
+variable behavior makes reaction field a less attractive method than
+the Ewald sum.  However, for very large systems, the computational
+benefit of reaction field is dramatic.
  
 We have also performed a companion set of simulations {\it without} a
 surrounding dielectric (i.e. using a simple cubic switching function
@@ -287,10 +287,10 @@ time steps is illustrated in figure
 \begin{center}
 \epsfxsize=6in
 \epsfbox{timeStep.epsi}
-\caption{Energy conservation using both quaternion-based integration and
-the {\sc dlm} method 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.}
+\caption{Energy conservation using both quaternion-based integration and the
+{\sc dlm} method 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}
@@ -299,38 +299,38 @@ handle orientational motion. At time steps of 0.1 and 
 steps for both the {\sc dlm} 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
+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 {\sc dlm} method are
+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 {\sc dlm} 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 simulations using {\sc dlm} 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
+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
+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$
-crystals were formed by first arranging the centers of mass of the
-SSD particles into a ``hexagonal'' ice lattice of 1024
-particles. Because of the crystal structure of $I_h$ ice, the
-simulation box assumed an orthorhombic shape with an edge length ratio
-of approximately 1.00$\times$1.06$\times$1.23. The particles were then
-allowed to orient freely about fixed positions with angular momenta
-randomized at 400 K for varying times. The rotational temperature was
-then scaled down in stages to slowly cool the crystals to 25 K. The
-particles were then allowed to translate with fixed orientations at a
-constant pressure of 1 atm for 50 ps at 25 K. Finally, all constraints
-were removed and the ice crystals were allowed to equilibrate for 50
-ps at 25 K and a constant pressure of 1 atm.  This procedure resulted
-in structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
+crystals were formed by first arranging the centers of mass of the SSD
+particles into a ``hexagonal'' ice lattice of 1024 particles. Because
+of the crystal structure of $I_h$ ice, the simulation box assumed an
+orthorhombic shape with an edge length ratio of approximately
+1.00$\times$1.06$\times$1.23. The particles were then allowed to
+orient freely about fixed positions with angular momenta randomized at
+400~K for varying times. The rotational temperature was then scaled
+down in stages to slowly cool the crystals to 25~K. The particles were
+then allowed to translate with fixed orientations at a constant
+pressure of 1 atm for 50~ps at 25~K. Finally, all constraints were
+removed and the ice crystals were allowed to equilibrate for 50~ps at
+25~K and a constant pressure of 1~atm.  This procedure resulted in
+structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
 rules.\cite{Bernal33,Rahman72} This method was also utilized in the
 making of diamond lattice $I_c$ ice crystals, with each cubic
 simulation box consisting of either 512 or 1000 particles. Only
@@ -347,32 +347,32 @@ for 100 ps prior to a 200 ps data collection run at ea
 supercooled regime. An ensemble average from five separate melting
 simulations was acquired, each starting from different ice crystals
 generated as described previously. All simulations were equilibrated
-for 100 ps prior to a 200 ps data collection run at each temperature
-setting. The temperature range of study spanned from 25 to 400 K, with
-a maximum degree increment of 25 K. For regions of interest along this
-stepwise progression, the temperature increment was decreased from 25
-K to 10 and 5 K.  The above equilibration and production times were
+for 100~ps prior to a 200~ps data collection run at each temperature
+setting. The temperature range of study spanned from 25 to 400~K, with
+a maximum degree increment of 25~K. For regions of interest along this
+stepwise progression, the temperature increment was decreased from
+25~K to 10 and 5~K.  The above equilibration and production times were
 sufficient in that fluctuations in the volume autocorrelation function
-were damped out in all simulations in under 20 ps.
+were damped out in all simulations in under 20~ps.
 
 \subsection{Density Behavior}
 
 Our initial simulations focused on the original SSD water model,
 and an average density versus temperature plot is shown in figure
 \ref{dense1}. Note that the density maximum when using a reaction
-field appears between 255 and 265 K.  There were smaller fluctuations
-in the density at 260 K than at either 255 or 265, so we report this
+field appears between 255 and 265~K.  There were smaller fluctuations
+in the density at 260~K than at either 255 or 265~K, so we report this
 value as the location of the density maximum. Figure \ref{dense1} was
 constructed using ice $I_h$ crystals for the initial configuration;
 though not pictured, the simulations starting from ice $I_c$ crystal
 configurations showed similar results, with a liquid-phase density
-maximum in this same region (between 255 and 260 K). 
+maximum in this same region (between 255 and 260~K). 
 
 \begin{figure}
 \begin{center}
 \epsfxsize=6in
 \epsfbox{denseSSDnew.eps}
-\caption{Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}],
+\caption{ Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}],
  TIP3P [Ref. \citen{Jorgensen98b}], SPC/E [Ref. \citen{Clancy94}], SSD
  without Reaction Field, SSD, and experiment [Ref. \citen{CRC80}]. The
  arrows indicate the change in densities observed when turning off the
@@ -398,8 +398,8 @@ cutoff radius of 12.0 \AA\ to complement the previous 
 dependent on the cutoff radius used both with and without the use of
 reaction field.\cite{Berendsen98} In order to address the possible
 effect of cutoff radius, simulations were performed with a dipolar
-cutoff radius of 12.0 \AA\ to complement the previous SSD
-simulations, all performed with a cutoff of 9.0 \AA. All of the
+cutoff radius of 12.0~\AA\ to complement the previous SSD
+simulations, all performed with a cutoff of 9.0~\AA. All of the
 resulting densities overlapped within error and showed no significant
 trend toward lower or higher densities as a function of cutoff radius,
 for simulations both with and without reaction field. These results
@@ -412,7 +412,7 @@ apparent at temperatures greater than 300 K. Lower tha
 scaling of SSD relative to other common models at any given
 temperature. SSD assumes a lower density than any of the other
 listed models at the same pressure, behavior which is especially
-apparent at temperatures greater than 300 K. Lower than expected
+apparent at temperatures greater than 300~K. Lower than expected
 densities have been observed for other systems using a reaction field
 for long-range electrostatic interactions, so the most likely reason
 for the significantly lower densities seen in these simulations is the
@@ -425,11 +425,11 @@ however, a shift in the density maximum location, down
 freezing point of liquid water. The shape of the curve is similar to
 the curve produced from SSD simulations using reaction field,
 specifically the rapidly decreasing densities at higher temperatures;
-however, a shift in the density maximum location, down to 245 K, is
+however, a shift in the density maximum location, down to 245~K, is
 observed. This is a more accurate comparison to the other listed water
 models, in that no long range corrections were applied in those
 simulations.\cite{Clancy94,Jorgensen98b} However, even without the
-reaction field, the density around 300 K is still significantly lower
+reaction field, the density around 300~K is still significantly lower
 than experiment and comparable water models. This anomalous behavior
 was what lead Tan {\it et al.} to recently reparameterize
 SSD.\cite{Ichiye03} Throughout the remainder of the paper our
@@ -444,8 +444,8 @@ underwent 50 ps of temperature scaling and 50 ps of co
 constant energy (NVE) simulations were performed at the average
 density obtained by the NPT simulations at an identical target
 temperature. Simulations started with randomized velocities and
-underwent 50 ps of temperature scaling and 50 ps of constant energy
-equilibration before a 200 ps data collection run. Diffusion constants
+underwent 50~ps of temperature scaling and 50~ps of constant energy
+equilibration before a 200~ps data collection run. Diffusion constants
 were calculated via linear fits to the long-time behavior of the
 mean-square displacement as a function of time. The averaged results
 from five sets of NVE simulations are displayed in figure
@@ -456,7 +456,7 @@ results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} 
 \begin{center}
 \epsfxsize=6in
 \epsfbox{betterDiffuse.epsi} 
-\caption{Average self-diffusion constant as a function of temperature for 
+\caption{ Average self-diffusion constant as a function of temperature for 
 SSD, SPC/E [Ref. \citen{Clancy94}], and TIP5P
 [Ref. \citen{Jorgensen01}] compared with experimental data
 [Refs. \citen{Gillen72} and \citen{Holz00}]. Of the three water models
@@ -471,7 +471,7 @@ reproducing values similar to experiment around 290 K;
 strengths of the SSD model. Of the three models shown, the SSD model
 has the most accurate depiction of self-diffusion in both the
 supercooled and liquid regimes.  SPC/E does a respectable job by
-reproducing values similar to experiment around 290 K; however, it
+reproducing values similar to experiment around 290~K; however, it
 deviates at both higher and lower temperatures, failing to predict the
 correct thermal trend. TIP5P and SSD both start off low at colder
 temperatures and tend to diffuse too rapidly at higher temperatures.
@@ -490,26 +490,17 @@ at 245 K, indicating a first order phase transition fo
 capacity (C$_\text{p}$) was monitored to locate the melting transition
 in each of the simulations. In the melting simulations of the 1024
 particle ice $I_h$ simulations, a large spike in C$_\text{p}$ occurs
-at 245 K, indicating a first order phase transition for the melting of
+at 245~K, indicating a first order phase transition for the melting of
 these ice crystals. When the reaction field is turned off, the melting
-transition occurs at 235 K.  These melting transitions are
+transition occurs at 235~K.  These melting transitions are
 considerably lower than the experimental value.
-
-\begin{figure}
-\begin{center}
-\epsfxsize=6in
-\epsfbox{corrDiag.eps} 
-\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
-\label{corrAngle}
-\end{center}
-\end{figure}
 
 \begin{figure}
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{fullContours.eps} 
-\caption{Contour plots of 2D angular pair correlation functions for
-512 SSD molecules at 100 K (A \& B) and 300 K (C \& D). Dark areas
+\caption{ Contour plots of 2D angular pair correlation functions for
+512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
 signify regions of enhanced density while light areas signify
 depletion relative to the bulk density. White areas have pair
 correlation values below 0.5 and black areas have values above 1.5.}
@@ -517,6 +508,15 @@ Additional analysis of the melting process was perform
 \end{center}
 \end{figure} 
 
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{corrDiag.eps} 
+\caption{ An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
+\label{corrAngle}
+\end{center}
+\end{figure}
+
 Additional analysis of the melting process was performed using
 two-dimensional structure and dipole angle correlations. Expressions
 for these correlations are as follows:
@@ -555,25 +555,25 @@ this split character alters to show the leading 4 \AA\
 $g_\mathrm{OO}(r)$.\cite{Ichiye96} At low temperatures, the second
 solvation shell peak appears to have two distinct components that
 blend together to form one observable peak. At higher temperatures,
-this split character alters to show the leading 4 \AA\ peak dominated
+this split character alters to show the leading 4~\AA\ peak dominated
 by equatorial anti-parallel dipole orientations. There is also a
 tightly bunched group of axially arranged dipoles that most likely
 consist of the smaller fraction of aligned dipole pairs. The trailing
-component of the split peak at 5 \AA\ is dominated by aligned dipoles
+component of the split peak at 5~\AA\ is dominated by aligned dipoles
 that assume hydrogen bond arrangements similar to those seen in the
 first solvation shell. This evidence indicates that the dipole pair
 interaction begins to dominate outside of the range of the dipolar
 repulsion term.  The energetically favorable dipole arrangements
 populate the region immediately outside this repulsion region (around
-4 \AA), while arrangements that seek to satisfy both the sticky and
+4~\AA), while arrangements that seek to satisfy both the sticky and
 dipole forces locate themselves just beyond this initial buildup
-(around 5 \AA).
+(around 5~\AA).
 
 From these findings, the split second peak is primarily the product of
 the dipolar repulsion term of the sticky potential. In fact, the inner
 peak can be pushed out and merged with the outer split peak just by
 extending the switching function ($s^\prime(r_{ij})$) from its normal
-4.0 \AA\ cutoff to values of 4.5 or even 5 \AA. This type of
+4.0~\AA\ cutoff to values of 4.5 or even 5~\AA. This type of
 correction is not recommended for improving the liquid structure,
 since the second solvation shell would still be shifted too far
 out. In addition, this would have an even more detrimental effect on
@@ -608,11 +608,11 @@ the liquid structure in simulations without a long-ran
 
 \begin{table}
 \begin{center}
-\caption{Parameters for the original and adjusted models}
+\caption{ Parameters for the original and adjusted models}
 \begin{tabular}{ l  c  c  c  c } 
 \hline \\[-3mm]
 \ \ \ Parameters\ \ \  & \ \ \ SSD [Ref. \citen{Ichiye96}] \ \ \
-& \ SSD1 [Ref. \citen{Ichiye03}]\ \  & \ SSD/E\ \  & \ SSD/RF \\
+& \ SSD1 [Ref. \citen{Ichiye03}]\ \  & \ SSD/E\ \  & \ \ SSD/RF \\
 \hline \\[-3mm]
 \ \ \ $\sigma$ (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
 \ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
@@ -632,7 +632,7 @@ the liquid structure in simulations without a long-ran
 \begin{center} 
 \epsfxsize=5in
 \epsfbox{GofRCompare.epsi} 
-\caption{Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with 
+\caption{ Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with 
 SSD/E and SSD1 without reaction field (top), as well as
 SSD/RF and SSD1 with reaction field turned on
 (bottom). The insets show the respective first peaks in detail. Note
@@ -646,7 +646,7 @@ peak of SSD/E and SSD/RF.}
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{dualsticky_bw.eps} 
-\caption{Positive and negative isosurfaces of the sticky potential for
+\caption{ Positive and negative 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 component.}
@@ -688,27 +688,26 @@ persistence of full dipolar character below the previo
 particles feel the pull of the ``hydrogen bonds''. Aside from
 improving the shape of the first peak in the g(\emph{r}), this
 modification improves the densities considerably by allowing the
-persistence of full dipolar character below the previous 4.0 \AA\
+persistence of full dipolar character below the previous 4.0~\AA\
 cutoff.
 
 While adjusting the location and shape of the first peak of $g(r)$
 improves the densities, these changes alone are insufficient to bring
 the system densities up to the values observed experimentally.  To
 further increase the densities, the dipole moments were increased in
-both of our adjusted models. Since SSD is a dipole based model,
-the structure and transport are very sensitive to changes in the
-dipole moment. The original SSD simply used the dipole moment
-calculated from the TIP3P water model, which at 2.35 D is
-significantly greater than the experimental gas phase value of 1.84
-D. The larger dipole moment is a more realistic value and improves the
-dielectric properties of the fluid. Both theoretical and experimental
-measurements indicate a liquid phase dipole moment ranging from 2.4 D
-to values as high as 3.11 D, providing a substantial range of
-reasonable values for a dipole
-moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
-increasing the dipole moments to 2.42 and 2.48 D for SSD/E and
-SSD/RF, respectively, leads to significant changes in the
-density and transport of the water models.
+both of our adjusted models. Since SSD is a dipole based model, the
+structure and transport are very sensitive to changes in the dipole
+moment. The original SSD simply used the dipole moment calculated from
+the TIP3P water model, which at 2.35~D is significantly greater than
+the experimental gas phase value of 1.84~D. The larger dipole moment
+is a more realistic value and improves the dielectric properties of
+the fluid. Both theoretical and experimental measurements indicate a
+liquid phase dipole moment ranging from 2.4~D to values as high as
+3.11~D, providing a substantial range of reasonable values for a
+dipole moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
+increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF,
+respectively, leads to significant changes in the density and
+transport of the water models.
 
 In order to demonstrate the benefits of these reparameterizations, a
 series of NPT and NVE simulations were performed to probe the density
@@ -719,7 +718,7 @@ simulation was equilibrated for 100 ps before a 200 ps
 results are obtained from five separate simulations of 1024 particle
 systems, and the melting sequences were started from different ice
 $I_h$ crystals constructed as described previously. Each NPT
-simulation was equilibrated for 100 ps before a 200 ps data collection
+simulation was equilibrated for 100~ps before a 200~ps data collection
 run at each temperature step, and the final configuration from the
 previous temperature simulation was used as a starting point. All NVE
 simulations had the same thermalization, equilibration, and data
@@ -729,7 +728,7 @@ collection times as stated previously.
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdeDense.epsi} 
-\caption{Comparison of densities calculated with SSD/E to 
+\caption{ Comparison of densities calculated with SSD/E to 
 SSD1 without a reaction field, TIP3P [Ref. \citen{Jorgensen98b}],
 TIP5P [Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}] and
 experiment [Ref. \citen{CRC80}]. The window shows a expansion around
@@ -740,7 +739,7 @@ Fig. \ref{ssdedense} shows the density profile for the
 \end{center}
 \end{figure} 
 
-Fig. \ref{ssdedense} shows the density profile for the SSD/E
+Figure \ref{ssdedense} shows the density profile for the SSD/E
 model in comparison to SSD1 without a reaction field, other
 common water models, and experimental results. The calculated
 densities for both SSD/E and SSD1 have increased
@@ -753,20 +752,20 @@ improved the structuring of the liquid (as seen in fig
 better than the SSD value of 0.967$\pm$0.003 g/cm$^3$. The
 changes to the dipole moment and sticky switching functions have
 improved the structuring of the liquid (as seen in figure
-\ref{grcompare}, but they have shifted the density maximum to much
+\ref{grcompare}), but they have shifted the density maximum to much
 lower temperatures. This comes about via an increase in the liquid
 disorder through the weakening of the sticky potential and
 strengthening of the dipolar character. However, this increasing
 disorder in the SSD/E model has little effect on the melting
 transition. By monitoring $C_p$ throughout these simulations, the
-melting transition for SSD/E was shown to occur at 235 K.  The
+melting transition for SSD/E was shown to occur at 235~K.  The
 same transition temperature observed with SSD and SSD1.
 
 \begin{figure}
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdrfDense.epsi} 
-\caption{Comparison of densities calculated with SSD/RF to 
+\caption{ Comparison of densities calculated with SSD/RF to 
 SSD1 with a reaction field, TIP3P [Ref. \citen{Jorgensen98b}],
 TIP5P [Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}], and
 experiment [Ref. \citen{CRC80}]. The inset shows the necessity of
@@ -790,17 +789,17 @@ which observed at 245 K for SSD/RF, is identical to SS
 further emphasize the importance of reparameterization in order to
 model the density properly under different simulation conditions.
 Again, these changes have only a minor effect on the melting point,
-which observed at 245 K for SSD/RF, is identical to SSD and only 5 K
+which observed at 245~K for SSD/RF, is identical to SSD and only 5~K
 lower than SSD1 with a reaction field. Additionally, the difference in
 density maxima is not as extreme, with SSD/RF showing a density
-maximum at 255 K, fairly close to the density maxima of 260 K and 265
-K, shown by SSD and SSD1 respectively.
+maximum at 255~K, fairly close to the density maxima of 260~K and
+265~K, shown by SSD and SSD1 respectively.
 
 \begin{figure}
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdeDiffuse.epsi} 
-\caption{The diffusion constants calculated from SSD/E and 
+\caption{ The diffusion constants calculated from SSD/E and 
 SSD1 (both without a reaction field) along with experimental results
 [Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations were
 performed at the average densities observed in the 1 atm NPT
@@ -823,8 +822,8 @@ K). Both models follow the shape of the experimental c
 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
+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
@@ -841,7 +840,7 @@ conditions.
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdrfDiffuse.epsi} 
-\caption{The diffusion constants calculated from SSD/RF and 
+\caption{ The diffusion constants calculated from SSD/RF and 
 SSD1 (both with an active reaction field) along with
 experimental results [Refs. \citen{Gillen72} and \citen{Holz00}]. The
 NVE calculations were performed at the average densities observed in
@@ -860,7 +859,7 @@ temperatures greater than 330 K.  As stated above, thi
 throughout most of the temperature range shown and exhibiting only a
 slight increasing trend at higher temperatures. SSD1 tends to diffuse
 more slowly at low temperatures and deviates to diffuse too rapidly at
-temperatures greater than 330 K.  As stated above, this deviation away
+temperatures greater than 330~K.  As stated above, this deviation away
 from the ideal trend is due to a rapid decrease in density at higher
 temperatures. SSD/RF does not suffer from this problem as much as SSD1
 because the calculated densities are closer to the experimental
@@ -871,27 +870,25 @@ reparameterization when using an altered long-range co
 \begin{minipage}{\linewidth}
 \renewcommand{\thefootnote}{\thempfootnote}
 \begin{center}
-\caption{Properties of the single-point water models compared with
-experimental data at ambient conditions} 
+\caption{ Properties of the single-point water models compared with
+experimental data at ambient conditions. Deviations of the of the
+averages are given in parentheses.}
 \begin{tabular}{ l  c  c  c  c  c } 
 \hline \\[-3mm]
-\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ 
-\ & \ SSD/RF \ \ \ & \ Expt. \\
+\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ \ SSD/E \ \ \ & \ \ SSD1 (RF) \ \ 
+\ & \ \ SSD/RF \ \ \ & \ \ Expt. \\
  \hline \\[-3mm]
-\ \ \ $\rho$ (g/cm$^3$) & 0.999 $\pm$0.001 & 0.996 $\pm$0.001 & 0.972 $\pm$0.002 & 0.997 $\pm$0.001 & 0.997 \\
-\ \ \ $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\cite{Mills73} \\
-\ \ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 &
+\ \ $\rho$ (g/cm$^3$) & 0.999(0.001) & 0.996(0.001) & 0.972(0.002) & 0.997(0.001) & 0.997 \\
+\ \ $C_p$ (cal/mol K) & 28.80(0.11) & 25.45(0.09) & 28.28(0.06) & 23.83(0.16) & 17.98 \\
+\ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78(0.7) & 2.51(0.18) & 2.00(0.17) & 2.32(0.06) & 2.299\cite{Mills73} \\
+\ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 &
 4.7\footnote{Calculated by integrating $g_{\text{OO}}(r)$ in
 Ref. \citen{Head-Gordon00_1}} \\ 
-\ \ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 &
+\ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 &
 3.5\footnote{Calculated by integrating $g_{\text{OH}}(r)$ in
 Ref. \citen{Soper86}}  \\
-\ \ \ $\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\footnote{Calculated for 298 K from data in Ref. \citen{Eisenberg69}} \\
-\ \ \ $\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\footnote{Calculated for 298 K from data in
+\ \ $\tau_1$ (ps) & 10.9(0.6) & 7.3(0.4) & 7.5(0.7) & 7.2(0.4) & 5.7\footnote{Calculated for 298 K from data in Ref. \citen{Eisenberg69}} \\
+\ \ $\tau_2$ (ps) & 4.7(0.4) & 3.1(0.2) & 3.5(0.3) & 3.2(0.2) & 2.3\footnote{Calculated for 298 K from data in
 Ref. \citen{Krynicki66}}
 \end{tabular}
 \label{liquidproperties}
@@ -943,11 +940,11 @@ dependence of $\tau_2$, it is necessary to recalculate
 the commonly cited value of 1.9 ps for $\tau_2$ was determined from
 the NMR data in Ref. \citen{Krynicki66} at a 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
+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 Ref. \citen{Krynicki66}. Similarly, $\tau_1$ was
-recomputed for 298 K from the data in Ref. \citen{Eisenberg69}.
+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
@@ -962,7 +959,7 @@ can be attributed to the use of the Ewald sum.\cite{Ic
 \begin{center}
 \epsfxsize=6in
 \epsfbox{icei_bw.eps} 
-\caption{The most stable crystal structure assumed by the SSD family
+\caption{ The most stable crystal structure assumed by the SSD family
 of water models.  We refer to this structure as Ice-{\it i} to
 indicate its origins in computer simulation.  This image was taken of
 the (001) face of the crystal.}
@@ -973,9 +970,9 @@ water.  After melting at 235 K, two of five systems un
 While performing a series of melting simulations on an early iteration
 of SSD/E not discussed in this paper, we observed
 recrystallization into a novel structure not previously known for
-water.  After melting at 235 K, two of five systems underwent
-crystallization events near 245 K.  The two systems remained
-crystalline up to 320 and 330 K, respectively.  The crystal exhibits
+water.  After melting at 235~K, two of five systems underwent
+crystallization events near 245~K.  The two systems remained
+crystalline up to 320 and 330~K, respectively.  The crystal exhibits
 an expanded zeolite-like structure that does not correspond to any
 known form of ice.  This appears to be an artifact of the point
 dipolar models, so to distinguish it from the experimentally observed
@@ -1009,12 +1006,12 @@ family, Ice-{\it i} had the lowest calculated enthalpy
 
 \begin{table}
 \begin{center}
-\caption{Enthalpies of Formation (in kcal / mol) of the three crystal
+\caption{ Enthalpies of Formation (in kcal / mol) of the three crystal
 structures (at 1 K) exhibited by the SSD family of water models} 
 \begin{tabular}{ l  c  c  c  } 
 \hline \\[-3mm]
-\ \ \ Water Model \ \ \  & \ \ \ Ice-$I_h$ \ \ \ & \ Ice-$I_c$\ \  & \
-Ice-{\it i} \\
+\ \ \ Water Model \ \ \  & \ \ \ Ice-$I_h$ \ \ \ & \ \ \ Ice-$I_c$ \ \ \  & 
+\ \ \ \ Ice-{\it i} \\
 \hline \\[-3mm]
 \ \ \ SSD/E & -72.444 & -72.450 & -73.748 \\
 \ \ \ SSD/RF & -73.093 & -73.075 & -74.180 \\
@@ -1026,7 +1023,7 @@ performed with ice-{\it i} as the initial configuratio
 \end{table}
 
 In addition to these energetic comparisons, melting simulations were
-performed with ice-{\it i} as the initial configuration using SSD/E,
+performed with Ice-{\it i} as the initial configuration using SSD/E,
 SSD/RF, and SSD1 both with and without a reaction field. The melting
 transitions for both SSD/E and SSD1 without reaction field occurred at
 temperature in excess of 375~K.  SSD/RF and SSD1 with a reaction field
@@ -1084,8 +1081,7 @@ DMR-0079647.
 \newpage
 
 \bibliographystyle{jcp}
-\bibliography{nptSSD} 
+\bibliography{nptSSD}
 
-%\pagebreak 
 
 \end{document}
