trunk/ssdePaper/nptSSD.tex

\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
%\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
\usepackage{amsmath}
\usepackage{graphicx}

%\usepackage{endfloat}
%\usepackage{berkeley}
%\usepackage{epsf}
%\usepackage[ref]{overcite}
%\usepackage{setspace}
%\usepackage{tabularx}
%\usepackage{graphicx}
%\usepackage{curves}
%\usepackage{amsmath}
%\pagestyle{plain}
%\pagenumbering{arabic}
%\oddsidemargin 0.0cm \evensidemargin 0.0cm
%\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{\thefootnote}
\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}

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

\date{\today}

\begin{abstract}
NVE and NPT molecular dynamics simulations were performed in order to
investigate the density maximum and temperature dependent transport
for the SSD water model, both with and without the use of reaction
field. The constant pressure simulations of the melting of both $I_h$
and $I_c$ ice showed a density maximum near 260 K. In most cases, the
calculated densities were significantly lower than the densities
calculated in simulations of other water models. Analysis of particle
diffusion showed SSD to capture the transport properties of
experimental very well in both the normal and super-cooled liquid
regimes. In order to correct the density behavior, SSD was
reparameterized for use both with and without a long-range interaction
correction, SSD/RF and SSD/E respectively. In addition to correcting
the abnormally low densities, these new versions were show to maintain
or improve upon the transport and structural features of the original
water model.
\end{abstract}

\maketitle

%\narrowtext 


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

\section{Introduction}

One of the most important tasks in simulations of biochemical systems
is the proper depiction of water and water solvation. In fact, the
bulk of the calculations performed in solvated simulations are of
interactions with or between solvent molecules. Thus, the outcomes of
these types of simulations are highly dependent on the physical
properties of water, both as individual molecules and in
groups/bulk. Due to the fact that explicit solvent accounts for a
massive portion of the calculations, it necessary to simplify the
solvent to some extent in order to complete simulations in a
reasonable amount of time. In the case of simulating water in
bio-molecular studies, the balance between accurate properties and
computational efficiency is especially delicate, and it has resulted
in a variety of different water
models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these models
get specific properties correct or better than their predecessors, but
this is often at a cost of some other properties or of computer
time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P succeeds
in improving the structural and transport properties over its
predecessors, yet this comes at a greater than 50\% increase in
computational cost.\cite{Jorgensen01,Jorgensen00} One recently
developed model that succeeds in both retaining accuracy of system
properties and simplifying calculations to increase computational
efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96}

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}
\begin{split}
u_{ij}^{sp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) 
&=
\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
& \quad \ +
s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
\end{split}
\end{equation}
where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
$\nu_0$ scales the strength of the overall sticky potential, $s$ and
$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
functions take the following forms,
\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 a
calculation speed up of up to an order of magnitude over other
comparable models while maintaining the structural behavior of
water.\cite{Ichiye96} In the original molecular dynamics studies, 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, a specific
interest within our group.

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. Towards the end, we suggest
alterations to the parameters that result in more water-like
behavior. It should be noted that in a recent publication, some the
original investigators of the SSD water model have put forth
adjustments to the original SSD water model to address abnormal
density behavior (also observed here), calling the corrected model
SSD1.\cite{Ichiye03} This study will consider this new model's
behavior as well, and hopefully improve upon its depiction of water
under conditions without the Ewald Sum.

\section{Methods}

