--- trunk/ssdePaper/nptSSD.tex	2004/02/05 20:47:50	1029
+++ trunk/ssdePaper/nptSSD.tex	2004/02/06 19:15:44	1033
@@ -1,63 +1,67 @@
 %\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
-\documentclass[11pt]{article}
-\usepackage{endfloat}
+\documentclass[preprint,aps]{revtex4}
+%\documentclass[11pt]{article}
+%\usepackage{endfloat}
 \usepackage{amsmath}
 \usepackage{epsf}
 \usepackage{berkeley}
 \usepackage{setspace}
 \usepackage{tabularx}
 \usepackage{graphicx}
-\usepackage[ref]{overcite}
-%\usepackage{berkeley}
+%\usepackage[ref]{overcite}
 %\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
+%\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
+\newcounter{captions}
+\newcounter{figs}
 
 \begin{document}
 
 \title{On the structural and transport properties of the soft sticky
-dipole ({\sc ssd}) and related single point water models}
+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}}
+
+\affiliation{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
 Notre Dame, Indiana 46556}
 
 \date{\today}
 
-\maketitle
 
 \begin{abstract}
 The density maximum and temperature dependence of the self-diffusion
-constant were investigated for the soft sticky dipole ({\sc ssd}) water
-model and two related re-parameterizations of this single-point model.
+constant were investigated for the soft sticky dipole (SSD) water
+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 {\sc ssd} model captures the transport properties of
+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 {\sc ssd}
+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 {\sc ssd/rf} and {\sc ssd/e}
+electrostatic corrections.  These are called the SSD/RF and SSD/E
 models respectively.  These modified models were shown to maintain or
 improve upon the experimental agreement with the structural and
-transport properties that can be obtained with either the original {\sc ssd}
-or the density corrected version of the original model ({\sc ssd1}).
+transport properties that can be obtained with either the original SSD
+or the density corrected version of the original model (SSD1).
 Additionally, a novel low-density ice structure is presented
-which appears to be the most stable ice structure for the entire {\sc ssd}
+which appears to be the most stable ice structure for the entire SSD
 family.
 \end{abstract}
 
+\maketitle
+
 \newpage
 
 %\narrowtext 
@@ -89,20 +93,19 @@ cost is the Soft Sticky Dipole ({\sc ssd}) water
 
 One recently developed model that largely succeeds in retaining the
 accuracy of bulk properties while greatly reducing the computational
-cost is the Soft Sticky Dipole ({\sc ssd}) water
-model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The {\sc ssd} model was
-developed by Ichiye \emph{et al.} as a modified form of the
+cost is the Soft Sticky Dipole (SSD) water
+model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The SSD model
+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} {\sc ssd} is a {\it single point} model which
-has an interaction site that is both a point dipole along with a
+Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model
+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
 is also added.  The sticky potential directs the molecules to assume
-the proper hydrogen bond orientation in the first solvation
-shell.  
+the proper hydrogen bond orientation in the first solvation shell.
 
-The interaction between two {\sc ssd} water molecules \emph{i} and \emph{j}
+The interaction between two SSD water molecules \emph{i} and \emph{j}
 is given by the potential
 \begin{equation}
 u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
@@ -159,60 +162,61 @@ can be found in the original {\sc ssd}
 enhances the tetrahedral geometry for hydrogen bonded structures),
 while $w^\prime$ is a purely empirical function.  A more detailed
 description of the functional parts and variables in this potential
-can be found in the original {\sc ssd}
+can be found in the original SSD
 articles.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03}
 
-Since {\sc ssd} is a single-point {\it dipolar} model, the force
+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
-{\sc 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, {\sc ssd} maintains reasonable agreement with the Soper
+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 {\sc ssd} agree with experiment better than those of more
+exhibited by SSD agree with experiment better than those of more
 computationally expensive models (like TIP3P and
 SPC/E).\cite{Ichiye99} The combination of speed and accurate depiction
-of solvent properties makes {\sc ssd} a very attractive model for the
+of solvent properties makes SSD a very attractive model for the
 simulation of large scale biochemical simulations.
 
-One feature of the {\sc ssd} model is that it was parameterized for use with
-the Ewald sum to handle long-range interactions.  This would normally
-be the best way of handling long-range interactions in systems that
-contain other point charges.  However, our group has recently become
-interested in systems with point dipoles as mimics for neutral, but
-polarized regions on molecules (e.g. the zwitterionic head group
-regions of phospholipids).  If the system of interest does not contain
-point charges, the Ewald sum and even particle-mesh Ewald become
-computational bottlenecks.  Their respective ideal $N^\frac{3}{2}$ and
-$N\log N$ calculation scaling orders for $N$ particles can become
-prohibitive when $N$ becomes large.\cite{Darden99} In applying this
-water model in these types of systems, it would be useful to know its
-properties and behavior under the more computationally efficient
-reaction field (RF) technique, or even with a simple cutoff. This
-study addresses these issues by looking at the structural and
-transport behavior of {\sc ssd} over a variety of temperatures with the
-purpose of utilizing the RF correction technique.  We then suggest
-modifications to the parameters that result in more realistic bulk
-phase behavior.  It should be noted that in a recent publication, some
-of the original investigators of the {\sc ssd} water model have suggested
-adjustments to the {\sc ssd} water model to address abnormal density
-behavior (also observed here), calling the corrected model
-{\sc ssd1}.\cite{Ichiye03} In what follows, we compare our
-reparamaterization of {\sc ssd} with both the original {\sc ssd} and {\sc ssd1} models
-with the goal of improving the bulk phase behavior of an {\sc ssd}-derived
-model in simulations utilizing the Reaction Field.
+One feature of the SSD model is that it was parameterized for
+use with the Ewald sum to handle long-range interactions.  This would
+normally be the best way of handling long-range interactions in
+systems that contain other point charges.  However, our group has
+recently become interested in systems with point dipoles as mimics for
+neutral, but polarized regions on molecules (e.g. the zwitterionic
+head group regions of phospholipids).  If the system of interest does
+not contain point charges, the Ewald sum and even particle-mesh Ewald
+become computational bottlenecks.  Their respective ideal
+$N^\frac{3}{2}$ and $N\log N$ calculation scaling orders for $N$
+particles can become prohibitive when $N$ becomes
+large.\cite{Darden99} In applying this water model in these types of
+systems, it would be useful to know its properties and behavior under
+the more computationally efficient reaction field (RF) technique, or
+even with a simple cutoff. This study addresses these issues by
+looking at the structural and transport behavior of SSD over a
+variety of temperatures with the purpose of utilizing the RF
+correction technique.  We then suggest modifications to the parameters
+that result in more realistic bulk phase behavior.  It should be noted
+that in a recent publication, some of the original investigators of
+the SSD water model have suggested adjustments to the SSD
+water model to address abnormal density behavior (also observed here),
+calling the corrected model SSD1.\cite{Ichiye03} In what
+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.
 
 \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,18 +230,18 @@ 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
 at the cutoff radius), and as a result we have two reparamaterizations
