--- trunk/ssdePaper/nptSSD.tex	2003/09/03 21:36:09	743
+++ trunk/ssdePaper/nptSSD.tex	2004/02/05 20:47:50	1029
@@ -1,57 +1,64 @@
-\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
-%\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
+%\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
+\documentclass[11pt]{article}
+\usepackage{endfloat}
 \usepackage{amsmath}
+\usepackage{epsf}
+\usepackage{berkeley}
+\usepackage{setspace}
+\usepackage{tabularx}
 \usepackage{graphicx}
-
-%\usepackage{endfloat}
+\usepackage[ref]{overcite}
 %\usepackage{berkeley}
-%\usepackage{epsf}
-%\usepackage[ref]{overcite}
-%\usepackage{setspace}
-%\usepackage{tabularx}
-%\usepackage{graphicx}
 %\usepackage{curves}
-%\usepackage{amsmath}
-%\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
 
 \begin{document}
 
-\title{On the temperature dependent structural and transport properties of the soft sticky dipole (SSD) and related single point water models}
+\title{On the structural and transport properties of the soft sticky
+dipole ({\sc ssd}) and related single point water models}
 
-\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
-\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
-
-\address{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
+
 \begin{abstract}
-NVE and NPT molecular dynamics simulations were performed in order to
-investigate the density maximum and temperature dependent transport
-for the SSD water model, both with and without the use of reaction
-field. The constant pressure simulations of the melting of both $I_h$
-and $I_c$ ice showed a density maximum near 260 K. In most cases, the
-calculated densities were significantly lower than the densities
-calculated in simulations of other water models. Analysis of particle
-diffusion showed SSD to capture the transport properties of
-experimental very well in both the normal and super-cooled liquid
-regimes. In order to correct the density behavior, SSD was
-reparameterized for use both with and without a long-range interaction
-correction, SSD/RF and SSD/E respectively. In addition to correcting
-the abnormally low densities, these new versions were show to maintain
-or improve upon the transport and structural features of the original
-water model.
+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.
+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
+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
+experimental water very well in both the normal and super-cooled
+liquid regimes.  We also present our re-parameterized versions of {\sc 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}
+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}).
+Additionally, a novel low-density ice structure is presented
+which appears to be the most stable ice structure for the entire {\sc ssd}
+family.
 \end{abstract}
 
-\maketitle
+\newpage
 
 %\narrowtext 
 
@@ -62,232 +69,268 @@ One of the most important tasks in simulations of bioc
 
 \section{Introduction}
 
-One of the most important tasks in simulations of biochemical systems
-is the proper depiction of water and water solvation. In fact, the
-bulk of the calculations performed in solvated simulations are of
-interactions with or between solvent molecules. Thus, the outcomes of
-these types of simulations are highly dependent on the physical
-properties of water, both as individual molecules and in
-groups/bulk. Due to the fact that explicit solvent accounts for a
-massive portion of the calculations, it necessary to simplify the
-solvent to some extent in order to complete simulations in a
-reasonable amount of time. In the case of simulating water in
-bio-molecular studies, the balance between accurate properties and
-computational efficiency is especially delicate, and it has resulted
-in a variety of different water
-models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these models
-get specific properties correct or better than their predecessors, but
-this is often at a cost of some other properties or of computer
-time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P succeeds
-in improving the structural and transport properties over its
-predecessors, yet this comes at a greater than 50\% increase in
-computational cost.\cite{Jorgensen01,Jorgensen00} One recently
-developed model that succeeds in both retaining accuracy of system
-properties and simplifying calculations to increase computational
-efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96}
+One of the most important tasks in the simulation of biochemical
+systems is the proper depiction of the aqueous environment of the
+molecules of interest.  In some cases (such as in the simulation of 
+phospholipid bilayers), the majority of the calculations that are
+performed involve interactions with or between solvent molecules.
+Thus, the properties one may observe in biochemical simulations are
+going to be highly dependent on the physical properties of the water
+model that is chosen.
 
-The Soft Sticky Dipole (SSD)\ water 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} SSD
-consists of a single point dipole with a Lennard-Jones core and a
-sticky potential that directs the particles to assume the proper
-hydrogen bond orientation in the first solvation shell. Thus, the
-interaction between two SSD water molecules \emph{i} and \emph{j} is
-given by the potential
+There is an especially delicate balance between computational
+efficiency and the ability of the water model to accurately predict
+the properties of bulk
+water.\cite{Jorgensen83,Berendsen87,Jorgensen00} For example, the
+TIP5P model improves on the structural and transport properties of
+water relative to the previous TIP models, yet this comes at a greater
+than 50\% increase in computational
+cost.\cite{Jorgensen01,Jorgensen00}
+
+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
+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
+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 interaction between two {\sc 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}
-(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
+({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)\ +
 u_{ij}^{sp}
-(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
+({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j),
 \end{equation}
-where the $\mathbf{r}_{ij}$ is the position vector between molecules
-\emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and
-$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
-orientations of the respective molecules. The Lennard-Jones, dipole,
-and sticky parts of the potential are giving by the following
-equations,
+where the ${\bf r}_{ij}$ is the position vector between molecules
+\emph{i} and \emph{j} with magnitude $r_{ij}$, and
+${\bf \Omega}_i$ and ${\bf \Omega}_j$ represent the orientations of
+the two molecules. The Lennard-Jones and dipole interactions are given
+by the following familiar forms:
 \begin{equation}
-u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
+u_{ij}^{LJ}(r_{ij}) = 4\epsilon
+\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]
+\ ,
 \end{equation}
+and
 \begin{equation}
-u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ ,
+u_{ij}^{dp} = \frac{|\mu_i||\mu_j|}{4 \pi \epsilon_0 r_{ij}^3} \left(
+\hat{\bf u}_i \cdot \hat{\bf u}_j - 3(\hat{\bf u}_i\cdot\hat{\bf
+r}_{ij})(\hat{\bf u}_j\cdot\hat{\bf r}_{ij}) \right)\ ,
 \end{equation}
+where $\hat{\bf u}_i$ and $\hat{\bf u}_j$ are the unit vectors along
+the dipoles of molecules $i$ and $j$ respectively. $|\mu_i|$ and
+$|\mu_j|$ are the strengths of the dipole moments, and $\hat{\bf
+r}_{ij}$ is the unit vector pointing from molecule $j$ to molecule
+$i$.
+
+The sticky potential is somewhat less familiar:
 \begin{equation}
-\begin{split}
 u_{ij}^{sp}
-(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) 
-&=
-\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
-& \quad \ +
-s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
-\end{split}
+({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
+\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)
++ s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf
+\Omega}_j)]\ .
+\label{stickyfunction}
 \end{equation}
-where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
-unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
-$\nu_0$ scales the strength of the overall sticky potential, $s$ and
-$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
-functions take the following forms,
+Here, $\nu_0$ is a strength parameter for the sticky potential, and
+$s$ and $s^\prime$ are cubic switching functions which turn off the
+sticky interaction beyond the first solvation shell. The $w$ function
+can be thought of as an attractive potential with tetrahedral
+geometry:
 \begin{equation}
-w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
+w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
 \end{equation}
+while the $w^\prime$ function counters the normal aligned and
+anti-aligned structures favored by point dipoles:
 \begin{equation}
-w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
+w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^\circ,
 \end{equation}
-where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
-term that promotes hydrogen bonding orientations within the first
-solvation shell, and $w^\prime$ is a dipolar repulsion term that
-repels unrealistic dipolar arrangements within the first solvation
-shell. A more detailed description of the functional parts and
-variables in this potential can be found in other
-articles.\cite{Ichiye96,Ichiye99}
+It should be noted that $w$ is proportional to the sum of the $Y_3^2$
+and $Y_3^{-2}$ spherical harmonics (a linear combination which
+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}
+articles.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03}
 