As stated previously, in this study the long-range dipole-dipole
interactions were accounted for using the reaction field method. The
magnitude of the reaction field acting on dipole \emph{i} is given by
\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 system, the computational cost benefit of reaction field
is dramatic. To address some of the dynamical property alterations due
to the use of reaction field, simulations were also performed without
a surrounding dielectric and suggestions are proposed on how to make
SSD more compatible with a reaction field.

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 a option, being that the rotation matrix for a single body is
nine elements long as opposed to 3 or 4 elements for Euler angles and
quaternions respectively. System memory has become much less of an
issue in recent times, and this has resulted in substantial benefits
in energy conservation. There is still the issue of 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 the integration
method, the orientational propagation involves a sequence of matrix
evaluations to update the rotation matrix.\cite{Dullweber1997} These
matrix rotations end up being more costly computationally than the
simpler arithmetic quaternion propagation. 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}
\includegraphics[width=61mm, angle=-90]{timeStep.epsi}
\caption{Energy conservation using quaternion based integration versus 
the symplectic step method proposed by Dullweber \emph{et al.} with
increasing time step. For each time step, the dotted line is total
energy using the symplectic step integrator, and the solid line comes
from the quaternion integrator. The larger time step plots are shifted
up from the true energy baseline for clarity.}
\label{timestep}
\end{figure}

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

Energy drift in these SSD particle simulations was unnoticeable for
time steps up to three femtoseconds. A slight energy drift on the
order of 0.012 kcal/mol per nanosecond was observed at a time step of
four femtoseconds, and as expected, this drift increases dramatically
with increasing time step. To insure accuracy in the constant energy
simulations, time steps were set at 2 fs and kept at this value for
constant pressure simulations as well.

Ice crystals in both the $I_h$ and $I_c$ lattices were generated as
starting points for all the simulations. The $I_h$ crystals were
formed by first arranging the center of masses of the SSD particles
into a ``hexagonal'' ice lattice of 1024 particles. Because of the
crystal structure of $I_h$ ice, the simulation box assumed a
rectangular shape with a edge length ratio of approximately
1.00$\times$1.06$\times$1.23. The particles were then allowed to
orient freely about fixed positions with angular momenta randomized at
400 K for varying times. The rotational temperature was then scaled
down in stages to slowly cool the crystals down to 25 K. The particles
were then allowed translate with fixed orientations at a constant
pressure of 1 atm for 50 ps at 25 K. Finally, all constraints were
removed and the ice crystals were allowed to equilibrate for 50 ps at
25 K and a constant pressure of 1 atm.  This procedure resulted in
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. This involved an initial
randomization of velocities about the starting temperature of 25 K for
varying amounts of time. The systems were all equilibrated for 100 ps
prior to a 200 ps data collection run at each temperature setting,
ranging from 25 to 400 K, with a maximum degree increment of 25 K. For
regions of interest along this stepwise progression, the temperature
increment was decreased from 25 K to 10 and then 5 K. The above
equilibration and production times were sufficient in that the system
volume fluctuations dampened out in all but the very cold simulations
(below 225 K). In order to further improve statistics, an ensemble
average was accumulated from five separate simulation progressions,
each starting from a different ice crystal.

\subsection{Density Behavior}
In the initial average density versus temperature plot, the density
maximum clearly appears between 255 and 265 K. The calculated
densities within this range were nearly indistinguishable, as can be
seen in the zoom of this region of interest, shown in figure
\ref{dense1}. The greater certainty of the average value at 260 K makes 
a good argument for the actual density maximum residing at this
midpoint value. Figure \ref{dense1} was constructed using ice $I_h$
crystals for the initial configuration; and though not pictured, the
simulations starting from ice $I_c$ crystal configurations showed
similar results, with a liquid-phase density maximum in this same
region (between 255 and 260 K). In addition, the $I_c$ crystals are
more fragile than the $I_h$ crystals, leading them to deform into a
dense glassy state at lower temperatures. This resulted in an overall
low temperature density maximum at 200 K, but they still retained a
common liquid state density maximum with the $I_h$ simulations.

\begin{figure}
\includegraphics[width=65mm,angle=-90]{dense2.eps}
\caption{Density versus temperature for TIP4P\cite{Jorgensen98b}, 
TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction
Field, SSD, and Experiment\cite{CRC80}. }
\label{dense2}
\end{figure}

