--- trunk/ssdePaper/nptSSD.tex	2003/09/10 15:38:43	757
+++ trunk/ssdePaper/nptSSD.tex	2004/02/06 19:15:44	1033
@@ -1,17 +1,15 @@
-\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
-%\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
+%\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
+\documentclass[preprint,aps]{revtex4}
+%\documentclass[11pt]{article}
+%\usepackage{endfloat}
 \usepackage{amsmath}
+\usepackage{epsf}
+\usepackage{berkeley}
+\usepackage{setspace}
+\usepackage{tabularx}
 \usepackage{graphicx}
-
-%\usepackage{endfloat}
-%\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
@@ -20,39 +18,52 @@
 %\brokenpenalty=10000
 %\renewcommand{\baselinestretch}{1.2}
 %\renewcommand\citemid{\ } % no comma in optional reference note
+\newcounter{captions}
+\newcounter{figs}
 
 \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 (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}}
+\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
 
-\address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
+\affiliation{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
 Notre Dame, Indiana 46556}
 
 \date{\today}
 
+
 \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 (SSD) water
+model and two related reparameterizations of this single-point model.
+A combination of microcanonical and isobaric-isothermal molecular
+dynamics simulations were used to calculate these properties, both
+with and without the use of reaction field to handle long-range
+electrostatics.  The isobaric-isothermal (NPT) simulations of the
+melting of both ice-$I_h$ and ice-$I_c$ showed a density maximum near
+260~K.  In most cases, the use of the reaction field resulted in
+calculated densities which were were significantly lower than
+experimental densities.  Analysis of self-diffusion constants shows
+that the original SSD model captures the transport properties of
+experimental water very well in both the normal and super-cooled
+liquid regimes.  We also present our reparameterized versions of SSD
+for use both with the reaction field or without any long-range
+electrostatic corrections.  These are called the SSD/RF and SSD/E
+models respectively.  These modified models were shown to maintain or
+improve upon the experimental agreement with the structural and
+transport properties that can be obtained with either the original SSD
+or the density corrected version of the original model (SSD1).
+Additionally, a novel low-density ice structure is presented
+which appears to be the most stable ice structure for the entire SSD
+family.
 \end{abstract}
 
 \maketitle
 
+\newpage
+
 %\narrowtext 
 
 
@@ -62,239 +73,269 @@ 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 (SSD) water
+model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The SSD model
+was developed by Ichiye \emph{et al.} as a modified form of the
+hard-sphere water model proposed by Bratko, Blum, and
+Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model
+which has an interaction site that is both a point dipole and a
+Lennard-Jones core.  However, since the normal aligned and
+anti-aligned geometries favored by point dipoles are poor mimics of
+local structure in liquid water, a short ranged ``sticky'' potential
+is also added.  The sticky potential directs the molecules to assume
+the proper hydrogen bond orientation in the first solvation shell.
+
+The interaction between two SSD water molecules \emph{i} and \emph{j}
+is given by the potential
 \begin{equation}
 u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