-Being that this is a one-site point dipole model, the actual force
-calculations are simplified significantly. In the original Monte Carlo
-simulations using this model, Ichiye \emph{et al.} reported a
-calculation speed up of up to an order of magnitude over other
-comparable models while maintaining the structural behavior of
-water.\cite{Ichiye96} In the original molecular dynamics studies of
-SSD, it was shown that it actually improves upon the prediction of
-water's dynamical properties 3 and 4-point models.\cite{Ichiye99} This
-attractive combination of speed and accurate depiction of solvent
-properties makes SSD a model of interest for the simulation of large
-scale biological systems, such as membrane phase behavior, a specific
-interest within our group.
+Since {\sc 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
+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
+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
+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
+simulation of large scale biochemical simulations.
 
-Up to this point, a detailed look at the model's structure and ion
-solvation abilities has been performed.\cite{Ichiye96} In addition, a
-thorough investigation of the dynamic properties of SSD was performed
-by Chandra and Ichiye focusing on translational and orientational
-properties at 298 K.\cite{Ichiye99} This study focuses on determining
-the density maximum for SSD utilizing both microcanonical and
-isobaric-isothermal ensemble molecular dynamics, while using the
-reaction field method for handling long-ranged dipolar interactions. A
-reaction field method has been previously implemented in Monte Carlo
-simulations by Liu and Ichiye in order to study the static dielectric
-constant for the model.\cite{Ichiye96b} This paper will expand the
-scope of these original simulations to look on how the reaction field
-affects the physical and dynamic properties of SSD systems.
+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.
 
 \section{Methods}
 
-As stated previously, in this study the long-range dipole-dipole
-interactions were accounted for using the reaction field method. The
-magnitude of the reaction field acting on dipole \emph{i} is given by
+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
 \begin{equation}
 \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
-\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} \boldsymbol{\mu}_{j} f(r_{ij})\  ,
+\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}),
 \label{rfequation}
 \end{equation}
 where $\mathcal{R}$ is the cavity defined by the cutoff radius
 ($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
-system (80 in this case), $\boldsymbol{\mu}_{j}$ is the dipole moment
-vector of particle \emph{j}, and $f(r_{ij})$ is a cubic switching
+system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole
+moment vector of particle $j$, and $s(r_{ij})$ is a cubic switching
 function.\cite{AllenTildesley} The reaction field contribution to the
-total energy by particle \emph{i} is given by
-$-\frac{1}{2}\boldsymbol{\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque
-on dipole \emph{i} by
-$\boldsymbol{\mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley} Use
-of reaction field is known to alter the orientational dynamic
-properties, such as the dielectric relaxation time, based on changes
-in the length of the cutoff radius.\cite{Berendsen98} This variable
-behavior makes reaction field a less attractive method than other
-methods, like the Ewald summation; however, for the simulation of
-large-scale system, the computational cost benefit of reaction field
-is dramatic. To address some of the dynamical property alterations due
-to the use of reaction field, simulations were also performed without
-a surrounding dielectric and suggestions are proposed on how to make
-SSD more compatible with a reaction field.
+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.
  
-Simulations were performed in both the isobaric-isothermal and
-microcanonical ensembles. The constant pressure simulations were
+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.
+
+Simulations to obtain the preferred densities of the models were
+performed in the isobaric-isothermal (NPT) ensemble, while all
+dynamical properties were obtained from microcanonical (NVE)
+simulations done at densities matching the NPT density for a
+particular target temperature.  The constant pressure simulations were
 implemented using an integral thermostat and barostat as outlined by
-Hoover.\cite{Hoover85,Hoover86} For the constant pressure
-simulations, the \emph{Q} parameter for the was set to 5.0 amu
-\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at
-100 ps.
+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.
 
 Integration of the equations of motion was carried out using the
-symplectic splitting method proposed by Dullweber \emph{et
-al.}.\cite{Dullweber1997} The reason for this integrator selection
-deals with poor energy conservation of rigid body systems using
-quaternions. While quaternions work well for orientational motion in
-alternate ensembles, the microcanonical ensemble has a constant energy
-requirement that is actually quite sensitive to errors in the
-equations of motion. The original implementation of this code utilized
-quaternions for rotational motion propagation; however, a detailed
-investigation showed that they resulted in a steady drift in the total
-energy, something that has been observed by others.\cite{Laird97}
+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}
 
 The key difference in the integration method proposed by Dullweber
 \emph{et al.} is that the entire rotation matrix is propagated from
-one time step to the next. In the past, this would not have been as
-feasible a option, being that the rotation matrix for a single body is
-nine elements long as opposed to 3 or 4 elements for Euler angles and
-quaternions respectively. System memory has become much less of an
-issue in recent times, and this has resulted in substantial benefits
-in energy conservation. There is still the issue of an additional 5 or
-6 additional elements for describing the orientation of each particle,
-which will increase dump files substantially. Simply translating the
-rotation matrix into its component Euler angles or quaternions for
-storage purposes relieves this burden.
+one time step to the next.  The additional memory required by the
+algorithm is inconsequential on modern computers, and translating the
+rotation matrix into quaternions for storage purposes makes trajectory
+data quite compact.
 
-The symplectic splitting method allows for Verlet style integration of
-both linear and angular motion of rigid bodies. In the integration
-method, the orientational propagation involves a sequence of matrix
-evaluations to update the rotation matrix.\cite{Dullweber1997} These
-matrix rotations end up being more costly computationally than the
-simpler arithmetic quaternion propagation. On average, a 1000 SSD
-particle simulation shows a 7\% increase in simulation time using the
-symplectic step method in place of quaternions. This cost is more than
-justified when comparing the energy conservation of the two methods as
-illustrated in figure \ref{timestep}.
+The {\sc dlm} method allows for Verlet style integration of both
+translational and orientational motion of rigid bodies. In this
+integration method, the orientational propagation involves a sequence
+of matrix evaluations to update the rotation
+matrix.\cite{Dullweber1997} These matrix rotations are more costly
+than the simpler arithmetic quaternion propagation. With the same time
+step, a 1000 {\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
+\ref{timestep}.
 
 \begin{figure}
-\includegraphics[width=61mm, angle=-90]{timeStep.epsi}
-\caption{Energy conservation using quaternion based integration versus 
-the symplectic step method proposed by Dullweber \emph{et al.} with
-increasing time step. For each time step, the dotted line is total
-energy using the symplectic step integrator, and the solid line comes
-from the quaternion integrator. The larger time step plots are shifted
-up from the true energy baseline for clarity.}
+\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 symplectic step and quaternion integration schemes
-is compared. All of the 1000 SSD particle simulations started with the
-same configuration, and the only difference was the method for
-handling rotational motion. At time steps of 0.1 and 0.5 fs, both
-methods for propagating particle rotation conserve energy fairly well,
-with the quaternion method showing a slight energy drift over time in
-the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
-energy conservation benefits of the symplectic step method are clearly
-demonstrated.
-
-Energy drift in these SSD particle simulations was unnoticeable for
-time steps up to three femtoseconds. A slight energy drift on the
-order of 0.012 kcal/mol per nanosecond was observed at a time step of
-four femtoseconds, and as expected, this drift increases dramatically
-with increasing time step. To insure accuracy in the constant energy
-simulations, time steps were set at 2 fs and kept at this value for
-constant pressure simulations as well.
+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.
 
-Ice crystals in both the $I_h$ and $I_c$ lattices were generated as
-starting points for all the simulations. The $I_h$ crystals were
-formed by first arranging the center of masses 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 a
-rectangular shape with a edge length ratio of approximately
+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}
+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 down to 25 K. The particles
-were then allowed translate with fixed orientations at a constant
+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
+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
 isotropic volume fluctuations were performed under constant pressure,
@@ -297,608 +340,744 @@ constant pressure and temperature dynamics. This invol
 \section{Results and discussion}
 
 Melting studies were performed on the randomized ice crystals using
-constant pressure and temperature dynamics. This involved an initial
-randomization of velocities about the starting temperature of 25 K for
-varying amounts of time. The systems were all equilibrated for 100 ps
-prior to a 200 ps data collection run at each temperature setting,
-ranging 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 then 5 K. The above
-equilibration and production times were sufficient in that the system
-volume fluctuations dampened out in all but the very cold simulations
-(below 225 K). In order to further improve statistics, five separate
-simulation progressions were performed, and the averaged results from
-the $I_h$ melting simulations are shown in figure \ref{dense1}.
+isobaric-isothermal (NPT) dynamics. During melting simulations, the
+melting transition and the density maximum can both be observed,
+provided that the density maximum occurs in the liquid and not the
+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
+sufficient in that fluctuations in the volume autocorrelation function
+were damped out in all simulations in under 20 ps.
 
-\begin{figure}
-\includegraphics[width=65mm, angle=-90]{1hdense.epsi}
-\caption{Average density of SSD water at increasing temperatures 
-starting from ice $I_h$ lattice.}
-\label{dense1}
-\end{figure}
-
 \subsection{Density Behavior}
-In the initial average density versus temperature plot, the density
-maximum clearly appears between 255 and 265 K. The calculated
-densities within this range were nearly indistinguishable, as can be
-seen in the zoom of this region of interest, shown in figure
-\ref{dense1}. The greater certainty of the average value at 260 K makes 
-a good argument for the actual density maximum residing at this
-midpoint value. Figure \ref{dense1} was constructed using ice $I_h$
-crystals for the initial configuration; and 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). In addition, the $I_c$ crystals are
-more fragile than the $I_h$ crystals, leading them to deform into a
-dense glassy state at lower temperatures. This resulted in an overall
-low temperature density maximum at 200 K, but they still retained a
-common liquid state density maximum with the $I_h$ simulations.
 