The density maximum for SSD actually compares quite favorably to other
simple water models. Figure \ref{dense2} shows a plot of these
findings with the density progression of several other models and
experiment obtained from other
sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
models, SSD has results closest to the experimentally observed water
density maximum. Of the listed water models, TIP4P has a density
maximum behavior most like that seen in SSD. Though not shown, it is
useful to note that TIP5P has a water density maximum nearly identical
to experiment.

Possibly of more importance is the density scaling of SSD relative to
other common models at any given temperature (Fig. \ref{dense2}). Note
that the SSD model assumes a lower density than any of the other
listed models at the same pressure, behavior which is especially
apparent at temperatures greater than 300 K. Lower than expected
densities have been observed for other systems with the use of a
reaction field for long-range electrostatic interactions, so the most
likely reason for these significantly lower densities in these
simulations is the presence of the reaction field.\cite{Berendsen98}
In order to test the effect of the reaction field on the density of
the systems, the simulations were repeated for the temperature region
of interest without a reaction field present. The results of these
simulations are also displayed in figure \ref{dense2}. Without
reaction field, these densities increase considerably to more
experimentally reasonable values, especially around the freezing point
of liquid water. The shape of the curve is similar to the curve
produced from SSD simulations using reaction field, specifically the
rapidly decreasing densities at higher temperatures; however, a slight
shift in the density maximum location, down to 245 K, is
observed. This is probably a more accurate comparison to the other
listed water models in that no long range corrections were applied in
those simulations.\cite{Clancy94,Jorgensen98b}

It has been observed that densities are dependent on the cutoff radius
used for a variety of water models in simulations both with and
without the use of reaction field.\cite{Berendsen98} In order to
address the possible affect of cutoff radius, simulations were
performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the
previous SSD simulations, all performed with a cutoff of 9.0 \AA. All
the resulting densities overlapped within error and showed no
significant trend in lower or higher densities as a function of cutoff
radius, both for simulations with and without reaction field. These
results indicate that there is no major benefit in choosing a longer
cutoff radius in simulations using SSD. This is comforting in that the
use of a longer cutoff radius results in a near doubling of the time
required to compute a single trajectory.

\subsection{Transport Behavior}
Of importance in these types of studies are the transport properties
of the particles and how they change when altering the environmental
conditions. In order to probe transport, constant energy simulations
were performed about the average density uncovered by the constant
pressure simulations. Simulations started with randomized velocities
and underwent 50 ps of temperature scaling and 50 ps of constant
energy equilibration before obtaining a 200 ps trajectory. Diffusion
constants were calculated via root-mean square deviation analysis. The
averaged results from 5 sets of these NVE simulations is displayed in
figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} 

\begin{figure}
\includegraphics[width=65mm, angle=-90]{betterDiffuse.epsi} 
\caption{Average diffusion coefficient over increasing temperature for 
SSD, SPC/E\cite{Clancy94}, TIP5P\cite{Jorgensen01}, and Experimental
data from Gillen \emph{et al.}\cite{Gillen72}, and from
Mills\cite{Mills73}.}
\label{diffuse}
\end{figure} 

The observed values for the diffusion constant point out one of the
strengths of the SSD model. Of the three experimental models shown,
the SSD model has the most accurate depiction of the diffusion trend
seen in experiment in both the supercooled and normal regimes. SPC/E
does a respectable job by getting similar values as SSD and experiment
around 290 K; however, it deviates at both higher and lower
temperatures, failing to predict the experimental trend. TIP5P and SSD
both start off low at the colder temperatures and tend to diffuse too
rapidly at the higher temperatures. This type of trend at the higher
temperatures is not surprising in that the densities of both TIP5P and
SSD are lower than experimental water at temperatures higher than room
temperature. When calculating the diffusion coefficients for SSD at
experimental densities, the resulting values fall more in line with
experiment at these temperatures, albeit not at standard
pressure. Results under these conditions can be found later in this
paper.

