trunk/ssdePaper/nptSSD.tex

%\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
\documentclass[11pt]{article}
\usepackage{endfloat}
\usepackage{amsmath}
\usepackage{epsf}
\usepackage{berkeley}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage[ref]{overcite}
%\usepackage{berkeley}
%\usepackage{curves}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

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

\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry\\ University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle

\begin{abstract}
NVE and NPT molecular dynamics simulations were performed in order to
investigate the density maximum and temperature dependent transport
for SSD and related water models, 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 water 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. Compared to the density
corrected version of SSD (SSD1), these modified models were shown to
maintain or improve upon the structural and transport properties.
\end{abstract}

\newpage

%\narrowtext 


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

\section{Introduction}

One of the most important tasks in the simulation 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 clusters or
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 biomolecular 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 predict specific properties more accurately than their
predecessors, but often at the cost of other properties or of computer
time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P improves
upon 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 succeeds in both retaining the accuracy of system
properties and simplifying calculations to increase computational
efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96}

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
\begin{equation}
u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
u_{ij}^{sp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\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:
\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}
\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}\ ,
\end{equation}
\begin{equation}
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) + s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
\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, and $s$
and $s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
functions take the following forms:
\begin{equation}
w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
\end{equation}
\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,
\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}

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 an
increase in calculation efficiency 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, it
was shown that SSD improves on the prediction of many of water's
dynamical properties over TIP3P and SPC/E.\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.

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}
In applying this water model in these types of systems, it would be
useful to know its properties and behavior with the more
computationally efficient reaction field (RF) technique, and even with
a cutoff that lacks any form of long-range correction. 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 alterations to
the parameters that result in more water-like behavior. It should be
noted that in a recent publication, some of the original investigators of
the SSD water model have put forth adjustments to the SSD water model
to address abnormal density behavior (also observed here), calling the
corrected model SSD1.\cite{Ichiye03} This study will make comparisons
with SSD1's behavior with the goal of improving upon the
depiction of water under conditions without the Ewald Sum.

\section{Methods}

As stated previously, the long-range dipole-dipole interactions were
accounted for in this study by using the reaction field method. The
magnitude of the reaction field acting on dipole \emph{i} is given by
\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})\  ,
\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
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 systems, 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 presented on how to make
SSD more accurate both with and without a reaction field.

Simulations were performed in both the isobaric-isothermal and
microcanonical ensembles. The constant pressure simulations were
implemented using an integral thermostat and barostat as outlined by
Hoover.\cite{Hoover85,Hoover86} All particles 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 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}

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

The symplectic splitting method allows for Verlet style integration of
both linear and angular 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 computationally 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 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}.

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{timeStep.epsi}
\caption{Energy conservation using quaternion based integration versus 
the symplectic step method proposed by Dullweber \emph{et al.} with
increasing time step. The larger time step plots are shifted up 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 used to
handle 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. 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 symplectic step simulations was unnoticeable for
time steps up to three femtoseconds. A slight energy drift on the
order of 0.012 kcal/mol per nanosecond was observed at a time step of
four femtoseconds, and as expected, this drift increases dramatically
with increasing time step. To insure accuracy in the constant energy
simulations, time steps were set at 2 fs and kept at this value for
constant pressure simulations as well.

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 a rectangular 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
constant pressure and temperature 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 the system volume
fluctuations dampened out in all but the very cold simulations (below
225 K).

\subsection{Density Behavior}
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, where the calculated densities
within this range were nearly indistinguishable. 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; 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 to deformation into a dense glassy state at lower
temperatures. This resulted in an overall low temperature density
maximum at 200 K, while still retaining a liquid state density maximum
in common with the $I_h$ simulations.

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{denseSSD.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} 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.\cite{Ichiye03}}
\label{dense1}
\end{center}
\end{figure}