+Our initial simulations focused on the original {\sc 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
+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). 
+
 \begin{figure}
-\includegraphics[width=65mm,angle=-90]{dense2.eps}
-\caption{Density versus temperature for TIP4P\cite{Jorgensen98b}, 
-TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction
-Field, SSD, and Experiment\cite{CRC80}. }
-\label{dense2}
+\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}
 
-The density maximum for SSD actually compares quite favorably to other
-simple water models. Figure \ref{dense2} shows a plot of these
-findings with the density progression of several other models and
-experiment obtained from other
+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
 sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
-models, SSD has results closest to the experimentally observed water
-density maximum. Of the listed water models, TIP4P has a density
-maximum behavior most like that seen in SSD. Though not shown, it is
-useful to note that TIP5P has a water density maximum nearly identical
-to experiment.
+models, {\sc 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
+experimentally measured temperature.
 
-Possibly of more importance is the density scaling of SSD relative to
-other common models at any given temperature (Fig. \ref{dense2}). Note
-that the SSD model 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 with the use of a
-reaction field for long-range electrostatic interactions, so the most
-likely reason for these significantly lower densities in these
-simulations is the presence of the reaction field.\cite{Berendsen98}
-In order to test the effect of the reaction field on the density of
-the systems, the simulations were repeated for the temperature region
-of interest without a reaction field present. The results of these
-simulations are also displayed in figure \ref{dense2}. Without
-reaction field, these densities increase considerably 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 SSD simulations using reaction field, specifically the
-rapidly decreasing densities at higher temperatures; however, a slight
-shift in the density maximum location, down to 245 K, is
-observed. This is probably a more accurate comparison to the other
-listed water models in that no long range corrections were applied in
-those simulations.\cite{Clancy94,Jorgensen98b}
+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.
 
-It has been observed that densities are dependent on the cutoff radius
-used for a variety of water models in simulations both with and
-without the use of reaction field.\cite{Berendsen98} In order to
-address the possible affect of cutoff radius, simulations were
-performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the
-previous SSD simulations, all performed with a cutoff of 9.0 \AA. All
-the resulting densities overlapped within error and showed no
-significant trend in lower or higher densities as a function of cutoff
-radius, both for simulations 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 comforting in that the
-use of a longer cutoff radius results in a near doubling of the time
-required to compute 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
+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
+present. The results of these simulations are also displayed in figure
+\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,
+specifically the rapidly decreasing densities at higher temperatures;
+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
+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}
+model.
 
 \subsection{Transport Behavior}