-of {\sc ssd} which could be used either with or without the reaction field
-turned on.
+of SSD which could be used either with or without the reaction
+field turned on.
 
 Simulations to obtain the preferred densities of the models were
 performed in the isobaric-isothermal (NPT) ensemble, while all
@@ -247,18 +251,19 @@ simulations of {\sc ssd}, because there are no explici
 implemented using an integral thermostat and barostat as outlined by
 Hoover.\cite{Hoover85,Hoover86} All molecules were treated as
 non-linear rigid bodies. Vibrational constraints are not necessary in
-simulations of {\sc ssd}, because there are no explicit hydrogen atoms, and
-thus no molecular vibrational modes need to be considered.
+simulations of SSD, because there are no explicit hydrogen
+atoms, and thus no molecular vibrational modes need to be considered.
 
 Integration of the equations of motion was carried out using the
 symplectic splitting method proposed by Dullweber, Leimkuhler, and
-McLachlan ({\sc 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 {\sc dlm}
-method (fig. \ref{timestep}).  This steady drift in the total energy
-has also been observed by Kol {\it et al.}\cite{Laird97}
+McLachlan ({\sc 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 {\sc 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
@@ -273,62 +278,62 @@ step, a 1000 {\sc ssd} particle simulation shows an av
 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 {\sc ssd} particle simulation shows an average 7\% increase in
-computation time using the {\sc dlm} 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 {\sc dlm}
-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
+step, a 1000 SSD particle simulation shows an average 7\%
+increase in computation time using the {\sc dlm} 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 {\sc dlm} 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}
-\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.}
-\label{timestep}
-\end{center}
-\end{figure}
+%\begin{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.}
+%\label{timestep}
+%\end{center}
+%\end{figure}
 
 In figure \ref{timestep}, the resulting energy drift at various time
-steps for both the {\sc dlm} and quaternion integration schemes is compared.
-All of the 1000 {\sc 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 {\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.
+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
+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
+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 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 {\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
+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$
-crystals were formed by first arranging the centers of mass of the {\sc ssd}
+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
+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
+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
@@ -346,75 +351,75 @@ 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 {\sc ssd} water model, and
-an average density versus temperature plot is shown in figure
+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{denseSSD.eps}
-\caption{Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}],
- TIP3P [Ref. \citen{Jorgensen98b}], SPC/E [Ref. \citen{Clancy94}], {\sc ssd}
- without Reaction Field, {\sc ssd}, and experiment [Ref. \citen{CRC80}]. The
- arrows indicate the change in densities observed when turning off the
- reaction field. The the lower than expected densities for the {\sc ssd}
- model were what prompted the original reparameterization of {\sc ssd1}
- [Ref. \citen{Ichiye03}].}
-\label{dense1}
-\end{center}
-\end{figure}
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{denseSSDnew.eps}
+%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
+% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
+% arrows indicate the change in densities observed when turning off the
+% reaction field. The the lower than expected densities for the SSD
+% model were what prompted the original reparameterization of SSD1
+% [Ref. \onlinecite{Ichiye03}].}
+%\label{dense1}
+%\end{center}
+%\end{figure}
 
-The density maximum for {\sc ssd} compares quite favorably to other simple
-water models. Figure \ref{dense1} also shows calculated densities of
-several other models and experiment obtained from other
+The density maximum for SSD compares quite favorably to other
+simple water models. Figure \ref{dense1} also shows calculated
+densities of several other models and experiment obtained from other
 sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
-models, {\sc ssd} has a temperature closest to the experimentally observed
-density maximum. Of the {\it charge-based} models in
+models, SSD has a temperature closest to the experimentally
+observed density maximum. Of the {\it charge-based} models in
 Fig. \ref{dense1}, TIP4P has a density maximum behavior most like that
-seen in {\sc ssd}. Though not included in this plot, it is useful
-to note that TIP5P has a density maximum nearly identical to the
+seen in SSD. Though not included in this plot, it is useful to
+note that TIP5P has a density maximum nearly identical to the
 experimentally measured temperature.
 
 It has been observed that liquid state densities in water are
 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 {\sc 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 indicate that
-there is no major benefit in choosing a longer cutoff radius in
-simulations using {\sc ssd}. This is advantageous in that the use of a
-longer cutoff radius results in a significant increase in the time
-required to obtain a single trajectory.
+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
+indicate that there is no major benefit in choosing a longer cutoff
+radius in simulations using SSD. This is advantageous in that
+the use of a longer cutoff radius results in a significant increase in
+the time required to obtain a single trajectory.
 
 The key feature to recognize in figure \ref{dense1} is the density
-scaling of {\sc ssd} relative to other common models at any given
-temperature. {\sc 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 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
+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
+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
 presence of the reaction field.\cite{Berendsen98,Nezbeda02} In order
 to test the effect of the reaction field on the density of the
 systems, the simulations were repeated without a reaction field
@@ -422,17 +427,17 @@ the curve produced from {\sc ssd} simulations using re
 \ref{dense1}. Without the reaction field, the densities increase
 to more experimentally reasonable values, especially around the
 freezing point of liquid water. The shape of the curve is similar to
-the curve produced from {\sc ssd} simulations using reaction field,
+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
-{\sc ssd}.\cite{Ichiye03} Throughout the remainder of the paper our
-reparamaterizations of {\sc ssd} will be compared with their newer {\sc ssd1}
+SSD.\cite{Ichiye03} Throughout the remainder of the paper our
+reparamaterizations of SSD will be compared with their newer SSD1
 model.
 
 \subsection{Transport Behavior}
@@ -443,43 +448,43 @@ 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
 \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
 results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} 
 
-\begin{figure}
-\begin{center}
-\epsfxsize=6in
-\epsfbox{betterDiffuse.epsi} 
-\caption{Average self-diffusion constant as a function of temperature for 
-{\sc 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
-shown, {\sc ssd} has the least deviation from the experimental values. The
-rapidly increasing diffusion constants for TIP5P and {\sc ssd} correspond to
-significant decreases in density at the higher temperatures.}
-\label{diffuse}
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{betterDiffuse.epsi} 
+%\caption{Average self-diffusion constant as a function of temperature for 
+%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
+%shown, SSD has the least deviation from the experimental values. The
+%rapidly increasing diffusion constants for TIP5P and SSD correspond to
+%significant decreases in density at the higher temperatures.}
+%\label{diffuse}
+%\end{center}
+%\end{figure} 
 
 The observed values for the diffusion constant point out one of the
-strengths of the {\sc ssd} model. Of the three models shown, the {\sc ssd} model
+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 {\sc ssd} both start off low at colder
+correct thermal trend. TIP5P and SSD both start off low at colder
 temperatures and tend to diffuse too rapidly at higher temperatures.
 This behavior at higher temperatures is not particularly surprising
-since the densities of both TIP5P and {\sc ssd} are lower than experimental
+since the densities of both TIP5P and SSD are lower than experimental
 water densities at higher temperatures.  When calculating the
-diffusion coefficients for {\sc ssd} at experimental densities (instead of
-the densities from the NPT simulations), the resulting values fall
-more in line with experiment at these temperatures.
+diffusion coefficients for SSD at experimental densities
+(instead of the densities from the NPT simulations), the resulting
+values fall more in line with experiment at these temperatures.
 
 \subsection{Structural Changes and Characterization}
 
@@ -489,32 +494,32 @@ 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{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 {\sc 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.}
-\label{contour} 
-\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
+%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.}
+%\label{contour} 
+%\end{center}
+%\end{figure} 
 
 Additional analysis of the melting process was performed using
 two-dimensional structure and dipole angle correlations. Expressions
@@ -554,38 +559,38 @@ 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
 the system densities, leading to a liquid with a more open structure
-and a density considerably lower than the already low {\sc ssd} density.  A
-better correction would be to include the quadrupole-quadrupole
-interactions for the water particles outside of the first solvation
-shell, but this would remove the simplicity and speed advantage of
-{\sc ssd}.
+and a density considerably lower than the already low SSD
+density.  A better correction would be to include the
+quadrupole-quadrupole interactions for the water particles outside of
+the first solvation shell, but this would remove the simplicity and
+speed advantage of SSD.
 
-\subsection{Adjusted Potentials: {\sc ssd/rf} and {\sc ssd/e}}
+\subsection{Adjusted Potentials: SSD/RF and SSD/E}
 
-The propensity of {\sc ssd} to adopt lower than expected densities under
+The propensity of SSD to adopt lower than expected densities under
 varying conditions is troubling, especially at higher temperatures. In
 order to correct this model for use with a reaction field, it is
 necessary to adjust the force field parameters for the primary
@@ -601,8 +606,8 @@ Sticky Dipole / Reaction Field ({\sc ssd/rf} - for use
 function cutoffs ($r_l$, $r_u$ and $r_l^\prime$, $r_u^\prime$
 respectively). The results of the reparameterizations are shown in
 table \ref{params}. We are calling these reparameterizations the Soft
-Sticky Dipole / Reaction Field ({\sc ssd/rf} - for use with a reaction
-field) and Soft Sticky Dipole Extended ({\sc ssd/e} - an attempt to improve
+Sticky Dipole / Reaction Field (SSD/RF - for use with a reaction
+field) and Soft Sticky Dipole Extended (SSD/E - an attempt to improve
 the liquid structure in simulations without a long-range correction).
 
 \begin{table}
@@ -610,8 +615,8 @@ the liquid structure in simulations without a long-ran
 \caption{Parameters for the original and adjusted models}
 \begin{tabular}{ l  c  c  c  c } 
 \hline \\[-3mm]
-\ \ \ Parameters\ \ \  & \ \ \ {\sc ssd} [Ref. \citen{Ichiye96}] \ \ \
-& \ {\sc ssd1} [Ref. \citen{Ichiye03}]\ \  & \ {\sc ssd/e}\ \  & \ {\sc ssd/rf} \\
+\ \ \ Parameters\ \ \  & \ \ \ SSD [Ref. \onlinecite{Ichiye96}] \ \ \
+& \ SSD1 [Ref. \onlinecite{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\\
@@ -627,43 +632,44 @@ the liquid structure in simulations without a long-ran
 \end{center}
 \end{table}
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=5in
-\epsfbox{GofRCompare.epsi} 
-\caption{Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with {\sc ssd/e} 
-and {\sc ssd1} without reaction field (top), as well as {\sc ssd/rf} and {\sc ssd1} with
-reaction field turned on (bottom). The insets show the respective
-first peaks in detail. Note how the changes in parameters have lowered
-and broadened the first peak of {\sc ssd/e} and {\sc ssd/rf}.}
-\label{grcompare} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=5in
+%\epsfbox{GofRCompare.epsi} 
+%\caption{Plots comparing experiment [Ref. \onlinecite{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
+%how the changes in parameters have lowered and broadened the first
+%peak of SSD/E and SSD/RF.}
+%\label{grcompare} 
+%\end{center}
+%\end{figure} 
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=6in
-\epsfbox{dualsticky_bw.eps} 
-\caption{Positive and negative isosurfaces of the sticky potential for
-{\sc ssd1} (left) and {\sc ssd/e} \& {\sc ssd/rf} (right). Light areas correspond to the
-tetrahedral attractive component, and darker areas correspond to the
-dipolar repulsive component.}
-\label{isosurface} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{dualsticky_bw.eps} 
+%\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.}
+%\label{isosurface} 
+%\end{center}
+%\end{figure} 
 
-In the original paper detailing the development of {\sc ssd}, Liu and Ichiye
+In the original paper detailing the development of SSD, Liu and Ichiye
 placed particular emphasis on an accurate description of the first
 solvation shell. This resulted in a somewhat tall and narrow first
 peak in $g(r)$ that integrated to give similar coordination numbers to
 the experimental data obtained by Soper and
 Phillips.\cite{Ichiye96,Soper86} New experimental x-ray scattering
 data from the Head-Gordon lab indicates a slightly lower and shifted
-first peak in the g$_\mathrm{OO}(r)$, so our adjustments to {\sc ssd} were
+first peak in the g$_\mathrm{OO}(r)$, so our adjustments to SSD were
 made after taking into consideration the new experimental
 findings.\cite{Head-Gordon00_1} Figure \ref{grcompare} shows the
 relocation of the first peak of the oxygen-oxygen $g(r)$ by comparing
-the revised {\sc ssd} model ({\sc ssd1}), {\sc ssd/e}, and {\sc ssd/rf} to the new
+the revised SSD model (SSD1), SSD/E, and SSD/RF to the new
 experimental results. Both modified water models have shorter peaks
 that match more closely to the experimental peak (as seen in the
 insets of figure \ref{grcompare}).  This structural alteration was
@@ -676,188 +682,190 @@ density for the overall system.  This change in intera
 see how altering the cutoffs changes the repulsive and attractive
 character of the particles. With a reduced repulsive surface (darker
 region), the particles can move closer to one another, increasing the
-density for the overall system.  This change in interaction cutoff also
-results in a more gradual orientational motion by allowing the
+density for the overall system.  This change in interaction cutoff
+also results in a more gradual orientational motion by allowing the
 particles to maintain preferred dipolar arrangements before they begin
 to feel the pull of the tetrahedral restructuring. As the particles
 move closer together, the dipolar repulsion term becomes active and
 excludes unphysical nearest-neighbor arrangements. This compares with
-how {\sc ssd} and {\sc ssd1} exclude preferred dipole alignments before the
+how SSD and SSD1 exclude preferred dipole alignments before the
 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 {\sc ssd} is a dipole based model, the
+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 {\sc 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
+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
+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 {\sc ssd/e} and {\sc ssd/rf},
+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
 and transport properties of the adapted models and compare the results
-to the original {\sc ssd} model. This comparison involved full NPT melting
-sequences for both {\sc ssd/e} and {\sc ssd/rf}, as well as NVE transport
+to the original SSD model. This comparison involved full NPT melting
+sequences for both SSD/E and SSD/RF, as well as NVE transport
 calculations at the calculated self-consistent densities. Again, the
 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
 collection times as stated previously.
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=6in
-\epsfbox{ssdeDense.epsi} 
-\caption{Comparison of densities calculated with {\sc ssd/e} to {\sc 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
-300 K with error bars included to clarify this region of
-interest. Note that both {\sc ssd1} and {\sc ssd/e} show good agreement with
-experiment when the long-range correction is neglected.}
-\label{ssdedense} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdeDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/E to 
+%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+%300 K with error bars included to clarify this region of
+%interest. Note that both SSD1 and SSD/E show good agreement with
+%experiment when the long-range correction is neglected.}
+%\label{ssdedense} 
+%\end{center}
+%\end{figure} 
 
-Fig. \ref{ssdedense} shows the density profile for the {\sc ssd/e} model
-in comparison to {\sc ssd1} without a reaction field, other common water
-models, and experimental results. The calculated densities for both
-{\sc ssd/e} and {\sc ssd1} have increased significantly over the original {\sc ssd}
-model (see fig. \ref{dense1}) and are in better agreement with the
-experimental values. At 298 K, the densities of {\sc ssd/e} and {\sc ssd1} without
+Fig. \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
+significantly over the original SSD model (see
+fig. \ref{dense1}) and are in better agreement with the experimental
+values. At 298 K, the densities of SSD/E and SSD1 without
 a long-range correction are 0.996$\pm$0.001 g/cm$^3$ and
 0.999$\pm$0.001 g/cm$^3$ respectively.  These both compare well with
 the experimental value of 0.997 g/cm$^3$, and they are considerably
-better than the {\sc 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 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 {\sc ssd/e} model has
-little effect on the melting transition. By monitoring $C_p$
-throughout these simulations, the melting transition for {\sc ssd/e} was
-shown to occur at 235 K.  The same transition temperature observed
-with {\sc ssd} and {\sc ssd1}.
+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
+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
+same transition temperature observed with SSD and SSD1.
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=6in
-\epsfbox{ssdrfDense.epsi} 
-\caption{Comparison of densities calculated with {\sc ssd/rf} to {\sc 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
-reparameterization when utilizing a reaction field long-ranged
-correction - {\sc ssd/rf} provides significantly more accurate densities
-than {\sc ssd1} when performing room temperature simulations.}
-\label{ssdrfdense} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdrfDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/RF to 
+%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+%reparameterization when utilizing a reaction field long-ranged
+%correction - SSD/RF provides significantly more accurate
+%densities than SSD1 when performing room temperature
+%simulations.}
+%\label{ssdrfdense} 
+%\end{center}
+%\end{figure} 
 
 Including the reaction field long-range correction in the simulations
 results in a more interesting comparison.  A density profile including
-{\sc ssd/rf} and {\sc ssd1} with an active reaction field is shown in figure
+SSD/RF and SSD1 with an active reaction field is shown in figure
 \ref{ssdrfdense}.  As observed in the simulations without a reaction
-field, the densities of {\sc ssd/rf} and {\sc ssd1} show a dramatic increase over
-normal {\sc ssd} (see figure \ref{dense1}). At 298 K, {\sc ssd/rf} has a density
+field, the densities of SSD/RF and SSD1 show a dramatic increase over
+normal SSD (see figure \ref{dense1}). At 298 K, SSD/RF has a density
 of 0.997$\pm$0.001 g/cm$^3$, directly in line with experiment and
-considerably better than the original {\sc ssd} value of 0.941$\pm$0.001
-g/cm$^3$ and the {\sc ssd1} value of 0.972$\pm$0.002 g/cm$^3$. These results
+considerably better than the original SSD value of 0.941$\pm$0.001
+g/cm$^3$ and the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results
 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 {\sc ssd/rf}, is identical to {\sc ssd} and only 5 K
-lower than {\sc ssd1} with a reaction field. Additionally, the difference in
-density maxima is not as extreme, with {\sc ssd/rf} showing a density
-maximum at 255 K, fairly close to the density maxima of 260 K and 265
-K, shown by {\sc ssd} and {\sc ssd1} respectively.
+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.
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=6in
-\epsfbox{ssdeDiffuse.epsi} 
-\caption{The diffusion constants calculated from {\sc ssd/e} and {\sc 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
-simulations for the respective models. {\sc ssd/e} is slightly more mobile
-than experiment at all of the temperatures, but it is closer to
-experiment at biologically relevant temperatures than {\sc ssd1} without a
-long-range correction.}
-\label{ssdediffuse} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdeDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/E and 
+%SSD1 (both without a reaction field) along with experimental results
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+%performed at the average densities observed in the 1 atm NPT
+%simulations for the respective models. SSD/E is slightly more mobile
+%than experiment at all of the temperatures, but it is closer to
+%experiment at biologically relevant temperatures than SSD1 without a
+%long-range correction.}
+%\label{ssdediffuse} 
+%\end{center}
+%\end{figure} 
 
-The reparameterization of the {\sc ssd} water model, both for use with and
+The reparameterization of the SSD water model, both for use with and
 without an applied long-range correction, brought the densities up to
 what is expected for simulating liquid water. In addition to improving
-the densities, it is important that the diffusive behavior of {\sc ssd} be
+the densities, it is important that the diffusive behavior of SSD be
 maintained or improved. Figure \ref{ssdediffuse} compares the
-temperature dependence of the diffusion constant of {\sc ssd/e} to {\sc ssd1}
+temperature dependence of the diffusion constant of SSD/E to SSD1
 without an active reaction field at the densities calculated from
 their respective NPT simulations at 1 atm. The diffusion constant for
-{\sc ssd/e} is consistently higher than experiment, while {\sc ssd1} remains lower
+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 {\sc 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 {\sc ssd1} and {\sc ssd/e} show this
+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 {\sc ssd1}
+implications on the diffusive behavior of the models.  While SSD1
 shows more experimentally accurate diffusive behavior in the high
-temperature regimes, {\sc ssd/e} shows more accurate behavior in the
+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 {\sc ssd/e} relative to {\sc ssd1} under the most commonly simulated
+of SSD/E relative to SSD1 under the most commonly simulated
 conditions.
 
-\begin{figure}
-\begin{center} 
-\epsfxsize=6in
-\epsfbox{ssdrfDiffuse.epsi} 
-\caption{The diffusion constants calculated from {\sc ssd/rf} and {\sc 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 the 1 atm NPT
-simulations for both of the models. {\sc ssd/rf} simulates the diffusion of
-water throughout this temperature range very well. The rapidly
-increasing diffusion constants at high temperatures for both models
-can be attributed to lower calculated densities than those observed in
-experiment.}
-\label{ssdrfdiffuse} 
-\end{center}
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdrfDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/RF and 
+%SSD1 (both with an active reaction field) along with
+%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+%NVE calculations were performed at the average densities observed in
+%the 1 atm NPT simulations for both of the models. SSD/RF
+%simulates the diffusion of water throughout this temperature range
+%very well. The rapidly increasing diffusion constants at high
+%temperatures for both models can be attributed to lower calculated
+%densities than those observed in experiment.}
+%\label{ssdrfdiffuse} 
+%\end{center}
+%\end{figure} 
 
-In figure \ref{ssdrfdiffuse}, the diffusion constants for {\sc ssd/rf} are
-compared to {\sc ssd1} with an active reaction field. Note that {\sc ssd/rf}
+In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
+compared to SSD1 with an active reaction field. Note that SSD/RF
 tracks the experimental results quantitatively, identical within error
 throughout most of the temperature range shown and exhibiting only a
-slight increasing trend at higher temperatures. {\sc ssd1} tends to diffuse
+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. {\sc ssd/rf} does not suffer from this problem as much as {\sc ssd1}
+temperatures. SSD/RF does not suffer from this problem as much as SSD1
 because the calculated densities are closer to the experimental
 values. These results again emphasize the importance of careful
 reparameterization when using an altered long-range correction.
@@ -867,27 +875,25 @@ experimental data at ambient conditions} 
 \renewcommand{\thefootnote}{\thempfootnote}
 \begin{center}
 \caption{Properties of the single-point water models compared with
-experimental data at ambient conditions} 
+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]
-\ \ \ \ \ \  & \ \ \ {\sc ssd1} \ \ \ & \ {\sc ssd/e} \ \ \ & \ {\sc ssd1} (RF) \ \ 
-\ & \ {\sc 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 &
-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 &
+\ \ $\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. \onlinecite{Head-Gordon00_1}} \\ 
+\ \ 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
-Ref. \citen{Krynicki66}}
+Ref. \onlinecite{Soper86}}  \\
+\ \ $\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. \onlinecite{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. \onlinecite{Krynicki66}}
 \end{tabular}
 \label{liquidproperties}
 \end{center}
@@ -909,15 +915,16 @@ but the coordination numbers of {\sc ssd/e} and {\sc s
 $a$ and $b$ are the radial locations of the minima following the first
 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 {\sc ssd/e} and {\sc ssd/rf} show an improvement
-over {\sc ssd1}.  This improvement is primarily due to extension of the
-first solvation shell in the new parameter sets.  Because $n_H$ and
-$n_C$ are nearly identical in {\sc ssd1}, it appears that all molecules in
-the first solvation shell are involved in hydrogen bonds.  Since $n_H$
-and $n_C$ differ in the newly parameterized models, the orientations
-in the first solvation shell are a bit more ``fluid''.  Therefore {\sc ssd1}
-overstructures the first solvation shell and our reparameterizations
-have returned this shell to more realistic liquid-like behavior.
+but the coordination numbers of SSD/E and SSD/RF show an
+improvement over SSD1.  This improvement is primarily due to
+extension of the first solvation shell in the new parameter sets.
+Because $n_H$ and $n_C$ are nearly identical in SSD1, it appears
+that all molecules in the first solvation shell are involved in
+hydrogen bonds.  Since $n_H$ and $n_C$ differ in the newly
+parameterized models, the orientations in the first solvation shell
+are a bit more ``fluid''.  Therefore SSD1 overstructures the
+first solvation shell and our reparameterizations have returned this
+shell to more realistic liquid-like behavior.
 
 The time constants for the orientational autocorrelation functions
 are also displayed in Table \ref{liquidproperties}. The dipolar
@@ -935,44 +942,45 @@ the NMR data in Ref. \citen{Krynicki66} at a temperatu
 averaged over five detailed NVE simulations performed at the ambient
 conditions for each of the respective models. It should be noted that
 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
+the NMR data in Ref. \onlinecite{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}.
-Again, {\sc ssd/e} and {\sc ssd/rf} show improved behavior over {\sc ssd1}, both with
+fashion described in Ref. \onlinecite{Krynicki66}. Similarly, $\tau_1$ was
+recomputed for 298~K from the data in Ref. \onlinecite{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 {\sc ssd1}; however, these results
+leads to much improved time constants for SSD1; however, these results
 also include a corresponding decrease in system density.
-Orientational relaxation times published in the original {\sc ssd} dynamics
+Orientational relaxation times published in the original SSD dynamics
 papers are smaller than the values observed here, and this difference
 can be attributed to the use of the Ewald sum.\cite{Ichiye99}
 
 \subsection{Additional Observations}
 
-\begin{figure}
-\begin{center}
-\epsfxsize=6in
-\epsfbox{icei_bw.eps} 
-\caption{The most stable crystal structure assumed by the {\sc 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.}
-\label{weirdice} 
-\end{center}
-\end{figure}
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{icei_bw.eps} 
+%\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.}
+%\label{weirdice} 
+%\end{center}
+%\end{figure}
 
 While performing a series of melting simulations on an early iteration
-of {\sc 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 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 forms of ice, we have denoted the structure
+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
+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
+forms of ice, we have denoted the structure
 Ice-$\sqrt{\smash[b]{-\text{I}}}$ (Ice-{\it i}).  A large enough
 portion of the sample crystallized that we have been able to obtain a
 near ideal crystal structure of Ice-{\it i}. Figure \ref{weirdice}
@@ -987,63 +995,62 @@ structure for at least the {\sc ssd/e} model. To verif
 favorable enough to stabilize an entire crystal generated around them.
 
 Initial simulations indicated that Ice-{\it i} is the preferred ice
-structure for at least the {\sc ssd/e} model. To verify this, a
-comparison was made between near ideal crystals of ice~$I_h$,
-ice~$I_c$, and Ice-{\it i} at constant pressure with {\sc ssd/e}, {\sc
-ssd/rf}, and {\sc ssd1}. Near-ideal versions of the three types of
-crystals were cooled to 1 K, and enthalpies of formation of each were
-compared using all three water models.  Enthalpies were estimated from
-the isobaric-isothermal simulations using $H=U+P_{\text ext}V$ where
-$P_{\text ext}$ is the applied pressure.  A constant value of
--60.158 kcal / mol has been added to place our zero for the
-enthalpies of formation for these systems at the traditional state
-(elemental forms at standard temperature and pressure).  With every
-model in the {\sc ssd} family, Ice-{\it i} had the lowest calculated
-enthalpy of formation.
+structure for at least the SSD/E model. To verify this, a comparison
+was made between near ideal crystals of ice~$I_h$, ice~$I_c$, and
+Ice-{\it i} at constant pressure with SSD/E, SSD/RF, and
+SSD1. Near-ideal versions of the three types of crystals were cooled
+to 1 K, and enthalpies of formation of each were compared using all
+three water models.  Enthalpies were estimated from the
+isobaric-isothermal simulations using $H=U+P_{\text ext}V$ where
+$P_{\text ext}$ is the applied pressure.  A constant value of -60.158
+kcal / mol has been added to place our zero for the enthalpies of
+formation for these systems at the traditional state (elemental forms
+at standard temperature and pressure).  With every model in the SSD
+family, Ice-{\it i} had the lowest calculated enthalpy of formation.
 
 \begin{table}
 \begin{center}
 \caption{Enthalpies of Formation (in kcal / mol) of the three crystal
-structures (at 1 K) exhibited by the {\sc ssd} family of water models} 
+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]
-\ \ \ {\sc ssd/e} & -12.286 & -12.292 & -13.590 \\
-\ \ \ {\sc ssd/rf} & -12.935 & -12.917 & -14.022 \\
-\ \ \ {\sc ssd1} & -12.496 & -12.411 & -13.417 \\
-\ \ \ {\sc ssd1} (RF) & -12.504 & -12.411 & -13.134 \\
+\ \ \ SSD/E & -72.444 & -72.450 & -73.748 \\
+\ \ \ SSD/RF & -73.093 & -73.075 & -74.180 \\
+\ \ \ SSD1 & -72.654 & -72.569 & -73.575 \\
+\ \ \ SSD1 (RF) & -72.662 & -72.569 & -73.292 \\
 \end{tabular}
 \label{iceenthalpy}
 \end{center}
 \end{table}
 
 In addition to these energetic comparisons, melting simulations were
-performed with ice-{\it i} as the initial configuration using {\sc ssd/e},
-{\sc ssd/rf}, and {\sc ssd1} both with and without a reaction field. The melting
-transitions for both {\sc ssd/e} and {\sc ssd1} without reaction field occurred at
-temperature in excess of 375~K.  {\sc ssd/rf} and {\sc ssd1} with a reaction field
+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
 showed more reasonable melting transitions near 325~K.  These melting
-point observations clearly show that all of the {\sc ssd}-derived models
+point observations clearly show that all of the SSD-derived models
 prefer the ice-{\it i} structure.
 
 \section{Conclusions}
 
 The density maximum and temperature dependence of the self-diffusion
-constant were studied for the {\sc ssd} water model, both with and without
-the use of reaction field, via a series of NPT and NVE
+constant were studied for the SSD water model, both with and
+without the use of reaction field, via a series of NPT and NVE
 simulations. The constant pressure simulations showed a density
 maximum near 260 K. In most cases, the calculated densities were
 significantly lower than the densities obtained from other water
-models (and experiment). Analysis of self-diffusion showed {\sc ssd} to
-capture the transport properties of water well in both the liquid and
-supercooled liquid regimes.
+models (and experiment). Analysis of self-diffusion showed SSD
+to capture the transport properties of water well in both the liquid
+and supercooled liquid regimes.
 
-In order to correct the density behavior, the original {\sc ssd} model was
-reparameterized for use both with and without a reaction field ({\sc ssd/rf}
-and {\sc ssd/e}), and comparisons were made with {\sc ssd1}, Ichiye's density
-corrected version of {\sc ssd}. Both models improve the liquid structure,
+In order to correct the density behavior, the original SSD model was
+reparameterized for use both with and without a reaction field (SSD/RF
+and SSD/E), and comparisons were made with SSD1, Ichiye's density
+corrected version of SSD. Both models improve the liquid structure,
 densities, and diffusive properties under their respective simulation
 conditions, indicating the necessity of reparameterization when
 changing the method of calculating long-range electrostatic
@@ -1052,9 +1059,9 @@ by the {\sc ssd} family of water models is somewhat tr
 simulations of biochemical systems.
 
 The existence of a novel low-density ice structure that is preferred
-by the {\sc ssd} family of water models is somewhat troubling, since liquid
-simulations on this family of water models at room temperature are
-effectively simulations of supercooled or metastable liquids.  One
+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 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
@@ -1078,8 +1085,286 @@ DMR-0079647.
 \newpage
 
 \bibliographystyle{jcp}
-\bibliography{nptSSD} 
+\bibliography{nptSSD}
 
-%\pagebreak 
+\newpage
 
+\begin{list}
+  {Figure \arabic{captions}: }{\usecounter{captions}
+	\setlength{\rightmargin}{\leftmargin}}
+	
+\item 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.
+
+\item Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E
+[Ref. \onlinecite{Clancy94}], SSD without Reaction Field, SSD, and
+experiment [Ref. \onlinecite{CRC80}]. The arrows indicate the change
+in densities observed when turning off the reaction field. The the
+lower than expected densities for the SSD model were what prompted the
+original reparameterization of SSD1 [Ref. \onlinecite{Ichiye03}].
+
+\item Average self-diffusion constant as a function of temperature for 
+SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three
+water models shown, SSD has the least deviation from the experimental
+values. The rapidly increasing diffusion constants for TIP5P and SSD
+correspond to significant decreases in density at the higher
+temperatures.
+
+\item An illustration of angles involved in the correlations observed in
+Fig. \ref{contour}.
+
+\item 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.
+
+\item Plots comparing experiment [Ref. \onlinecite{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 how the changes in parameters
+have lowered and broadened the first peak of SSD/E and SSD/RF.
+
+\item 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.
+
+\item Comparison of densities calculated with SSD/E to 
+SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+300 K with error bars included to clarify this region of
+interest. Note that both SSD1 and SSD/E show good agreement with
+experiment when the long-range correction is neglected.
+
+\item Comparison of densities calculated with SSD/RF to 
+SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+reparameterization when utilizing a reaction field long-ranged
+correction - SSD/RF provides significantly more accurate
+densities than SSD1 when performing room temperature
+simulations.
+
+\item The diffusion constants calculated from SSD/E and 
+SSD1 (both without a reaction field) along with experimental results
+[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+performed at the average densities observed in the 1 atm NPT
+simulations for the respective models. SSD/E is slightly more mobile
+than experiment at all of the temperatures, but it is closer to
+experiment at biologically relevant temperatures than SSD1 without a
+long-range correction.
+
+\item The diffusion constants calculated from SSD/RF and 
+SSD1 (both with an active reaction field) along with
+experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+NVE calculations were performed at the average densities observed in
+the 1 atm NPT simulations for both of the models. SSD/RF
+simulates the diffusion of water throughout this temperature range
+very well. The rapidly increasing diffusion constants at high
+temperatures for both models can be attributed to lower calculated
+densities than those observed in experiment.
+
+\item 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.
+\end{list}
+
+\newpage
+
+\begin{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.}
+\label{timestep}
+\end{center}
+\end{figure}
+
+\newpage
+
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{denseSSDnew.eps}
+%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
+% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
+% arrows indicate the change in densities observed when turning off the
+% reaction field. The the lower than expected densities for the SSD
+% model were what prompted the original reparameterization of SSD1
+% [Ref. \onlinecite{Ichiye03}].}
+\label{dense1}
+\end{center}
+\end{figure}
+
+\newpage
+
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{betterDiffuse.epsi} 
+%\caption{Average self-diffusion constant as a function of temperature for 
+%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
+%shown, SSD has the least deviation from the experimental values. The
+%rapidly increasing diffusion constants for TIP5P and SSD correspond to
+%significant decreases in density at the higher temperatures.}
+\label{diffuse}
+\end{center}
+\end{figure}
+
+\newpage
+
+\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}
+
+\newpage
+
+\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
+%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.}
+\label{contour} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{GofRCompare.epsi} 
+%\caption{Plots comparing experiment [Ref. \onlinecite{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
+%how the changes in parameters have lowered and broadened the first
+%peak of SSD/E and SSD/RF.}
+\label{grcompare} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=7in
+\epsfbox{dualsticky_bw.eps} 
+%\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.}
+\label{isosurface} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdeDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/E to 
+%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+%300 K with error bars included to clarify this region of
+%interest. Note that both SSD1 and SSD/E show good agreement with
+%experiment when the long-range correction is neglected.}
+\label{ssdedense} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdrfDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/RF to 
+%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+%reparameterization when utilizing a reaction field long-ranged
+%correction - SSD/RF provides significantly more accurate
+%densities than SSD1 when performing room temperature
+%simulations.}
+\label{ssdrfdense} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+ \begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdeDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/E and 
+%SSD1 (both without a reaction field) along with experimental results
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+%performed at the average densities observed in the 1 atm NPT
+%simulations for the respective models. SSD/E is slightly more mobile
+%than experiment at all of the temperatures, but it is closer to
+%experiment at biologically relevant temperatures than SSD1 without a
+%long-range correction.}
+\label{ssdediffuse} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdrfDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/RF and 
+%SSD1 (both with an active reaction field) along with
+%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+%NVE calculations were performed at the average densities observed in
+%the 1 atm NPT simulations for both of the models. SSD/RF
+%simulates the diffusion of water throughout this temperature range
+%very well. The rapidly increasing diffusion constants at high
+%temperatures for both models can be attributed to lower calculated
+%densities than those observed in experiment.}
+\label{ssdrfdiffuse} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{icei_bw.eps} 
+%\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.}
+\label{weirdice} 
+\end{center}
+\end{figure} 
+
 \end{document}