-(\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 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 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, Liu and Ichiye reported that using SSD
+decreased computer time by a factor of 6-7 compared to other
+models.\cite{Ichiye96} What is most impressive is that this savings
+did not come at the expense of accurate depiction of the liquid state
+properties.  Indeed, SSD maintains reasonable agreement with the Soper
+data for the structural features of liquid
+water.\cite{Soper86,Ichiye96} Additionally, the dynamical properties
+exhibited by SSD agree with experiment better than those of more
+computationally expensive models (like TIP3P and
+SPC/E).\cite{Ichiye99} The combination of speed and accurate depiction
+of solvent properties makes SSD a very attractive model for the
+simulation of large scale biochemical simulations.
 
-One of the key limitations of this water model, however, is that it
-has been parameterized for use with the Ewald Sum technique for the
-handling of long-ranged interactions.  When studying very large
-systems, the Ewald summation and even particle-mesh Ewald become
-computational burdens with their respective ideal $N^\frac{3}{2}$ and
-$N\log N$ calculation scaling orders for $N$ particles.\cite{Darden99}
-
-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 SSD model is that it was parameterized for
+use with the Ewald sum to handle long-range interactions.  This would
+normally be the best way of handling long-range interactions in
+systems that contain other point charges.  However, our group has
+recently become interested in systems with point dipoles as mimics for
+neutral, but polarized regions on molecules (e.g. the zwitterionic
+head group regions of phospholipids).  If the system of interest does
+not contain point charges, the Ewald sum and even particle-mesh Ewald
+become computational bottlenecks.  Their respective ideal
+$N^\frac{3}{2}$ and $N\log N$ calculation scaling orders for $N$
+particles can become prohibitive when $N$ becomes
+large.\cite{Darden99} In applying this water model in these types of
+systems, it would be useful to know its properties and behavior under
+the more computationally efficient reaction field (RF) technique, or
+even with a simple cutoff. This study addresses these issues by
+looking at the structural and transport behavior of SSD over a
+variety of temperatures with the purpose of utilizing the RF
+correction technique.  We then suggest modifications to the parameters
+that result in more realistic bulk phase behavior.  It should be noted
+that in a recent publication, some of the original investigators of
+the SSD water model have suggested adjustments to the SSD
+water model to address abnormal density behavior (also observed here),
+calling the corrected model SSD1.\cite{Ichiye03} In what
+follows, we compare our reparamaterization of SSD with both the
+original SSD and SSD1 models with the goal of improving
+the bulk phase behavior of an SSD-derived model in simulations
+utilizing the reaction field.
 
 \section{Methods}
 
-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 technique or by resorting to a
+simple cubic switching function at a cutoff radius.  One of the early
+applications of a reaction field was actually in Monte Carlo
+simulations of liquid water.\cite{Barker73} Under this method, the
+magnitude of the reaction field acting on dipole $i$ is
 \begin{equation}
 \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
-\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} \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 of simulated
+water, and there is particular sensitivity of these properties on
+changes in the length of the cutoff radius.\cite{Berendsen98} This
+variable behavior makes reaction field a less attractive method than
+the Ewald sum.  However, for very large systems, the computational
+benefit of reaction field is dramatic.
  
-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 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 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 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.}
-\label{timestep}
-\end{figure}
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{timeStep.epsi}
+%\caption{Energy conservation using both quaternion-based integration and
+%the {\sc dlm} method with increasing time step. The larger time step
+%plots are shifted from the true energy baseline (that of $\Delta t$ =
+%0.1~fs) for clarity.}
+%\label{timestep}
+%\end{center}
+%\end{figure}
 
 In figure \ref{timestep}, the resulting energy drift at various time
-steps for both the 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.
+steps for both the {\sc dlm} and quaternion integration schemes is
+compared.  All of the 1000 SSD particle simulations started with
+the same configuration, and the only difference was the method used to
+handle orientational motion. At time steps of 0.1 and 0.5~fs, both
+methods for propagating the orientational degrees of freedom conserve
+energy fairly well, with the quaternion method showing a slight energy
+drift over time in the 0.5~fs time step simulation. At time steps of 1
+and 2~fs, the energy conservation benefits of the {\sc dlm} method are
+clearly demonstrated. Thus, while maintaining the same degree of
+energy conservation, one can take considerably longer time steps,
+leading to an overall reduction in computation time.
 
-Energy drift in 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
+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.
 
-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
+Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
+were generated as starting points for all simulations. The $I_h$
+crystals were formed by first arranging the centers of mass of the SSD
+particles into a ``hexagonal'' ice lattice of 1024 particles. Because
+of the crystal structure of $I_h$ ice, the simulation box assumed an
+orthorhombic shape with an edge length ratio of approximately
 1.00$\times$1.06$\times$1.23. The particles were then allowed to
 orient freely about fixed positions with angular momenta randomized at
-400 K for varying times. The rotational temperature was then scaled
-down in stages to slowly cool the crystals down to 25 K. The particles
-were then allowed 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
+400~K for varying times. The rotational temperature was then scaled
+down in stages to slowly cool the crystals to 25~K. The particles were
+then allowed to translate with fixed orientations at a constant
+pressure of 1 atm for 50~ps at 25~K. Finally, all constraints were
+removed and the ice crystals were allowed to equilibrate for 50~ps at
+25~K and a constant pressure of 1~atm.  This procedure resulted in
 structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
-rules\cite{Bernal33,Rahman72}.  This method was also utilized in the
+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,
@@ -304,610 +345,1026 @@ 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.
 
-\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}
-\end{figure}
+Our initial simulations focused on the original SSD water model,
+and an average density versus temperature plot is shown in figure
+\ref{dense1}. Note that the density maximum when using a reaction
+field appears between 255 and 265~K.  There were smaller fluctuations
+in the density at 260~K than at either 255 or 265~K, so we report this
+value as the location of the density maximum. Figure \ref{dense1} was
+constructed using ice $I_h$ crystals for the initial configuration;
+though not pictured, the simulations starting from ice $I_c$ crystal
+configurations showed similar results, with a liquid-phase density
+maximum in this same region (between 255 and 260~K). 
 