-Of importance in these types of studies are the transport properties
-of the particles and how they change when altering the environmental
-conditions. In order to probe transport, constant energy simulations
-were performed about the average density uncovered by the constant
-pressure simulations. Simulations started with randomized velocities
-and underwent 50 ps of temperature scaling and 50 ps of constant
-energy equilibration before obtaining a 200 ps trajectory. Diffusion
-constants were calculated via root-mean square deviation analysis. The
-averaged results from 5 sets of these NVE simulations is displayed in
-figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
-results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} 
+
+Accurate dynamical properties of a water model are particularly
+important when using the model to study permeation or transport across
+biological membranes.  In order to probe transport in bulk water,
+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
+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}
-\includegraphics[width=65mm, angle=-90]{betterDiffuse.epsi} 
-\caption{Average diffusion coefficient over increasing temperature for 
-SSD, SPC/E\cite{Clancy94}, TIP5P\cite{Jorgensen01}, and Experimental
-data from Gillen \emph{et al.}\cite{Gillen72}, and from
-Mills\cite{Mills73}.}
+\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} 
 
 The observed values for the diffusion constant point out one of the
-strengths of the SSD model. Of the three experimental models shown,
-the SSD model has the most accurate depiction of the diffusion trend
-seen in experiment in both the supercooled and normal regimes. SPC/E
-does a respectable job by getting similar values as SSD and experiment
-around 290 K; however, it deviates at both higher and lower
-temperatures, failing to predict the experimental trend. TIP5P and SSD
-both start off low at the colder temperatures and tend to diffuse too
-rapidly at the higher temperatures. This type of trend at the higher
-temperatures is not surprising in that the densities of both TIP5P and
-SSD are lower than experimental water at temperatures higher than room
-temperature. When calculating the diffusion coefficients for SSD at
-experimental densities, the resulting values fall more in line with
-experiment at these temperatures, albeit not at standard
-pressure. Results under these conditions can be found later in this
-paper.
+strengths of the {\sc ssd} model. Of the three models shown, the {\sc 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
+deviates at both higher and lower temperatures, failing to predict the
+correct thermal trend. TIP5P and {\sc 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
+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.
 
 \subsection{Structural Changes and Characterization}
+
 By starting the simulations from the crystalline state, the melting
-transition and the ice structure can be studied along with the liquid
-phase behavior beyond the melting point. To locate the melting
-transition, the constant pressure heat capacity (C$_\text{p}$) was
-monitored 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 these ice crystals. When the reaction field is turned off,
-the melting transition occurs at 235 K.  These melting transitions are
-considerably lower than the experimental value, but this is not
-surprising in that SSD is a simple rigid body model with a fixed
-dipole.
+transition and the ice structure can be obtained along with the liquid
+phase behavior beyond the melting point. The constant pressure heat
+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
+these ice crystals. When the reaction field is turned off, the melting
+transition occurs at 235 K.  These melting transitions are
+considerably lower than the experimental value.
 
-\begin{figure} 
-\includegraphics[width=85mm]{fullContours.eps} 
-\caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at 
-100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for
-clarity: dark areas signify peaks while light areas signify
-depressions. White areas have g(\emph{r}) values below 0.5 and black
-areas have values above 1.5.}
-\label{contour} 
-\end{figure} 
-
-Additional analyses for understanding the melting phase-transition
-process were performed via two-dimensional structure and dipole angle
-correlations. Expressions for the correlations are as follows:
-
 \begin{figure}
-\includegraphics[width=45mm]{corrDiag.eps} 
-\caption{Two dimensional illustration of the angles involved in the 
-correlations observed in figure \ref{contour}.}
+\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{multline}
-g_{\text{AB}}(r,\cos\theta) = \\
-\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
-\end{multline}
-\begin{multline}
-g_{\text{AB}}(r,\cos\omega) = \\
-\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
-\end{multline}
-where $\theta$ and $\omega$ refer to the angles shown in the above
-illustration. By binning over both distance and the cosine of the
-desired angle between the two dipoles, the g(\emph{r}) can be
-dissected to determine the common dipole arrangements that constitute
-the peaks and troughs. Frames A and B of figure \ref{contour} show a
-relatively crystalline state of an ice $I_c$ simulation. The first
-peak of the g(\emph{r}) primarily consists of the preferred hydrogen
-bonding arrangements as dictated by the tetrahedral sticky potential,
-one peak for the donating and the other for the accepting hydrogen
-bonds. Due to the high degree of crystallinity of the sample, the
-second and third solvation shells show a repeated peak arrangement
+\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} 
+
+Additional analysis of the melting process was performed using
+two-dimensional structure and dipole angle correlations. Expressions
+for these correlations are as follows:
+
+\begin{equation}
+g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
+\end{equation}
+\begin{equation}
+g_{\text{AB}}(r,\cos\omega) = 
+\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
+\end{equation}
+where $\theta$ and $\omega$ refer to the angles shown in figure
+\ref{corrAngle}. By binning over both distance and the cosine of the
+desired angle between the two dipoles, the $g(r)$ can be analyzed to
+determine the common dipole arrangements that constitute the peaks and
+troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
+of figure \ref{contour} show results from an ice $I_c$ simulation. The
+first peak in the $g(r)$ consists primarily of the preferred hydrogen
+bonding arrangements as dictated by the tetrahedral sticky potential -
+one peak for the hydrogen bond donor and the other for the hydrogen
+bond acceptor.  Due to the high degree of crystallinity of the sample,
+the second and third solvation shells show a repeated peak arrangement
 which decays at distances around the fourth solvation shell, near the
 imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
