trunk/ssdePaper/nptSSD.tex

%\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[ref]{overcite}
%\usepackage{curves}
%\pagestyle{plain}
%\pagenumbering{arabic}
%\oddsidemargin 0.0cm \evensidemargin 0.0cm
%\topmargin -21pt \headsep 10pt
%\textheight 9.0in \textwidth 6.5in
%\brokenpenalty=10000
%\renewcommand{\baselinestretch}{1.2}
%\renewcommand\citemid{\ } % no comma in optional reference note
\newcounter{captions}
\newcounter{figs}

\begin{document}

\title{On the structural and transport properties of the soft sticky
dipole (SSD) and related single point water models}

\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu}}

\affiliation{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}


\begin{abstract}
The density maximum and temperature dependence of the self-diffusion
constant were investigated for the soft sticky dipole (SSD) water
model and two related 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 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

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.

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}
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)\ +
u_{ij}^{sp}
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j),
\end{equation}
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]
\ ,
\end{equation}
and
\begin{equation}
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}
u_{ij}^{sp}
({\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}
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({\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({\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}
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}

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 feature of the SSD model is that it was parameterized for
use with the Ewald sum to handle long-range interactions.  This would
normally be the best way of handling long-range interactions in
systems that contain other point charges.  However, our group has
recently become interested in systems with point dipoles as mimics for
neutral, but polarized regions on molecules (e.g. the zwitterionic
head group regions of phospholipids).  If the system of interest does
not contain point charges, the Ewald sum and even particle-mesh Ewald
become computational bottlenecks.  Their respective ideal
$N^\frac{3}{2}$ and $N\log N$ calculation scaling orders for $N$
particles can become prohibitive when $N$ becomes
large.\cite{Darden99} In applying this water model in these types of
systems, it would be useful to know its properties and behavior under
the more computationally efficient reaction field (RF) technique, or
even with a simple cutoff. This study addresses these issues by
looking at the structural and transport behavior of SSD over a
variety of temperatures with the purpose of utilizing the RF
correction technique.  We then suggest modifications to the parameters
that result in more realistic bulk phase behavior.  It should be noted
that in a recent publication, some of the original investigators of
the SSD water model have suggested adjustments to the SSD
water model to address abnormal density behavior (also observed here),
calling the corrected model SSD1.\cite{Ichiye03} In what
follows, we compare our reparamaterization of SSD with both the
original SSD and SSD1 models with the goal of improving
the bulk phase behavior of an SSD-derived model in simulations
utilizing the reaction field.

\section{Methods}

Long-range dipole-dipole interactions were accounted for in this study
by using either the reaction field technique or by resorting to a
simple cubic switching function at a cutoff radius.  One of the early
applications of a reaction field was actually in Monte Carlo
simulations of liquid water.\cite{Barker73} Under this method, the
magnitude of the reaction field acting on dipole $i$ is
\begin{equation}
\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}),
\label{rfequation}
\end{equation}
where $\mathcal{R}$ is the cavity defined by the cutoff radius
($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole
moment vector of particle $j$, and $s(r_{ij})$ is a cubic switching
function.\cite{AllenTildesley} The reaction field contribution to the
total energy by particle $i$ is given by $-\frac{1}{2}{\bf
\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf
\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.
 
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} 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, 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.  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 {\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}
%\begin{center}
%\epsfxsize=6in
%\epsfbox{timeStep.epsi}
%\caption{Energy conservation using both quaternion-based integration and
%the {\sc dlm} method with increasing time step. The larger time step
%plots are shifted from the true energy baseline (that of $\Delta t$ =
%0.1~fs) for clarity.}
%\label{timestep}
%\end{center}
%\end{figure}

In figure \ref{timestep}, the resulting energy drift at various time
steps for both the {\sc dlm} and quaternion integration schemes is
compared.  All of the 1000 SSD particle simulations started with
the same configuration, and the only difference was the method used to
handle orientational motion. At time steps of 0.1 and 0.5~fs, both
methods for propagating the orientational degrees of freedom conserve
energy fairly well, with the quaternion method showing a slight energy
drift over time in the 0.5~fs time step simulation. At time steps of 1
and 2~fs, the energy conservation benefits of the {\sc dlm} method are
clearly demonstrated. Thus, while maintaining the same degree of
energy conservation, one can take considerably longer time steps,
leading to an overall reduction in computation time.

Energy drift in the simulations using {\sc dlm} integration was
unnoticeable for time steps up to 3~fs. A slight energy drift on the
order of 0.012~kcal/mol per nanosecond was observed at a time step of
4~fs, and as expected, this drift increases dramatically with
increasing time step. To insure accuracy in our microcanonical
simulations, time steps were set at 2~fs and kept at this value for
constant pressure simulations as well.

Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
were generated as starting points for all simulations. The $I_h$
crystals were formed by first arranging the centers of mass of the SSD
particles into a ``hexagonal'' ice lattice of 1024 particles. Because
of the crystal structure of $I_h$ ice, the simulation box assumed an
orthorhombic shape with an edge length ratio of approximately
1.00$\times$1.06$\times$1.23. The particles were then allowed to
orient freely about fixed positions with angular momenta randomized at
400~K for varying times. The rotational temperature was then scaled
down in stages to slowly cool the crystals to 25~K. The particles were
then allowed to translate with fixed orientations at a constant
pressure of 1 atm for 50~ps at 25~K. Finally, all constraints were
removed and the ice crystals were allowed to equilibrate for 50~ps at
25~K and a constant pressure of 1~atm.  This procedure resulted in
structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
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,
so the ratio of edge lengths remained constant throughout the
simulations.

\section{Results and discussion}

Melting studies were performed on the randomized ice crystals using
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.

\subsection{Density Behavior}

Our initial simulations focused on the original SSD water model,
and an average density versus temperature plot is shown in figure
\ref{dense1}. Note that the density maximum when using a reaction
field appears between 255 and 265~K.  There were smaller fluctuations
in the density at 260~K than at either 255 or 265~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). 

%\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 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.

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 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.

\subsection{Transport Behavior}

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 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 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}
%\begin{center}
%\epsfxsize=6in
%\epsfbox{corrDiag.eps} 
%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
%\label{corrAngle}
%\end{center}
%\end{figure}