\subsection{Structural Changes and Characterization}
By starting the simulations from the crystalline state, the melting
transition and the ice structure can be studied along with the liquid
phase behavior beyond the melting point. To locate the melting
transition, the constant pressure heat capacity (C$_\text{p}$) was
monitored in each of the simulations. In the melting simulations of
the 1024 particle ice $I_h$ simulations, a large spike in C$_\text{p}$
occurs at 245 K, indicating a first order phase transition for the
melting of these ice crystals. When the reaction field is turned off,
the melting transition occurs at 235 K.  These melting transitions are
considerably lower than the experimental value, but this is not
surprising in that SSD is a simple rigid body model with a fixed
dipole.

\begin{figure} 
\includegraphics[width=85mm]{fullContours.eps} 
\caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at 
100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for
clarity: dark areas signify peaks while light areas signify
depressions. White areas have g(\emph{r}) values below 0.5 and black
areas have values above 1.5.}
\label{contour} 
\end{figure} 

Additional analyses for understanding the melting phase-transition
process were performed via two-dimensional structure and dipole angle
correlations. Expressions for the correlations are as follows:

\begin{figure}
\includegraphics[width=45mm]{corrDiag.eps} 
\caption{Two dimensional illustration of the angles involved in the 
correlations observed in figure \ref{contour}.}
\label{corrAngle}
\end{figure}

\begin{multline}
g_{\text{AB}}(r,\cos\theta) = \\
\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
\end{multline}
\begin{multline}
g_{\text{AB}}(r,\cos\omega) = \\
\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
\end{multline}
where $\theta$ and $\omega$ refer to the angles shown in the above
illustration. By binning over both distance and the cosine of the
desired angle between the two dipoles, the g(\emph{r}) can be
dissected to determine the common dipole arrangements that constitute
the peaks and troughs. Frames A and B of figure \ref{contour} show a
relatively crystalline state of an ice $I_c$ simulation. The first
peak of the g(\emph{r}) primarily consists of the preferred hydrogen
bonding arrangements as dictated by the tetrahedral sticky potential,
one peak for the donating and the other for the accepting hydrogen
bonds. Due to the high degree of crystallinity of the sample, the
second and third solvation shells show a repeated peak arrangement
which decays at distances around the fourth solvation shell, near the
imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
In the higher temperature simulation shown in frames C and D, the
repeated peak features are significantly blurred. The first solvation
shell still shows the strong effect of the sticky-potential, although
it covers a larger area, extending to include a fraction of aligned
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 parts that
blend together to form one observable peak. At higher temperatures,
this split character alters to show the leading 4 \AA\ peak dominated
by equatorial anti-parallel dipole orientations, and there is tightly
bunched group of axially arranged dipoles that most likely consist of
the smaller fraction aligned dipole pairs. The trailing part of the
split peak at 5 \AA\ is dominated by aligned dipoles that range
primarily within the axial to the chief hydrogen bond arrangements
similar to those seen in the first solvation shell. This evidence
indicates that the dipole pair interaction begins to dominate outside
of the range of the dipolar repulsion term, with the primary
energetically favorable dipole arrangements populating the region
immediately outside of it's range (around 4 \AA), and arrangements
that seek to ideally satisfy both the sticky and dipole forces locate
themselves just beyond this region (around 5 \AA).

From these findings, the split second peak is primarily the product of
the dipolar repulsion term of the sticky potential. In fact, the
leading of the two peaks can be pushed out and merged with the outer
split peak just by extending the switching function cutoff
($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even
5 \AA. This type of correction is not recommended for improving the
liquid structure, because the second solvation shell will still be
shifted too far out. In addition, this would have an even more
detrimental effect on the system densities, leading to a liquid with a
more open structure and a density considerably lower than the normal
SSD behavior shown previously. A better correction would be to include
the quadrupole-quadrupole interactions for the water particles outside
of the first solvation shell, but this reduces the simplicity and
speed advantage of SSD, so it is not the most desirable path to take.

\subsection{Adjusted Potentials: SSD/E and SSD/RF}
The propensity of SSD to adopt lower than expected densities under
varying conditions is troubling, especially at higher temperatures. In
order to correct this behavior, it's necessary to adjust the force
field parameters for the primary intermolecular interactions. In
undergoing a reparameterization, it is important not to focus on just
one property and neglect the other important properties. In this case,
it would be ideal to correct the densities while maintaining the
accurate transport properties.

The possible parameters for tuning include the $\sigma$ and $\epsilon$
Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
attractive and dipole repulsive terms with their respective
cutoffs. To alter the attractive and repulsive terms of the sticky
potential independently, it is necessary to separate the terms as
follows:
\begin{equation}
\begin{split}
u_{ij}^{sp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) &=
\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\\
& \quad \ + \frac{\nu_0^\prime}{2}
[s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)],
\end{split}
\end{equation}