-In the higher temperature simulation shown in frames C and D, the
-repeated peak features are significantly blurred. The first solvation
-shell still shows the strong effect of the sticky-potential, although
-it covers a larger area, extending to include a fraction of aligned
+In the higher temperature simulation shown in frames C and D, these
+long-range features deteriorate rapidly. The first solvation shell
+still shows the strong effect of the sticky-potential, although it
+covers a larger area, extending to include a fraction of aligned
 dipole peaks within the first solvation shell. The latter peaks lose
-definition as thermal motion and the competing dipole force overcomes
-the sticky potential's tight tetrahedral structuring of the fluid.
+due to thermal motion and as the competing dipole force overcomes the
+sticky potential's tight tetrahedral structuring of the crystal.
 
 This complex interplay between dipole and sticky interactions was
 remarked upon as a possible reason for the split second peak in the
-oxygen-oxygen g(\emph{r}).\cite{Ichiye96} At low temperatures, the
-second solvation shell peak appears to have two distinct parts that
+oxygen-oxygen pair correlation function,
+$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
-by equatorial anti-parallel dipole orientations, and there is tightly
-bunched group of axially arranged dipoles that most likely consist of
-the smaller fraction aligned dipole pairs. The trailing part of the
-split peak at 5 \AA\ is dominated by aligned dipoles that range
-primarily within the axial to the chief 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, with the primary
-energetically favorable dipole arrangements populating the region
-immediately outside of it's range (around 4 \AA), and arrangements
-that seek to ideally satisfy both the sticky and dipole forces locate
-themselves just beyond this region (around 5 \AA).
+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
+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
+dipole forces locate themselves just beyond this initial buildup
+(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
-leading of the two peaks can be pushed out and merged with the outer
-split peak just by extending the switching function cutoff
-($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even
-5 \AA. This type of correction is not recommended for improving the
-liquid structure, because the second solvation shell will 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 normal
-SSD behavior shown previously. A better correction would be to include
-the quadrupole-quadrupole interactions for the water particles outside
-of the first solvation shell, but this reduces the simplicity and
-speed advantage of SSD, so it is not the most desirable path to take.
+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
+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}.
 
-\subsection{Adjusted Potentials: SSD/E and SSD/RF}
-The propensity of SSD to adopt lower than expected densities under
+\subsection{Adjusted Potentials: {\sc ssd/rf} and {\sc ssd/e}}
+
+The propensity of {\sc ssd} to adopt lower than expected densities under
 varying conditions is troubling, especially at higher temperatures. In
-order to correct this behavior, it's necessary to adjust the force
-field parameters for the primary intermolecular interactions. In
-undergoing a reparameterization, it is important not to focus on just
-one property and neglect the other important properties. In this case,
-it would be ideal to correct the densities while maintaining the
-accurate transport properties.
+order to correct this model for use with a reaction field, it is
+necessary to adjust the force field parameters for the primary
+intermolecular interactions. In undergoing a reparameterization, it is
+important not to focus on just one property and neglect the other
+important properties. In this case, it would be ideal to correct the
+densities while maintaining the accurate transport behavior.
 
-The possible parameters for tuning include the $\sigma$ and $\epsilon$
-Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
-attractive and dipole repulsive terms with their respective
-cutoffs. To alter the attractive and repulsive terms of the sticky
-potential independently, it is necessary to separate the terms as
-follows:
-\begin{equation}
-\begin{split}
-u_{ij}^{sp}
-(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) &=
-\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\\
-& \quad \ + \frac{\nu_0^\prime}{2}
-[s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)],
-\end{split}
-\end{equation}
+The parameters available for tuning include the $\sigma$ and
+$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the
+strength of the sticky potential ($\nu_0$), and the cutoff distances
+for the sticky attractive and dipole repulsive cubic switching
+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
+the liquid structure in simulations without a long-range correction).
 
-where $\nu_0$ scales the strength of the tetrahedral attraction and
-$\nu_0^\prime$ acts in an identical fashion on the dipole repulsion
-term. For purposes of the reparameterization, the separation was
-performed, but the final parameters were adjusted so that it is
-unnecessary to separate the terms when implementing the adjusted water
-potentials. The results of the reparameterizations are shown in table
-\ref{params}. Note that both the tetrahedral attractive and dipolar 
-repulsive don't share the same lower cutoff ($r_l$) in the newly
-parameterized potentials - soft sticky dipole enhanced (SSD/E) and
-soft sticky dipole reaction field (SSD/RF).
-
 \begin{table}
+\begin{center}
 \caption{Parameters for the original and adjusted models}
-\begin{tabular}{ l  c  c  c } 
+\begin{tabular}{ l  c  c  c  c } 
 \hline \\[-3mm]
-\ Parameters & \ \ \  SSD$^\dagger$\ \ \ \  & \ SSD/E\ \  & \ SSD/RF\ \ \\
+\ \ \ Parameters\ \ \  & \ \ \ {\sc ssd} [Ref. \citen{Ichiye96}] \ \ \
+& \ {\sc ssd1} [Ref. \citen{Ichiye03}]\ \  & \ {\sc ssd/e}\ \  & \ {\sc ssd/rf} \\
 \hline \\[-3mm]