-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
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{denseSSDnew.eps}
+%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
+% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
+% arrows indicate the change in densities observed when turning off the
+% reaction field. The the lower than expected densities for the SSD
+% model were what prompted the original reparameterization of SSD1
+% [Ref. \onlinecite{Ichiye03}].}
+%\label{dense1}
+%\end{center}
+%\end{figure}
+
+The density maximum for 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, 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 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
+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 SSD
+simulations, all performed with a cutoff of 9.0~\AA. All of the
+resulting densities overlapped within error and showed no significant
+trend toward lower or higher densities as a function of cutoff radius,
+for simulations both with and without reaction field. These results
+indicate that there is no major benefit in choosing a longer cutoff
+radius in simulations using SSD. This is advantageous in that
+the use of a longer cutoff radius results in a significant increase in
+the time required to obtain a single trajectory.
+
+The key feature to recognize in figure \ref{dense1} is the density
+scaling of SSD relative to other common models at any given
+temperature. SSD assumes a lower density than any of the other
 listed models at the same pressure, behavior which is especially
-apparent at temperatures greater than 300 K. Lower than expected
-densities have been observed for other systems 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}
+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 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
+SSD.\cite{Ichiye03} Throughout the remainder of the paper our
+reparamaterizations of SSD will be compared with their newer SSD1
+model.
 
-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.
-
 \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} 
 
-\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}.}
-\label{diffuse}
-\end{figure} 
+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}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{betterDiffuse.epsi} 
+%\caption{Average self-diffusion constant as a function of temperature for 
+%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
+%shown, SSD has the least deviation from the experimental values. The
+%rapidly increasing diffusion constants for TIP5P and SSD correspond to
+%significant decreases in density at the higher temperatures.}
+%\label{diffuse}
+%\end{center}
+%\end{figure} 
+
 The observed values for the diffusion constant point out one of the
-strengths of the 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 SSD model. Of the three models shown, the SSD model
+has the most accurate depiction of self-diffusion in both the
+supercooled and liquid regimes.  SPC/E does a respectable job by
+reproducing values similar to experiment around 290~K; however, it
+deviates at both higher and lower temperatures, failing to predict the
+correct thermal trend. TIP5P and SSD both start off low at colder
+temperatures and tend to diffuse too rapidly at higher temperatures.
+This behavior at higher temperatures is not particularly surprising
+since the densities of both TIP5P and SSD are lower than experimental
+water densities at higher temperatures.  When calculating the
+diffusion coefficients for SSD at experimental densities
+(instead of the densities from the NPT simulations), the resulting
+values fall more in line with experiment at these temperatures.
 
 \subsection{Structural Changes and Characterization}
+
 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}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{corrDiag.eps} 
+%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
+%\label{corrAngle}
+%\end{center}
+%\end{figure}
 
-\begin{figure}
-\includegraphics[width=45mm]{corrDiag.eps} 
-\caption{Two dimensional illustration of the angles involved in the 
-correlations observed in figure \ref{contour}.}
-\label{corrAngle}
-\end{figure}
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{fullContours.eps} 
+%\caption{Contour plots of 2D angular pair correlation functions for
+%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
+%signify regions of enhanced density while light areas signify
+%depletion relative to the bulk density. White areas have pair
+%correlation values below 0.5 and black areas have values above 1.5.}
+%\label{contour} 
+%\end{center}
+%\end{figure} 
 
-\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
+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).
+this split character alters to show the leading 4~\AA\ peak dominated
+by equatorial anti-parallel dipole orientations. There is also a
+tightly bunched group of axially arranged dipoles that most likely
+consist of the smaller fraction of aligned dipole pairs. The trailing
+component of the split peak at 5~\AA\ is dominated by aligned dipoles
+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 SSD
+density.  A better correction would be to include the
+quadrupole-quadrupole interactions for the water particles outside of
+the first solvation shell, but this would remove the simplicity and
+speed advantage of SSD.
 
-\subsection{Adjusted Potentials: SSD/E and SSD/RF}
+\subsection{Adjusted Potentials: SSD/RF and SSD/E}
+
 The propensity of 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 (SSD/RF - for use with a reaction
+field) and Soft Sticky Dipole Extended (SSD/E - an attempt to improve
+the liquid structure in simulations without a long-range correction).
 