The density maximum for SSD actually 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 included in
this particular plot, it is useful to note that TIP5P has a water
density maximum nearly identical to experiment.

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
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
significant increases 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. Note that the SSD model assumes a lower density than any
of the other listed models at the same pressure, behavior which is
especially apparent at temperatures greater than 300 K. Lower than
expected densities have been observed for other systems 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
reaction field, the 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 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 a
reaction field, the density around 300 K is still significantly lower
than experiment and comparable water models. This anomalous behavior
was what lead Ichiye \emph{et al.} to recently reparameterize SSD and
make SSD1.\cite{Ichiye03} In discussing potential adjustments later in
this paper, all comparisons were performed with this new model.

\subsection{Transport Behavior}
Of importance in these types of studies are the transport properties
of the particles and their change in responce to altering
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 five sets of NVE
simulations are displayed in figure \ref{diffuse}, alongside
experimental, SPC/E, and TIP5P
results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01}

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{betterDiffuse.epsi} 
\caption{Average diffusion coefficient over increasing temperature for 
SSD, SPC/E,\cite{Clancy94} TIP5P,\cite{Jorgensen01} and Experimental
data.\cite{Gillen72,Mills73} 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 decrease 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 liquid temperature
regimes. SPC/E does a respectable job by producing values similar to
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 colder temperatures and
tend to diffuse too rapidly at higher temperatures. This trend at
higher temperatures is not surprising in that the densities of both
TIP5P and SSD are lower than experimental water at these higher
temperatures. 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.

\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. 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, but this is not a
surprise considering the simplicity of the SSD model.

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

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{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{center}
\end{figure} 

Additional analysis of the melting phase-transition process was
performed by 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|\mathbf{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|\mathbf{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(\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}) consists primarily 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
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
longer-ranged repeated peak 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.

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 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. Primary energetically favorable dipole
arrangements populate the region immediately outside this repulsion
region (around 4 \AA), while arrangements that seek to ideally 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 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,
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 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.

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

The parameters available 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}
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)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)],
\end{equation}