-\ \ \ $\sigma$ (\AA)  & 3.051 & 3.035 & 3.019\\
-\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\
-\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\
-\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
-\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\
-\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\
-\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
-\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\
-\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\
-\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96}
+\ \ \ $\sigma$ (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
+\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
+\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\
+\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
+\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\
+\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
+\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
+\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\
+\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
 \end{tabular}
 \label{params}
+\end{center}
 \end{table}
 
-\begin{figure} 
-\includegraphics[width=85mm]{gofrCompare.epsi} 
-\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E 
-and SSD without reaction field (top), as well as SSD/RF and SSD with
+\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. Solid Line - experiment, dashed line - SSD/E
-and SSD/RF, and dotted line - SSD (with and without reaction field).}
+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} 
-\includegraphics[width=85mm]{dualsticky.ps} 
-\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \& 
-SSD/RF (right). Light areas correspond to the tetrahedral attractive
-part, and the darker areas correspond to the dipolar repulsive part.}
+\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} 
 
-In the 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 sharp first peak 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
-adjustments to SSD were made while 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(\emph{r}) by comparing the original SSD (with and without reaction
-field), SSD-E, and SSD-RF to the new experimental results. Both the
-modified water models have shorter peaks that are brought in more
-closely to the experimental peak (as seen in the insets of figure
-\ref{grcompare}). This structural alteration was accomplished by a
-reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
-of the sticky potential strength and cutoffs. The cutoffs for the
-tetrahedral attractive and dipolar repulsive terms were nearly swapped
-with each other. Isosurfaces of the original and modified sticky
-potentials are shown in figure \cite{isosurface}. In these
-isosurfaces, it is easy to see how altering the cutoffs changes the
-repulsive and attractive character of the particles. With a reduced
-repulsive surface (the 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 particles to maintain preferred
-dipolar arrangements before they begin to feel the pull of the
-tetrahedral restructuring. Upon moving closer together, the dipolar
-repulsion term becomes active and excludes the unphysical
-arrangements. This compares with the original SSD's excluding dipolar
-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
-improves the densities considerably by allowing the persistence of
-full dipolar character below the previous 4.0 \AA\ cutoff.
+In the original paper detailing the development of {\sc 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
+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
+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
+accomplished by the combined reduction in the Lennard-Jones $\sigma$
+variable and adjustment of the sticky potential strength and cutoffs.
+As can be seen in table \ref{params}, the cutoffs for the tetrahedral
+attractive and dipolar repulsive terms were nearly swapped with each
+other.  Isosurfaces of the original and modified sticky potentials are
+shown in figure \ref{isosurface}. In these isosurfaces, it is easy to
+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
+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
+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\
+cutoff.
 
-While adjusting the location and shape of the first peak of
-g(\emph{r}) improves the densities to some degree, these changes alone
-are insufficient to bring the system densities up to the values
-observed experimentally. To finish bringing up the densities, the
-dipole moments were increased in both the adjusted models. Being 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 improve 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, so there is quite a range available for
-adjusting the dipole
-moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of
-the dipole moments to 2.418 and 2.48 D for SSD/E and SSD/RF
-respectively is moderate in the range of the experimental values;
-however, it leads to significant changes in the density and transport
-of the water models.
+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
+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
+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 {\sc ssd/e} and {\sc ssd/rf},
+respectively, leads to significant changes in the density and
+transport of the water models.
 
-In order to demonstrate the benefits of this reparameterization, a
+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 SSD model. This comparison involved full NPT melting
-sequences for both SSD/E and SSD/RF, as well as NVE transport
-calculations at both self-consistent and experimental
-densities. Again, the results come from five separate simulations of
-1024 particle systems, and the melting sequences were started from
-different ice $I_h$ crystals constructed as stated previously. Like
-before, all of the NPT simulations were equilibrated for 100 ps before
-a 200 ps data collection run, and they used the previous temperature's
-final configuration as a starting point. All of the NVE simulations
-had the same thermalization, equilibration, and data collection times
-stated earlier in this paper.
+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
+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
+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} 
-\includegraphics[width=85mm]{ssdecompare.epsi} 
-\caption{Comparison of densities calculated with SSD/E to SSD without a 
-reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
-SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
-includes error bars, and the calculated results from the other
-references were removed for clarity.}
+\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} 
 
-Figure \ref{ssdedense} shows the density profile for the SSD/E water
-model in comparison to the original SSD without a reaction field,
-experiment, and the other common water models considered
-previously. The calculated densities have increased significantly over
-the original SSD model and match the experimental value just below 298
-K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
-compares well with the experimental value of 0.997 g/cm$^3$ and is
-considerably better than the SSD value of 0.967$\pm$0.003
-g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
-out the curve at higher temperatures, only the improvement is marginal
-at best. This steep drop in densities is due to the dipolar rather
-than charge based interactions which decay more rapidly at longer
-distances.
- 
-By monitoring C$\text{p}$ throughout these simulations, the melting
-transition for SSD/E was observed at 230 K, about 5 degrees lower than
-SSD. The resulting density maximum is located at 240 K, again about 5
-degrees lower than the SSD value of 245 K. Though there is a decrease
-in both of these values, the corrected densities near room temperature
-justify the modifications taken.
+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
+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}.
 
-\begin{figure} 
-\includegraphics[width=85mm]{ssdrfcompare.epsi} 
-\caption{Comparison of densities calculated with SSD/RF to SSD with a 
-reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
-SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
-includes error bars, and the calculated results from the other
-references were removed for clarity.}
+\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} 
 
-Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
-SSD with an active reaction field. Like in the simulations of SSD/E,
-the densities show a dramatic increase over normal SSD. At 298 K,
-SSD/RF has a density of 0.997$\pm$0.001 g/cm$^3$, right in line with
-experiment and considerably better than the SSD value of
-0.941$\pm$0.001 g/cm$^3$. The melting point is observed at 240 K,
-which is 5 degrees lower than SSD with a reaction field, and the
-density maximum at 255 K, again 5 degrees lower than SSD. The density
-at higher temperature still drops off more rapidly than the charge
-based models but is in better agreement than SSD/E.
+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
+\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
+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
+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.
 