-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\ \ \  & \ \ \ SSD [Ref. \onlinecite{Ichiye96}] \ \ \
+& \ SSD1 [Ref. \onlinecite{Ichiye03}]\ \  & \ SSD/E\ \  & \ \ 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
-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).}
-\label{grcompare} 
-\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.}
-\label{isosurface} 
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=5in
+%\epsfbox{GofRCompare.epsi} 
+%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with 
+%SSD/E and SSD1 without reaction field (top), as well as
+%SSD/RF and SSD1 with reaction field turned on
+%(bottom). The insets show the respective first peaks in detail. Note
+%how the changes in parameters have lowered and broadened the first
+%peak of SSD/E and SSD/RF.}
+%\label{grcompare} 
+%\end{center}
+%\end{figure} 
 
-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.
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{dualsticky_bw.eps} 
+%\caption{Positive and negative isosurfaces of the sticky potential for
+%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
+%correspond to the tetrahedral attractive component, and darker areas
+%correspond to the dipolar repulsive component.}
+%\label{isosurface} 
+%\end{center}
+%\end{figure} 
 
-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.
+In the original paper detailing the development of SSD, Liu and Ichiye
+placed particular emphasis on an accurate description of the first
+solvation shell. This resulted in a somewhat tall and narrow first
+peak in $g(r)$ that integrated to give similar coordination numbers to
+the experimental data obtained by Soper and
+Phillips.\cite{Ichiye96,Soper86} New experimental x-ray scattering
+data from the Head-Gordon lab indicates a slightly lower and shifted
+first peak in the g$_\mathrm{OO}(r)$, so our adjustments to 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 SSD model (SSD1), SSD/E, and SSD/RF to the new
+experimental results. Both modified water models have shorter peaks
+that match more closely to the experimental peak (as seen in the
+insets of figure \ref{grcompare}).  This structural alteration was
+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 SSD and SSD1 exclude preferred dipole alignments before the
+particles feel the pull of the ``hydrogen bonds''. Aside from
+improving the shape of the first peak in the g(\emph{r}), this
+modification improves the densities considerably by allowing the
+persistence of full dipolar character below the previous 4.0~\AA\
+cutoff.
 
-In order to demonstrate the benefits of this reparameterization, a
+While adjusting the location and shape of the first peak of $g(r)$
+improves the densities, these changes alone are insufficient to bring
+the system densities up to the values observed experimentally.  To
+further increase the densities, the dipole moments were increased in
+both of our adjusted models. Since SSD is a dipole based model, the
+structure and transport are very sensitive to changes in the dipole
+moment. The original SSD simply used the dipole moment calculated from
+the TIP3P water model, which at 2.35~D is significantly greater than
+the experimental gas phase value of 1.84~D. The larger dipole moment
+is a more realistic value and improves the dielectric properties of
+the fluid. Both theoretical and experimental measurements indicate a
+liquid phase dipole moment ranging from 2.4~D to values as high as
+3.11~D, providing a substantial range of reasonable values for a
+dipole moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
+increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF,
+respectively, leads to significant changes in the density and
+transport of the water models.
+
+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.
+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.}
-\label{ssdedense} 
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdeDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/E to 
+%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+%300 K with error bars included to clarify this region of
+%interest. Note that both SSD1 and SSD/E show good agreement with
+%experiment when the long-range correction is neglected.}
+%\label{ssdedense} 
+%\end{center}
+%\end{figure} 
 
-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 SSD/E
+model in comparison to SSD1 without a reaction field, other
+common water models, and experimental results. The calculated
+densities for both SSD/E and SSD1 have increased
+significantly over the original SSD model (see
+fig. \ref{dense1}) and are in better agreement with the experimental
+values. At 298 K, the densities of SSD/E and SSD1 without
+a long-range correction are 0.996$\pm$0.001 g/cm$^3$ and
+0.999$\pm$0.001 g/cm$^3$ respectively.  These both compare well with
+the experimental value of 0.997 g/cm$^3$, and they are considerably
+better than the SSD value of 0.967$\pm$0.003 g/cm$^3$. The
+changes to the dipole moment and sticky switching functions have
+improved the structuring of the liquid (as seen in figure
+\ref{grcompare}, but they have shifted the density maximum to much
+lower temperatures. This comes about via an increase in the liquid
+disorder through the weakening of the sticky potential and
+strengthening of the dipolar character. However, this increasing
+disorder in the SSD/E model has little effect on the melting
+transition. By monitoring $C_p$ throughout these simulations, the
+melting transition for SSD/E was shown to occur at 235~K.  The
+same transition temperature observed with SSD and SSD1.
 
-\begin{figure} 
-\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.}
-\label{ssdrfdense} 
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdrfDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/RF to 
+%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+%reparameterization when utilizing a reaction field long-ranged
+%correction - SSD/RF provides significantly more accurate
+%densities than SSD1 when performing room temperature
+%simulations.}
+%\label{ssdrfdense} 
+%\end{center}
+%\end{figure} 
 