%\begin{figure}
%\begin{center} 
%\epsfxsize=6in
%\epsfbox{fullContours.eps} 
%\caption{Contour plots of 2D angular pair correlation functions for
%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
%signify regions of enhanced density while light areas signify
%depletion relative to the bulk density. White areas have pair
%correlation values below 0.5 and black areas have values above 1.5.}
%\label{contour} 
%\end{center}
%\end{figure} 

Additional analysis of the melting process was performed using
two-dimensional structure and dipole angle correlations. Expressions
for these correlations are as follows:

\begin{equation}
g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
\end{equation}
\begin{equation}
g_{\text{AB}}(r,\cos\omega) = 
\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
\end{equation}
where $\theta$ and $\omega$ refer to the angles shown in figure
\ref{corrAngle}. By binning over both distance and the cosine of the
desired angle between the two dipoles, the $g(r)$ can be analyzed to
determine the common dipole arrangements that constitute the peaks and
troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
of figure \ref{contour} show results from an ice $I_c$ simulation. The
first peak in the $g(r)$ consists primarily of the preferred hydrogen
bonding arrangements as dictated by the tetrahedral sticky potential -
one peak for the hydrogen bond donor and the other for the hydrogen
bond acceptor.  Due to the high degree of crystallinity of the sample,
the second and third solvation shells show a repeated peak arrangement
which decays at distances around the fourth solvation shell, near the
imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
In the higher temperature simulation shown in frames C and D, 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
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 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. 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 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/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 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 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).

\begin{table}
\begin{center}
\caption{Parameters for the original and adjusted models}
\begin{tabular}{ l  c  c  c  c } 
\hline \\[-3mm]
\ \ \ Parameters\ \ \  & \ \ \ SSD [Ref. \onlinecite{Ichiye96}] \ \ \
& \ SSD1 [Ref. \onlinecite{Ichiye03}]\ \  & \ SSD/E\ \  & \ \ SSD/RF \\
\hline \\[-3mm]
\ \ \ $\sigma$ (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
\ \ \ $\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}
%\begin{center} 
%\epsfxsize=5in
%\epsfbox{GofRCompare.epsi} 
%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with 
%SSD/E and SSD1 without reaction field (top), as well as
%SSD/RF and SSD1 with reaction field turned on
%(bottom). The insets show the respective first peaks in detail. Note
%how the changes in parameters have lowered and broadened the first
%peak of SSD/E and SSD/RF.}
%\label{grcompare} 
%\end{center}
%\end{figure} 

%\begin{figure}
%\begin{center} 
%\epsfxsize=6in
%\epsfbox{dualsticky_bw.eps} 
%\caption{Positive and negative isosurfaces of the sticky potential for
%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
%correspond to the tetrahedral attractive component, and darker areas
%correspond to the dipolar repulsive component.}
%\label{isosurface} 
%\end{center}
%\end{figure} 

In the original paper detailing the development of 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.

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 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}
%\begin{center} 
%\epsfxsize=6in
%\epsfbox{ssdeDense.epsi} 
%\caption{Comparison of densities calculated with SSD/E to 
%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
%300 K with error bars included to clarify this region of
%interest. Note that both SSD1 and SSD/E show good agreement with
%experiment when the long-range correction is neglected.}
%\label{ssdedense} 
%\end{center}
%\end{figure} 

Fig. \ref{ssdedense} shows the density profile for the 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}
%\begin{center} 
%\epsfxsize=6in
%\epsfbox{ssdrfDense.epsi} 
%\caption{Comparison of densities calculated with SSD/RF to 
%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
%reparameterization when utilizing a reaction field long-ranged
%correction - SSD/RF provides significantly more accurate
%densities than SSD1 when performing room temperature
%simulations.}
%\label{ssdrfdense} 
%\end{center}
%\end{figure} 