-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 particle transport be maintained
-or improved. Figure \ref{ssdediffuse} compares the temperature
-dependence of the diffusion constant of SSD/E to SSD without an active
-reaction field, both at the densities calculated at 1 atm and at the
-experimentally calculated densities for super-cooled and liquid
-water. In the upper plot, the diffusion constant for SSD/E is
-consistently a little faster than experiment, while SSD starts off
-slower than experiment and crosses to merge with SSD/E at high
-temperatures. Both models follow the experimental trend well, but
-diffuse too rapidly at higher temperatures. This abnormally fast
-diffusion is caused by the decreased system density. Since the
-densities of SSD/E don't deviate as much from experiment as those of
-SSD, it follows the experimental trend more closely. This observation
-is backed up by looking at the lower plot. The diffusion constants for
-SSD/E track with the experimental values while SSD deviates on the low
-side of the trend with increasing temperature. This is again a product
-of SSD/E having densities closer to experiment, and not deviating to
-lower densities with increasing temperature as rapidly.
-
-\begin{figure} 
-\includegraphics[width=85mm]{ssdediffuse.epsi} 
-\caption{Plots of the diffusion constants calculated from SSD/E and SSD,
- both without a reaction field along with experimental results from
- Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
- upper plot is at densities calculated from the NPT simulations at a
- pressure of 1 atm, while the lower plot is at the experimentally
- calculated densities.}
+\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} 
-\includegraphics[width=85mm]{ssdrfdiffuse.epsi} 
-\caption{Plots of the diffusion constants calculated from SSD/RF and SSD,
- both with an active reaction field along with experimental results
- from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
- upper plot is at densities calculated from the NPT simulations at a
- pressure of 1 atm, while the lower plot is at the experimentally
- calculated densities.}
+The reparameterization of the {\sc 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
+maintained or improved. Figure \ref{ssdediffuse} compares the
+temperature dependence of the diffusion constant of {\sc ssd/e} to {\sc 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
+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
+the experimental values is caused by the rapidly decreasing system
+density with increasing temperature.  Both {\sc ssd1} and {\sc ssd/e} show this
+deviation in particle mobility, but this trend has different
+implications on the diffusive behavior of the models.  While {\sc ssd1}
+shows more experimentally accurate diffusive behavior in the high
+temperature regimes, {\sc 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
+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} 
 
-In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
-compared with SSD with an active reaction field. In the upper plot,
-SSD/RF tracks with the experimental results incredibly well, identical
-within error throughout the temperature range and only showing a
-slight increasing trend at higher temperatures. SSD also tracks
-experiment well, only it tends to diffuse a little more slowly at low
-temperatures and deviates to diffuse too rapidly at high
-temperatures. As was stated in the SSD/E comparisons, this deviation
-away from the ideal trend is due to a rapid decrease in density at
-higher temperatures. SSD/RF doesn't suffer from this problem as much
-as SSD, because the calculated densities are more true to
-experiment. This is again emphasized in the lower plot, where SSD/RF
-tracks the experimental diffusion exactly while SSD's diffusion
-constants are slightly too low due to its need for a lower density at
-the specified temperature.
+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}
+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
+more slowly at low temperatures and deviates to diffuse too rapidly at
+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}
+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.
 
-\subsection{Additional Observations}
+\begin{table}
+\begin{minipage}{\linewidth}
+\renewcommand{\thefootnote}{\thempfootnote}
+\begin{center}
+\caption{Properties of the single-point water models compared with
+experimental data at ambient conditions} 
+\begin{tabular}{ l  c  c  c  c  c } 
+\hline \\[-3mm]
+\ \ \ \ \ \  & \ \ \ {\sc ssd1} \ \ \ & \ {\sc ssd/e} \ \ \ & \ {\sc ssd1} (RF) \ \ 
+\ & \ {\sc 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 &
+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}}
+\end{tabular}
+\label{liquidproperties}
+\end{center}
+\end{minipage}
+\end{table}
 
-While performing the melting sequences of SSD/E, some interesting
-observations were made. After melting at 230 K, two of the systems
-underwent crystallization events near 245 K. As the heating process
-continued, the two systems remained crystalline until finally melting
-between 320 and 330 K. These simulations were excluded from the data
-set shown in figure \ref{ssdedense} and replaced with two additional
-melting sequences that did not undergo this anomalous phase
-transition, while this crystallization event was investigated
-separately.
+Table \ref{liquidproperties} gives a synopsis of the liquid state
+properties of the water models compared in this study along with the
+experimental values for liquid water at ambient conditions. The
+coordination number ($n_C$) and number of hydrogen bonds per particle
+($n_H$) were calculated by integrating the following relations:
+\begin{equation}
+n_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr,
+\end{equation}
+\begin{equation}
+n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr,
+\end{equation}
+where $\rho$ is the number density of the specified pair interactions,
+$a$ and $b$ are the radial locations of the minima following the first
+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.
 
-\begin{figure}
-\includegraphics[width=85mm]{povIce.ps} 
-\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
-\label{weirdice} 
-\end{figure}
+The time constants for the orientational autocorrelation functions
+are also displayed in Table \ref{liquidproperties}. The dipolar
+orientational time correlation functions ($C_{l}$) are described
+by:
+\begin{equation}
+C_{l}(t) = \langle P_l[\hat{\mathbf{u}}_j(0)\cdot\hat{\mathbf{u}}_j(t)]\rangle,
+\end{equation}
+where $P_l$ are Legendre polynomials of order $l$ and
+$\hat{\mathbf{u}}_j$ is the unit vector pointing along the molecular
+dipole.\cite{Rahman71} From these correlation functions, the
+orientational relaxation time of the dipole vector can be calculated
+from an exponential fit in the long-time regime ($t >
+\tau_l$).\cite{Rothschild84} Calculation of these time constants were
+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
+34$^\circ$C.\cite{Rahman71} Because of the strong temperature
+dependence of $\tau_2$, it is necessary to recalculate it at 298 K to
+make proper comparisons. The value shown in Table
+\ref{liquidproperties} was calculated from the same NMR data in the 
+fashion described in 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
+and without an active reaction field. Turning on the reaction field
+leads to much improved time constants for {\sc ssd1}; however, these results
+also include a corresponding decrease in system density.
+Orientational relaxation times published in the original {\sc 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}
 