-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
+SSD/RF and SSD1 with an active reaction field is shown in figure
+\ref{ssdrfdense}.  As observed in the simulations without a reaction
+field, the densities of SSD/RF and SSD1 show a dramatic increase over
+normal SSD (see figure \ref{dense1}). At 298 K, SSD/RF has a density
+of 0.997$\pm$0.001 g/cm$^3$, directly in line with experiment and
+considerably better than the original SSD value of 0.941$\pm$0.001
+g/cm$^3$ and the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results
+further emphasize the importance of reparameterization in order to
+model the density properly under different simulation conditions.
+Again, these changes have only a minor effect on the melting point,
+which observed at 245~K for SSD/RF, is identical to SSD and only 5~K
+lower than SSD1 with a reaction field. Additionally, the difference in
+density maxima is not as extreme, with SSD/RF showing a density
+maximum at 255~K, fairly close to the density maxima of 260~K and
+265~K, shown by SSD and SSD1 respectively.
 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdeDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/E and 
+%SSD1 (both without a reaction field) along with experimental results
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+%performed at the average densities observed in the 1 atm NPT
+%simulations for the respective models. SSD/E is slightly more mobile
+%than experiment at all of the temperatures, but it is closer to
+%experiment at biologically relevant temperatures than SSD1 without a
+%long-range correction.}
+%\label{ssdediffuse} 
+%\end{center}
+%\end{figure} 
+
 The reparameterization of the 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.
+the densities, it is important that the diffusive behavior of SSD be
+maintained or improved. Figure \ref{ssdediffuse} compares the
+temperature dependence of the diffusion constant of SSD/E to SSD1
+without an active reaction field at the densities calculated from
+their respective NPT simulations at 1 atm. The diffusion constant for
+SSD/E is consistently higher than experiment, while SSD1 remains lower
+than experiment until relatively high temperatures (around 360
+K). Both models follow the shape of the experimental curve well below
+300~K but tend to diffuse too rapidly at higher temperatures, as seen
+in SSD1's crossing above 360~K.  This increasing diffusion relative to
+the experimental values is caused by the rapidly decreasing system
+density with increasing temperature.  Both SSD1 and SSD/E show this
+deviation in particle mobility, but this trend has different
+implications on the diffusive behavior of the models.  While SSD1
+shows more experimentally accurate diffusive behavior in the high
+temperature regimes, SSD/E shows more accurate behavior in the
+supercooled and biologically relevant temperature ranges.  Thus, the
+changes made to improve the liquid structure may have had an adverse
+affect on the density maximum, but they improve the transport behavior
+of SSD/E relative to SSD1 under the most commonly simulated
+conditions.
 
-\begin{figure} 
-\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.}
-\label{ssdediffuse} 
-\end{figure} 
+%\begin{figure}
+%\begin{center} 
+%\epsfxsize=6in
+%\epsfbox{ssdrfDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/RF and 
+%SSD1 (both with an active reaction field) along with
+%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+%NVE calculations were performed at the average densities observed in
+%the 1 atm NPT simulations for both of the models. SSD/RF
+%simulates the diffusion of water throughout this temperature range
+%very well. The rapidly increasing diffusion constants at high
+%temperatures for both models can be attributed to lower calculated
+%densities than those observed in experiment.}
+%\label{ssdrfdiffuse} 
+%\end{center}
+%\end{figure} 
 
-\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.}
-\label{ssdrfdiffuse} 
-\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.
+compared to SSD1 with an active reaction field. Note that SSD/RF
+tracks the experimental results quantitatively, identical within error
+throughout most of the temperature range shown and exhibiting only a
+slight increasing trend at higher temperatures. 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. SSD/RF does not suffer from this problem as much as SSD1
+because the calculated densities are closer to the experimental
+values. These results again emphasize the importance of careful
+reparameterization when using an altered long-range correction.
 