where $\nu_0$ scales the strength of the tetrahedral attraction and
$\nu_0^\prime$ acts in an identical fashion on the dipole repulsion
term. For purposes of the reparameterization, the separation was
performed, but the final parameters were adjusted so that it is
unnecessary to separate the terms when implementing the adjusted water
potentials. The results of the reparameterizations are shown in table
\ref{params}. Note that both the tetrahedral attractive and dipolar 
repulsive don't share the same lower cutoff ($r_l$) in the newly
parameterized potentials - soft sticky dipole enhanced (SSD/E) and
soft sticky dipole reaction field (SSD/RF).

\begin{table}
\caption{Parameters for the original and adjusted models}
\begin{tabular}{ l  c  c  c } 
\hline \\[-3mm]
\ Parameters & \ \ \  SSD$^\dagger$\ \ \ \  & \ SSD/E\ \  & \ SSD/RF\ \ \\
\hline \\[-3mm]
\ \ \ $\sigma$ (\AA)  & 3.051 & 3.035 & 3.019\\
\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\
\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\
\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\
\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\
\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\
\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\
\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96}
\end{tabular}
\label{params}
\end{table}

\begin{figure} 
\includegraphics[width=85mm]{gofrCompare.epsi} 
\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E 
and SSD without reaction field (top), as well as SSD/RF and SSD with
reaction field turned on (bottom). The insets show the respective
first peaks in detail. Solid Line - experiment, dashed line - SSD/E
and SSD/RF, and dotted line - SSD (with and without reaction field).}
\label{grcompare} 
\end{figure} 

\begin{figure} 
\includegraphics[width=85mm]{dualsticky.ps} 
\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \& 
SSD/RF (right). Light areas correspond to the tetrahedral attractive
part, and the darker areas correspond to the dipolar repulsive part.}
\label{isosurface} 
\end{figure} 

In the paper detailing the development of SSD, Liu and Ichiye placed
particular emphasis on an accurate description of the first solvation
shell. This resulted in a somewhat tall and sharp first peak that
integrated to give similar coordination numbers to the experimental
data obtained by Soper and Phillips.\cite{Ichiye96,Soper86} New
experimental x-ray scattering data from the Head-Gordon lab indicates
a slightly lower and shifted first peak in the g$_\mathrm{OO}(r)$, so
adjustments to SSD were made while taking into consideration the new
experimental findings.\cite{Head-Gordon00_1} Figure \ref{grcompare}
shows the relocation of the first peak of the oxygen-oxygen
g(\emph{r}) by comparing the original SSD (with and without reaction
field), SSD-E, and SSD-RF to the new experimental results. Both the
modified water models have shorter peaks that are brought in more
closely to the experimental peak (as seen in the insets of figure
\ref{grcompare}). This structural alteration was accomplished by a
reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
of the sticky potential strength and cutoffs. The cutoffs for the
tetrahedral attractive and dipolar repulsive terms were nearly swapped
with each other. Isosurfaces of the original and modified sticky
potentials are shown in figure \cite{isosurface}. In these
isosurfaces, it is easy to see how altering the cutoffs changes the
repulsive and attractive character of the particles. With a reduced
repulsive surface (the darker region), the particles can move closer
to one another, increasing the density for the overall system. This
change in interaction cutoff also results in a more gradual
orientational motion by allowing the particles to maintain preferred
dipolar arrangements before they begin to feel the pull of the
tetrahedral restructuring. Upon moving closer together, the dipolar
repulsion term becomes active and excludes the unphysical
arrangements. This compares with the original SSD's excluding dipolar
before the particles feel the pull of the ``hydrogen bonds''. Aside
from improving the shape of the first peak in the g(\emph{r}), this
improves the densities considerably by allowing the persistence of
full dipolar character below the previous 4.0 \AA\ cutoff.