Including the reaction field long-range correction in the simulations
results in a more interesting comparison.  A density profile including
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 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}
%\begin{center} 
%\epsfxsize=6in
%\epsfbox{ssdrfDiffuse.epsi} 
%\caption{The diffusion constants calculated from SSD/RF and 
%SSD1 (both with an active reaction field) along with
%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
%NVE calculations were performed at the average densities observed in
%the 1 atm NPT simulations for both of the models. SSD/RF
%simulates the diffusion of water throughout this temperature range
%very well. The rapidly increasing diffusion constants at high
%temperatures for both models can be attributed to lower calculated
%densities than those observed in experiment.}
%\label{ssdrfdiffuse} 
%\end{center}
%\end{figure} 

In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
compared to SSD1 with an active reaction field. Note that SSD/RF
tracks the experimental results quantitatively, identical within error
throughout most of the temperature range shown and exhibiting only a
slight increasing trend at higher temperatures. 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}

%\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}
\begin{center}
\epsfxsize=6in
\epsfbox{timeStep.epsi}
%\caption{Energy conservation using both quaternion-based integration and
%the {\sc dlm} method with increasing time step. The larger time step
%plots are shifted from the true energy baseline (that of $\Delta t$ =
%0.1~fs) for clarity.}
\label{timestep}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{denseSSDnew.eps}
%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
% arrows indicate the change in densities observed when turning off the
% reaction field. The the lower than expected densities for the SSD
% model were what prompted the original reparameterization of SSD1
% [Ref. \onlinecite{Ichiye03}].}
\label{dense1}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{betterDiffuse.epsi} 
%\caption{Average self-diffusion constant as a function of temperature for 
%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
%shown, SSD has the least deviation from the experimental values. The
%rapidly increasing diffusion constants for TIP5P and SSD correspond to
%significant decreases in density at the higher temperatures.}
\label{diffuse}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{corrDiag.eps} 
%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
\label{corrAngle}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{fullContours.eps} 
%\caption{Contour plots of 2D angular pair correlation functions for
%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
%signify regions of enhanced density while light areas signify
%depletion relative to the bulk density. White areas have pair
%correlation values below 0.5 and black areas have values above 1.5.}
\label{contour} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{GofRCompare.epsi} 
%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with 
%SSD/E and SSD1 without reaction field (top), as well as
%SSD/RF and SSD1 with reaction field turned on
%(bottom). The insets show the respective first peaks in detail. Note
%how the changes in parameters have lowered and broadened the first
%peak of SSD/E and SSD/RF.}
\label{grcompare} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=7in
\epsfbox{dualsticky_bw.eps} 
%\caption{Positive and negative isosurfaces of the sticky potential for
%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
%correspond to the tetrahedral attractive component, and darker areas
%correspond to the dipolar repulsive component.}
\label{isosurface} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdeDense.epsi} 
%\caption{Comparison of densities calculated with SSD/E to 
%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
%300 K with error bars included to clarify this region of
%interest. Note that both SSD1 and SSD/E show good agreement with
%experiment when the long-range correction is neglected.}
\label{ssdedense} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdrfDense.epsi} 
%\caption{Comparison of densities calculated with SSD/RF to 
%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
%reparameterization when utilizing a reaction field long-ranged
%correction - SSD/RF provides significantly more accurate
%densities than SSD1 when performing room temperature
%simulations.}
\label{ssdrfdense} 
\end{center}
\end{figure} 

\newpage

 \begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdeDiffuse.epsi} 
%\caption{The diffusion constants calculated from SSD/E and 
%SSD1 (both without a reaction field) along with experimental results
%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
%performed at the average densities observed in the 1 atm NPT
%simulations for the respective models. SSD/E is slightly more mobile
%than experiment at all of the temperatures, but it is closer to
%experiment at biologically relevant temperatures than SSD1 without a
%long-range correction.}
\label{ssdediffuse} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdrfDiffuse.epsi} 
%\caption{The diffusion constants calculated from SSD/RF and 
%SSD1 (both with an active reaction field) along with
%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
%NVE calculations were performed at the average densities observed in
%the 1 atm NPT simulations for both of the models. SSD/RF
%simulates the diffusion of water throughout this temperature range
%very well. The rapidly increasing diffusion constants at high
%temperatures for both models can be attributed to lower calculated
%densities than those observed in experiment.}
\label{ssdrfdiffuse} 
\end{center}
\end{figure} 

\newpage

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{icei_bw.eps} 
%\caption{The most stable crystal structure assumed by the SSD family
%of water models.  We refer to this structure as Ice-{\it i} to
%indicate its origins in computer simulation.  This image was taken of
%the (001) face of the crystal.}
\label{weirdice} 
\end{center}
\end{figure} 

\end{document}
Revision:	1033
Committed:	Fri Feb 6 19:15:44 2004 UTC (20 years, 4 months ago) by chrisfen
Content type:	application/x-tex
File size:	66536 byte(s)
Log Message:	I hate revtex