+\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. Deviations of the of the
+averages are given in parentheses.}
+\begin{tabular}{ l  c  c  c  c  c } 
+\hline \\[-3mm]
+\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ \ SSD/E \ \ \ & \ \ SSD1 (RF) \ \ 
+\ & \ \ SSD/RF \ \ \ & \ \ Expt. \\
+ \hline \\[-3mm]
+\ \ $\rho$ (g/cm$^3$) & 0.999(0.001) & 0.996(0.001) & 0.972(0.002) & 0.997(0.001) & 0.997 \\
+\ \ $C_p$ (cal/mol K) & 28.80(0.11) & 25.45(0.09) & 28.28(0.06) & 23.83(0.16) & 17.98 \\
+\ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78(0.7) & 2.51(0.18) & 2.00(0.17) & 2.32(0.06) & 2.299\cite{Mills73} \\
+\ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 &
+4.7\footnote{Calculated by integrating $g_{\text{OO}}(r)$ in
+Ref. \onlinecite{Head-Gordon00_1}} \\ 
+\ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 &
+3.5\footnote{Calculated by integrating $g_{\text{OH}}(r)$ in
+Ref. \onlinecite{Soper86}}  \\
+\ \ $\tau_1$ (ps) & 10.9(0.6) & 7.3(0.4) & 7.5(0.7) & 7.2(0.4) & 5.7\footnote{Calculated for 298 K from data in Ref. \onlinecite{Eisenberg69}} \\
+\ \ $\tau_2$ (ps) & 4.7(0.4) & 3.1(0.2) & 3.5(0.3) & 3.2(0.2) & 2.3\footnote{Calculated for 298 K from data in
+Ref. \onlinecite{Krynicki66}}
+\end{tabular}
+\label{liquidproperties}
+\end{center}
+\end{minipage}
+\end{table}
+
+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 SSD/E and SSD/RF show an
+improvement over SSD1.  This improvement is primarily due to
+extension of the first solvation shell in the new parameter sets.
+Because $n_H$ and $n_C$ are nearly identical in SSD1, it appears
+that all molecules in the first solvation shell are involved in
+hydrogen bonds.  Since $n_H$ and $n_C$ differ in the newly
+parameterized models, the orientations in the first solvation shell
+are a bit more ``fluid''.  Therefore SSD1 overstructures the
+first solvation shell and our reparameterizations have returned this
+shell to more realistic liquid-like behavior.
+
+The time constants for the orientational autocorrelation functions
+are also displayed in Table \ref{liquidproperties}. The dipolar
+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. \onlinecite{Krynicki66} at a temperature near
+34$^\circ$C.\cite{Rahman71} Because of the strong temperature
+dependence of $\tau_2$, it is necessary to recalculate it at 298~K to
+make proper comparisons. The value shown in Table
+\ref{liquidproperties} was calculated from the same NMR data in the 
+fashion described in Ref. \onlinecite{Krynicki66}. Similarly, $\tau_1$ was
+recomputed for 298~K from the data in Ref. \onlinecite{Eisenberg69}.
+Again, SSD/E and SSD/RF show improved behavior over SSD1, both with
+and without an active reaction field. Turning on the reaction field
+leads to much improved time constants for SSD1; however, these results
+also include a corresponding decrease in system density.
+Orientational relaxation times published in the original SSD dynamics
+papers are smaller than the values observed here, and this difference
+can be attributed to the use of the Ewald sum.\cite{Ichiye99}
+
 \subsection{Additional Observations}
 
-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.
+%\begin{figure}
+%\begin{center}
+%\epsfxsize=6in
+%\epsfbox{icei_bw.eps} 
+%\caption{The most stable crystal structure assumed by the SSD family
+%of water models.  We refer to this structure as Ice-{\it i} to
+%indicate its origins in computer simulation.  This image was taken of
+%the (001) face of the crystal.}
+%\label{weirdice} 
+%\end{center}
+%\end{figure}
 