-The final configurations of these two melting sequences shows an
-expanded zeolite-like crystal structure that does not correspond to
-any known form of ice. For convenience and to help distinguish it from
-the experimentally observed forms of ice, this crystal structure will
-henceforth be referred to as ice-zero (ice 0). The crystallinity was
-extensive enough than a near ideal crystal structure could be
-obtained. Figure \ref{weirdice} shows the repeating crystal structure
-of a typical crystal at 5 K. The unit cell contains eight molecules,
-and figure \ref{unitcell} shows a unit cell built from the water
-particle center of masses that can be used to construct a repeating
-lattice, similar to figure \ref{weirdice}. Each molecule is hydrogen
-bonded to four other water molecules; however, the hydrogen bonds are
-flexed rather than perfectly straight. This results in a skewed
-tetrahedral geometry about the central molecule. Looking back at
-figure \ref{isosurface}, it is easy to see how these flexed hydrogen
-bonds are allowed in that the attractive regions are conical in shape,
-with the greatest attraction in the central region. Though not ideal,
-these flexed hydrogen bonds are favorable enough to stabilize an
-entire crystal generated around them. In fact, the imperfect ice 0
-crystals were so stable that they melted at greater than room
-temperature.
+\subsection{Additional Observations}
 
 \begin{figure}
-\includegraphics[width=65mm]{ice0cell.eps} 
-\caption{Simple unit cell for constructing ice 0. In this cell, $c$ is 
-equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
-\label{unitcell} 
+\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}
 
-The initial simulations indicated that ice 0 is the preferred ice
-structure for at least SSD/E. To verify this, a comparison was made
-between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at
-constant pressure with SSD/E, SSD/RF, and SSD. Near ideal versions of
-the three types of crystals were cooled to ~1 K, and the potential
-energies of each were compared using all three water models. With
-every water model, ice 0 turned out to have the lowest potential
-energy: 5\% lower than $I_h$ with SSD, 6.5\% lower with SSD/E, and
-7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
-was observed to be ~2\% less stable than ice $I_h$. In addition to
-having the lowest potential energy, ice 0 was the most expanded of the
-three ice crystals, ~5\% less dense than ice $I_h$ with all of the
-water models. In all three water models, ice $I_c$ was observed to be
-~2\% more dense than ice $I_h$. 
+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
+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}
+shows the repeating crystal structure of a typical crystal at 5
+K. Each water molecule is hydrogen bonded to four others; however, the
+hydrogen bonds are bent rather than perfectly straight. This results
+in a skewed tetrahedral geometry about the central molecule.  In
+figure \ref{isosurface}, it is apparent that these flexed hydrogen
+bonds are allowed due to the conical shape of the attractive regions,
+with the greatest attraction along the direct hydrogen bond
+configuration. Though not ideal, these flexed hydrogen bonds are
+favorable enough to stabilize an entire crystal generated around them.
 
-In addition to the low temperature comparisons, melting sequences were
-performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
-and SSD both with and without a reaction field. The melting
-transitions for both SSD/E and SSD without a reaction field occurred
-at temperature in excess of 375 K. SSD/RF and SSD with a reaction
-field had more reasonable melting transitions, down near 325 K. These
-melting point observations emphasize how preferred this crystal
-structure is over the most common types of ice when using these single
-point water models.
+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.
 
-Recognizing that the above tests show ice 0 to be both the most stable
-and lowest density crystal structure for these single point water
-models, it is interesting to speculate on the favorability of this
-crystal structure with the different charge based models. As a quick
-test, these 3 crystal types were converted from SSD type particles to
-TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
-minimizations were performed on all of these crystals to compare the
-system energies. Again, ice 0 was observed to have the lowest total
-system energy. The total energy of ice 0 was ~2\% lower than ice
-$I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial
-results, we would not be surprised if results from the other common
-water models show ice 0 to be the lowest energy crystal structure. A
-continuation on work studing ice 0 with multipoint water models will
-be published in a coming article.
+\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} 
+\begin{tabular}{ l  c  c  c  } 
+\hline \\[-3mm]
+\ \ \ 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 \\
+\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
+showed more reasonable melting transitions near 325~K.  These melting
+point observations clearly show that all of the {\sc ssd}-derived models
+prefer the ice-{\it i} structure.
+
 \section{Conclusions}
-The density maximum and temperature dependent transport for the SSD
-water model, both with and without the use of reaction field, were
-studied via a series of NPT and NVE simulations. The constant pressure
-simulations of the melting of both $I_h$ and $I_c$ ice showed a
-density maximum near 260 K. In most cases, the calculated densities
-were significantly lower than the densities calculated in simulations
-of other water models. Analysis of particle diffusion showed SSD to
-capture the transport properties of experimental very well in both the
-normal and super-cooled liquid regimes. In order to correct the
-density behavior, SSD was reparameterized for use both with and
-without a long-range interaction correction, SSD/RF and SSD/E
-respectively. In addition to correcting the abnormally low densities,
-these new versions were show to maintain or improve upon the transport
-and structural features of the original water model, all while
-maintaining the fast performance of the SSD water model. This work
-shows these simple water models, and in particular SSD/E and SSD/RF,
-to be excellent choices to represent explicit water in future
+
+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
+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.
+
+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,
+densities, and diffusive properties under their respective simulation
+conditions, indicating the necessity of reparameterization when
+changing the method of calculating long-range electrostatic
+interactions.  In general, however, these simple water models are
+excellent choices for representing explicit water in large scale
 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
+way to destabilize this unphysical ice structure would be to make the
+range of angles preferred by the attractive part of the sticky
+potential much narrower.  This would require extensive
+reparameterization to maintain the same level of agreement with the
+experiments.
+
+Additionally, our initial calculations show that the Ice-{\it i}
+structure may also be a preferred crystal structure for at least one
+other popular multi-point water model (TIP3P), and that much of the
+simulation work being done using this popular model could also be at
+risk for crystallization into this unphysical structure.  A future
+publication will detail the relative stability of the known ice
+structures for a wide range of popular water models.
+
 \section{Acknowledgments}
-The authors would like to thank the National Science Foundation for
-funding under grant CHE-0134881. Computation time was provided by the
-Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR
-00 79647.
+Support for this project was provided by the National Science
+Foundation under grant CHE-0134881. Computation time was provided by
+the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
+DMR-0079647.
 
-\bibliographystyle{jcp}
+\newpage
 
+\bibliographystyle{jcp}
 \bibliography{nptSSD} 
 
 %\pagebreak