While adjusting the location and shape of the first peak of
g(\emph{r}) improves the densities to some degree, these changes alone
are insufficient to bring the system densities up to the values
observed experimentally. To finish bringing up the densities, the
dipole moments were increased in both the adjusted models. Being a
dipole based model, the structure and transport are very sensitive to
changes in the dipole moment. The original SSD simply used the dipole
moment calculated from the TIP3P water model, which at 2.35 D is
significantly greater than the experimental gas phase value of 1.84
D. The larger dipole moment is a more realistic value and improve the
dielectric properties of the fluid. Both theoretical and experimental
measurements indicate a liquid phase dipole moment ranging from 2.4 D
to values as high as 3.11 D, so there is quite a range available for
adjusting the dipole
moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of
the dipole moments to 2.418 and 2.48 D for SSD/E and SSD/RF
respectively is moderate in the range of the experimental values;
however, it leads to significant changes in the density and transport
of the water models.

In order to demonstrate the benefits of this reparameterization, a
series of NPT and NVE simulations were performed to probe the density
and transport properties of the adapted models and compare the results
to the original SSD model. This comparison involved full NPT melting
sequences for both SSD/E and SSD/RF, as well as NVE transport
calculations at both self-consistent and experimental
densities. Again, the results come from five separate simulations of
1024 particle systems, and the melting sequences were started from
different ice $I_h$ crystals constructed as stated previously. Like
before, all of the NPT simulations were equilibrated for 100 ps before
a 200 ps data collection run, and they used the previous temperature's
final configuration as a starting point. All of the NVE simulations
had the same thermalization, equilibration, and data collection times
stated earlier in this paper.

\begin{figure} 
\includegraphics[width=85mm]{ssdecompare.epsi} 
\caption{Comparison of densities calculated with SSD/E to SSD without a 
reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
includes error bars, and the calculated results from the other
references were removed for clarity.}
\label{ssdedense} 
\end{figure} 

Figure \ref{ssdedense} shows the density profile for the SSD/E water
model in comparison to the original SSD without a reaction field,
experiment, and the other common water models considered
previously. The calculated densities have increased significantly over
the original SSD model and match the experimental value just below 298
K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
compares well with the experimental value of 0.997 g/cm$^3$ and is
considerably better than the SSD value of 0.967$\pm$0.003
g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
out the curve at higher temperatures, only the improvement is marginal
at best. This steep drop in densities is due to the dipolar rather
than charge based interactions which decay more rapidly at longer
distances.
 
By monitoring C$\text{p}$ throughout these simulations, the melting
transition for SSD/E was observed at 230 K, about 5 degrees lower than
SSD. The resulting density maximum is located at 240 K, again about 5
degrees lower than the SSD value of 245 K. Though there is a decrease
in both of these values, the corrected densities near room temperature
justify the modifications taken.

\begin{figure} 
\includegraphics[width=85mm]{ssdrfcompare.epsi} 
\caption{Comparison of densities calculated with SSD/RF to SSD with a 
reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
includes error bars, and the calculated results from the other
references were removed for clarity.}
\label{ssdrfdense} 
\end{figure} 

Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
SSD with an active reaction field. Like in the simulations of SSD/E,
the densities show a dramatic increase over normal SSD. At 298 K,
SSD/RF has a density of 0.997$\pm$0.001 g/cm$^3$, right in line with
experiment and considerably better than the SSD value of
0.941$\pm$0.001 g/cm$^3$. The melting point is observed at 240 K,
which is 5 degrees lower than SSD with a reaction field, and the
density maximum at 255 K, again 5 degrees lower than SSD. The density
at higher temperature still drops off more rapidly than the charge
based models but is in better agreement than SSD/E.