+While performing a series of melting simulations on an early iteration
+of 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.
+
+Initial simulations indicated that Ice-{\it i} is the preferred ice
+structure for at least the SSD/E model. To verify this, a comparison
+was made between near ideal crystals of ice~$I_h$, ice~$I_c$, and
+Ice-{\it i} at constant pressure with SSD/E, SSD/RF, and
+SSD1. Near-ideal versions of the three types of crystals were cooled
+to 1 K, and enthalpies of formation of each were compared using all
+three water models.  Enthalpies were estimated from the
+isobaric-isothermal simulations using $H=U+P_{\text ext}V$ where
+$P_{\text ext}$ is the applied pressure.  A constant value of -60.158
+kcal / mol has been added to place our zero for the enthalpies of
+formation for these systems at the traditional state (elemental forms
+at standard temperature and pressure).  With every model in the SSD
+family, Ice-{\it i} had the lowest calculated enthalpy of formation.
+
+\begin{table}
+\begin{center}
+\caption{Enthalpies of Formation (in kcal / mol) of the three crystal
+structures (at 1 K) exhibited by the 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]
+\ \ \ SSD/E & -72.444 & -72.450 & -73.748 \\
+\ \ \ SSD/RF & -73.093 & -73.075 & -74.180 \\
+\ \ \ SSD1 & -72.654 & -72.569 & -73.575 \\
+\ \ \ SSD1 (RF) & -72.662 & -72.569 & -73.292 \\
+\end{tabular}
+\label{iceenthalpy}
+\end{center}
+\end{table}
+
+In addition to these energetic comparisons, melting simulations were
+performed with Ice-{\it i} as the initial configuration using SSD/E,
+SSD/RF, and SSD1 both with and without a reaction field. The melting
+transitions for both SSD/E and SSD1 without reaction field occurred at
+temperature in excess of 375~K.  SSD/RF and SSD1 with a reaction field
+showed more reasonable melting transitions near 325~K.  These melting
+point observations clearly show that all of the SSD-derived models
+prefer the ice-{\it i} structure.
+
+\section{Conclusions}
+
+The density maximum and temperature dependence of the self-diffusion
+constant were studied for the 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 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 SSD model was
+reparameterized for use both with and without a reaction field (SSD/RF
+and SSD/E), and comparisons were made with SSD1, Ichiye's density
+corrected version of SSD. Both models improve the liquid structure,
+densities, and diffusive properties under their respective simulation
+conditions, indicating the necessity of reparameterization when
+changing the method of calculating long-range electrostatic
+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 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}
+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.
+
+\newpage
+
+\bibliographystyle{jcp}
+\bibliography{nptSSD}
+
+\newpage
+
+\begin{list}
+  {Figure \arabic{captions}: }{\usecounter{captions}
+	\setlength{\rightmargin}{\leftmargin}}
+	
+\item Energy conservation using both quaternion-based integration and
+the {\sc dlm} method with increasing time step. The larger time step
+plots are shifted from the true energy baseline (that of $\Delta t$ =
+0.1~fs) for clarity.
+
+\item Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E
+[Ref. \onlinecite{Clancy94}], SSD without Reaction Field, SSD, and
+experiment [Ref. \onlinecite{CRC80}]. The arrows indicate the change
+in densities observed when turning off the reaction field. The the
+lower than expected densities for the SSD model were what prompted the
+original reparameterization of SSD1 [Ref. \onlinecite{Ichiye03}].
+
+\item Average self-diffusion constant as a function of temperature for 
+SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three
+water models shown, SSD has the least deviation from the experimental
+values. The rapidly increasing diffusion constants for TIP5P and SSD
+correspond to significant decreases in density at the higher
+temperatures.
+
+\item An illustration of angles involved in the correlations observed in
+Fig. \ref{contour}.
+
+\item Contour plots of 2D angular pair correlation functions for
+512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
+signify regions of enhanced density while light areas signify
+depletion relative to the bulk density. White areas have pair
+correlation values below 0.5 and black areas have values above 1.5.
+
+\item Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with 
+SSD/E and SSD1 without reaction field (top), as well as SSD/RF and
+SSD1 with reaction field turned on (bottom). The insets show the
+respective first peaks in detail. Note how the changes in parameters
+have lowered and broadened the first peak of SSD/E and SSD/RF.
+
+\item Positive and negative isosurfaces of the sticky potential for
+SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
+correspond to the tetrahedral attractive component, and darker areas
+correspond to the dipolar repulsive component.
+
+\item Comparison of densities calculated with SSD/E to 
+SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+300 K with error bars included to clarify this region of
+interest. Note that both SSD1 and SSD/E show good agreement with
+experiment when the long-range correction is neglected.
+
+\item Comparison of densities calculated with SSD/RF to 
+SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+reparameterization when utilizing a reaction field long-ranged
+correction - SSD/RF provides significantly more accurate
+densities than SSD1 when performing room temperature
+simulations.
+
+\item The diffusion constants calculated from SSD/E and 
+SSD1 (both without a reaction field) along with experimental results
+[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+performed at the average densities observed in the 1 atm NPT
+simulations for the respective models. SSD/E is slightly more mobile
+than experiment at all of the temperatures, but it is closer to
+experiment at biologically relevant temperatures than SSD1 without a
+long-range correction.
+
+\item The diffusion constants calculated from SSD/RF and 
+SSD1 (both with an active reaction field) along with
+experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+NVE calculations were performed at the average densities observed in
+the 1 atm NPT simulations for both of the models. SSD/RF
+simulates the diffusion of water throughout this temperature range
+very well. The rapidly increasing diffusion constants at high
+temperatures for both models can be attributed to lower calculated
+densities than those observed in experiment.
+
+\item The most stable crystal structure assumed by the SSD family
+of water models.  We refer to this structure as Ice-{\it i} to
+indicate its origins in computer simulation.  This image was taken of
+the (001) face of the crystal.
+\end{list}
+
+\newpage
+
 \begin{figure}
-\includegraphics[width=85mm]{povIce.ps} 
-\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
-\label{weirdice} 
+\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}
 