where $\nu_0$ scales the strength of the tetrahedral attraction and
$\nu_0^\prime$ acts in an identical fashion on the dipole repulsion
term. The separation was performed for purposes of the
reparameterization, 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 the tetrahedral attractive and dipolar
repulsive terms do not share the same lower cutoff ($r_l$) in the
newly parameterized potentials - soft sticky dipole reaction field
(SSD/RF - for use with a reaction field) and soft sticky dipole
enhanced (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\cite{Ichiye96} \ \ \ & \ SSD1\cite{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\\
\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
\ \ \ $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\cite{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.ps} 
\caption{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 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 the
g(\emph{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 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 revised SSD model (SSD1), SSD-E, and SSD-RF to the new
experimental results. Both 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 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(\emph{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 the 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 earlier in this paper.

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdeDense.epsi} 
\caption{Comparison of densities calculated with SSD/E to SSD1 without a 
reaction field, TIP3P,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00}
SPC/E,\cite{Clancy94} and experiment.\cite{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 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 figure \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$\text{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,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00}
SPC/E,\cite{Clancy94} and experiment.\cite{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 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{Plots of the diffusion constants calculated from SSD/E and SSD1,
 both without a reaction field, along with experimental
 results.\cite{Gillen72,Mills73} 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 fluid than experiment
 at all of the temperatures, but it is closer 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 SSD1 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. The diffusion constant for SSD/E is consistently a little
higher than experiment, while SSD1 remains lower than experiment until
relatively high temperatures (greater than 330 K). Both models follow
the shape of the experimental curve well below 300 K but tend to
diffuse too rapidly at higher temperatures, something that is
especially apparent with SSD1. This accelerated increasing of
diffusion is caused by the rapidly decreasing system density with
increasing temperature. Though it is difficult to see in figure
\ref{ssdedense}, the densities of SSD1 decay more rapidly with
temperature than do those of SSD/E, leading to more visible deviation
from the experimental diffusion trend. 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.

\begin{figure}
\begin{center} 
\epsfxsize=6in
\epsfbox{ssdrfDiffuse.epsi} 
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1,
 both with an active reaction field, along with experimental
 results.\cite{Gillen72,Mills73} The NVE calculations were performed
 at the average densities observed in the 1 atm NPT simulations for
 both of the models. Note how accurately SSD/RF simulates the
 diffusion of water throughout this temperature range. The more
 rapidly increasing diffusion constants at high temperatures for both
 models is attributed to the significantly lower densities than
 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 incredibly well, identical within
error throughout the temperature range shown and with 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 in relation to SSD/E, 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 value. These results again emphasize the
importance of careful reparameterization when using an altered
long-range correction.

\subsection{Additional Observations}

\begin{figure}
\begin{center}
\epsfxsize=6in
\epsfbox{povIce.ps} 
\caption{A water lattice built from the crystal structure assumed by 
 SSD/E when undergoing an extremely restricted temperature NPT
 simulation. This form of ice is referred to as ice \emph{i} to
 emphasize its simulation origins. This image was taken of the (001)
 face of the crystal.}
\label{weirdice} 
\end{center}
\end{figure}

While performing restricted temperature melting sequences of SSD/E not
previously discussed, some interesting observations were made. After
melting at 235 K, two of five 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. The final
configurations of these two melting sequences show 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 $\sqrt{\smash[b]{-\text{I}}}$ (ice
\emph{i}). The crystallinity was extensive enough that a near ideal
crystal structure of ice \emph{i} could be obtained. 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 flexed rather than perfectly
straight. This results in a skewed tetrahedral geometry about the
central molecule. Referring to figure \ref{isosurface}, these flexed
hydrogen bonds are allowed due to the conical shape of the attractive
regions, with the greatest attraction along the direct hydrogen bond
configuration. Though not ideal, these flexed hydrogen bonds are
favorable enough to stabilize an entire crystal generated around
them. In fact, the imperfect ice \emph{i} crystals were so stable that
they melted at temperatures nearly 100 K greater than both ice I$_c$
and I$_h$.

These initial simulations indicated that ice \emph{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 0 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 the potential energies of each were compared using all
three water models. With every water model, ice \emph{i} turned out to
have the lowest potential energy: 5\% lower than $I_h$ with SSD1,
6.5\% lower with SSD/E, and 7.5\% lower with SSD/RF.

In addition to these low temperature comparisons, melting sequences
were performed with ice \emph{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 a 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 emphasize the preference for this
crystal structure over the most common types of ice when using these
single point water models.

Recognizing that the above tests show ice \emph{i} to be both the most
stable and lowest density crystal structure for these single point
water models, it is interesting to speculate on the relative stability
of this crystal structure with charge based water 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 the crystals to compare the system
energies. Again, ice \emph{i} was observed to have the lowest total
system energy. The total energy of ice \emph{i} was ~2\% lower than
ice $I_h$, which was in turn ~3\% lower than ice $I_c$. Based on these
initial studies, it would not be surprising if results from the other
common water models show ice \emph{i} to be the lowest energy crystal
structure. A continuation of this work studying ice \emph{i} with
multi-point water models will be published in a coming article.

\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 water well in both
the liquid and super-cooled 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
comparison simulations were performed with SSD1, the density corrected
version of SSD. Both models improve the liquid structure, density
values, and diffusive properties under their respective conditions,
indicating the necessity of reparameterization when altering the
long-range correction specifics. When taking into account the
appropriate considerations, these simple water models are excellent
choices for representing explicit water in large scale simulations of
biochemical systems.

\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 00 79647.


\newpage

\bibliographystyle{jcp}
\bibliography{nptSSD} 

%\pagebreak 

\end{document}
Revision:	863
Committed:	Thu Nov 13 15:55:20 2003 UTC (21 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	49961 byte(s)
Log Message:	Changed ice 0 to ice i.