The reparameterization of the SSD water model, both for use with and
without an applied long-range correction, brought the densities up to
what is expected for simulating liquid water. In addition to improving
the densities, it is important that particle transport be maintained
or improved. Figure \ref{ssdediffuse} compares the temperature
dependence of the diffusion constant of SSD/E to SSD without an active
reaction field, both at the densities calculated at 1 atm and at the
experimentally calculated densities for super-cooled and liquid
water. In the upper plot, the diffusion constant for SSD/E is
consistently a little faster than experiment, while SSD starts off
slower than experiment and crosses to merge with SSD/E at high
temperatures. Both models follow the experimental trend well, but
diffuse too rapidly at higher temperatures. This abnormally fast
diffusion is caused by the decreased system density. Since the
densities of SSD/E don't deviate as much from experiment as those of
SSD, it follows the experimental trend more closely. This observation
is backed up by looking at the lower plot. The diffusion constants for
SSD/E track with the experimental values while SSD deviates on the low
side of the trend with increasing temperature. This is again a product
of SSD/E having densities closer to experiment, and not deviating to
lower densities with increasing temperature as rapidly.

\begin{figure} 
\includegraphics[width=85mm]{ssdediffuse.epsi} 
\caption{Plots of the diffusion constants calculated from SSD/E and SSD,
 both without a reaction field along with experimental results from
 Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
 upper plot is at densities calculated from the NPT simulations at a
 pressure of 1 atm, while the lower plot is at the experimentally
 calculated densities.}
\label{ssdediffuse} 
\end{figure} 

\begin{figure} 
\includegraphics[width=85mm]{ssdrfdiffuse.epsi} 
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD,
 both with an active reaction field along with experimental results
 from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
 upper plot is at densities calculated from the NPT simulations at a
 pressure of 1 atm, while the lower plot is at the experimentally
 calculated densities.}
\label{ssdrfdiffuse} 
\end{figure} 

In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
compared with SSD with an active reaction field. In the upper plot,
SSD/RF tracks with the experimental results incredibly well, identical
within error throughout the temperature range and only showing a
slight increasing trend at higher temperatures. SSD also tracks
experiment well, only it tends to diffuse a little more slowly at low
temperatures and deviates to diffuse too rapidly at high
temperatures. As was stated in the SSD/E comparisons, this deviation
away from the ideal trend is due to a rapid decrease in density at
higher temperatures. SSD/RF doesn't suffer from this problem as much
as SSD, because the calculated densities are more true to
experiment. This is again emphasized in the lower plot, where SSD/RF
tracks the experimental diffusion exactly while SSD's diffusion
constants are slightly too low due to its need for a lower density at
the specified temperature.

\subsection{Additional Observations}

While performing the melting sequences of SSD/E, some interesting
observations were made. After melting at 230 K, two of the systems
underwent crystallization events near 245 K. As the heating process
continued, the two systems remained crystalline until finally melting
between 320 and 330 K. These simulations were excluded from the data
set shown in figure \ref{ssdedense} and replaced with two additional
melting sequences that did not undergo this anomalous phase
transition, while this crystallization event was investigated
separately.

\begin{figure}
\includegraphics[width=85mm]{povIce.ps} 
\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
\label{weirdice} 
\end{figure}