-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.
+\newpage
 
 \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{denseSSDnew.eps}
+%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
+% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
+% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
+% arrows indicate the change in densities observed when turning off the
+% reaction field. The the lower than expected densities for the SSD
+% model were what prompted the original reparameterization of SSD1
+% [Ref. \onlinecite{Ichiye03}].}
+\label{dense1}
+\end{center}
 \end{figure}
 
-The 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$. 
+\newpage
 
-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.
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{betterDiffuse.epsi} 
+%\caption{Average self-diffusion constant as a function of temperature for 
+%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
+%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
+%shown, SSD has the least deviation from the experimental values. The
+%rapidly increasing diffusion constants for TIP5P and SSD correspond to
+%significant decreases in density at the higher temperatures.}
+\label{diffuse}
+\end{center}
+\end{figure}
 
-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.
+\newpage
 
-\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
-simulations of biochemical systems.
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{corrDiag.eps} 
+%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
+\label{corrAngle}
+\end{center}
+\end{figure}
 
-\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.
+\newpage
 
-\bibliographystyle{jcp}
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{fullContours.eps} 
+%\caption{Contour plots of 2D angular pair correlation functions for
+%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
+%signify regions of enhanced density while light areas signify
+%depletion relative to the bulk density. White areas have pair
+%correlation values below 0.5 and black areas have values above 1.5.}
+\label{contour} 
+\end{center}
+\end{figure} 
 
-\bibliography{nptSSD} 
+\newpage
 
-%\pagebreak 
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{GofRCompare.epsi} 
+%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with 
+%SSD/E and SSD1 without reaction field (top), as well as
+%SSD/RF and SSD1 with reaction field turned on
+%(bottom). The insets show the respective first peaks in detail. Note
+%how the changes in parameters have lowered and broadened the first
+%peak of SSD/E and SSD/RF.}
+\label{grcompare} 
+\end{center}
+\end{figure} 
 
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=7in
+\epsfbox{dualsticky_bw.eps} 
+%\caption{Positive and negative isosurfaces of the sticky potential for
+%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
+%correspond to the tetrahedral attractive component, and darker areas
+%correspond to the dipolar repulsive component.}
+\label{isosurface} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdeDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/E to 
+%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
+%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
+%300 K with error bars included to clarify this region of
+%interest. Note that both SSD1 and SSD/E show good agreement with
+%experiment when the long-range correction is neglected.}
+\label{ssdedense} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdrfDense.epsi} 
+%\caption{Comparison of densities calculated with SSD/RF to 
+%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
+%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
+%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
+%reparameterization when utilizing a reaction field long-ranged
+%correction - SSD/RF provides significantly more accurate
+%densities than SSD1 when performing room temperature
+%simulations.}
+\label{ssdrfdense} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+ \begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdeDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/E and 
+%SSD1 (both without a reaction field) along with experimental results
+%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
+%performed at the average densities observed in the 1 atm NPT
+%simulations for the respective models. SSD/E is slightly more mobile
+%than experiment at all of the temperatures, but it is closer to
+%experiment at biologically relevant temperatures than SSD1 without a
+%long-range correction.}
+\label{ssdediffuse} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center} 
+\epsfxsize=6in
+\epsfbox{ssdrfDiffuse.epsi} 
+%\caption{The diffusion constants calculated from SSD/RF and 
+%SSD1 (both with an active reaction field) along with
+%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
+%NVE calculations were performed at the average densities observed in
+%the 1 atm NPT simulations for both of the models. SSD/RF
+%simulates the diffusion of water throughout this temperature range
+%very well. The rapidly increasing diffusion constants at high
+%temperatures for both models can be attributed to lower calculated
+%densities than those observed in experiment.}
+\label{ssdrfdiffuse} 
+\end{center}
+\end{figure} 
+
+\newpage
+
+\begin{figure}
+\begin{center}
+\epsfxsize=6in
+\epsfbox{icei_bw.eps} 
+%\caption{The most stable crystal structure assumed by the SSD family
+%of water models.  We refer to this structure as Ice-{\it i} to
+%indicate its origins in computer simulation.  This image was taken of
+%the (001) face of the crystal.}
+\label{weirdice} 
+\end{center}
+\end{figure} 
+
 \end{document}