The final configurations of these two melting sequences shows an
expanded zeolite-like crystal structure that does not correspond to
any known form of ice. For convenience and to help distinguish it from
the experimentally observed forms of ice, this crystal structure will
henceforth be referred to as ice-zero (ice 0). The crystallinity was
extensive enough than a near ideal crystal structure could be
obtained. Figure \ref{weirdice} shows the repeating crystal structure
of a typical crystal at 5 K. The unit cell contains eight molecules,
and figure \ref{unitcell} shows a unit cell built from the water
particle center of masses that can be used to construct a repeating
lattice, similar to figure \ref{weirdice}. Each molecule is hydrogen
bonded to four other water molecules; however, the hydrogen bonds are
flexed rather than perfectly straight. This results in a skewed
tetrahedral geometry about the central molecule. Looking back at
figure \ref{isosurface}, it is easy to see how these flexed hydrogen
bonds are allowed in that the attractive regions are conical in shape,
with the greatest attraction in the central region. Though not ideal,
these flexed hydrogen bonds are favorable enough to stabilize an
entire crystal generated around them. In fact, the imperfect ice 0
crystals were so stable that they melted at greater than room
temperature.

\begin{figure}
\includegraphics[width=65mm]{ice0cell.eps} 
\caption{Simple unit cell for constructing ice 0. In this cell, $c$ is 
equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
\label{unitcell} 
\end{figure}

The initial simulations indicated that ice 0 is the preferred ice
structure for at least SSD/E. To verify this, a comparison was made
between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at
constant pressure with SSD/E, SSD/RF, and SSD. Near ideal versions of
the three types of crystals were cooled to ~1 K, and the potential
energies of each were compared using all three water models. With
every water model, ice 0 turned out to have the lowest potential
energy: 5\% lower than $I_h$ with SSD, 6.5\% lower with SSD/E, and
7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
was observed to be ~2\% less stable than ice $I_h$. In addition to
having the lowest potential energy, ice 0 was the most expanded of the
three ice crystals, ~5\% less dense than ice $I_h$ with all of the
water models. In all three water models, ice $I_c$ was observed to be
~2\% more dense than ice $I_h$. 

In addition to the low temperature comparisons, melting sequences were
performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
and SSD both with and without a reaction field. The melting
transitions for both SSD/E and SSD without a reaction field occurred
at temperature in excess of 375 K. SSD/RF and SSD with a reaction
field had more reasonable melting transitions, down near 325 K. These
melting point observations emphasize how preferred this crystal
structure is over the most common types of ice when using these single
point water models.

Recognizing that the above tests show ice 0 to be both the most stable
and lowest density crystal structure for these single point water
models, it is interesting to speculate on the favorability of this
crystal structure with the different charge based models. As a quick
test, these 3 crystal types were converted from SSD type particles to
TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
minimizations were performed on all of these crystals to compare the
system energies. Again, ice 0 was observed to have the lowest total
system energy. The total energy of ice 0 was ~2\% lower than ice
$I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial
results, we would not be surprised if results from the other common
water models show ice 0 to be the lowest energy crystal structure. A
continuation on work studing ice 0 with multipoint water models will
be published in a coming article.

\section{Conclusions}
The density maximum and temperature dependent transport for the SSD
water model, both with and without the use of reaction field, were
studied via a series of NPT and NVE simulations. The constant pressure
simulations of the melting of both $I_h$ and $I_c$ ice showed a
density maximum near 260 K. In most cases, the calculated densities
were significantly lower than the densities calculated in simulations
of other water models. Analysis of particle diffusion showed SSD to
capture the transport properties of experimental very well in both the
normal and super-cooled liquid regimes. In order to correct the
density behavior, SSD was reparameterized for use both with and
without a long-range interaction correction, SSD/RF and SSD/E
respectively. In addition to correcting the abnormally low densities,
these new versions were show to maintain or improve upon the transport
and structural features of the original water model, all while
maintaining the fast performance of the SSD water model. This work
shows these simple water models, and in particular SSD/E and SSD/RF,
to be excellent choices to represent explicit water in future
simulations of biochemical systems.

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

\bibliographystyle{jcp}

\bibliography{nptSSD} 

%\pagebreak 

\end{document}
Revision:	777
Committed:	Fri Sep 19 19:29:24 2003 UTC (20 years, 9 months ago) by chrisfen
Content type:	application/x-tex
File size:	48772 byte(s)
Log Message:	I must have changed something, but for the life of me I can't recall it...