trunk/chrisDissertation/Water.tex

\chapter{\label{chap:water}SIMPLE MODELS FOR WATER}

One of the most important tasks in the simulation of biochemical
systems is the proper depiction of the aqueous environment around the
molecules of interest.  In some cases (such as in the simulation of
phospholipid bilayers), the majority of the calculations that are
performed involve interactions with or between solvent molecules.
Thus, the motion and behavior of molecules in biochemical simulations
are highly dependent on the properties of the water model that is
chosen as the solvent.
\begin{figure}
\includegraphics[width=\linewidth]{./figures/waterModels.pdf}
\caption{Partial-charge geometries for the TIP3P, TIP4P, TIP5P, and 
SPC/E rigid body water models.\cite{Jorgensen83,Mahoney00,Berendsen87}
In the case of the TIP models, the depiction of water improves with 
increasing number of point charges; however, computational performance 
simultaneously degrades due to the increasing number of distances and 
interactions that need to be calculated.}
\label{fig:waterModels}
\end{figure}

As discussed in the previous chapter, water is typically modeled with
fixed geometries of point charges shielded by the repulsive part of a
Lennard-Jones interaction. Some of the common water models are shown
in figure \ref{fig:waterModels}. The various models all have their
benefits and drawbacks, and these primarily focus on the balance
between computational efficiency and the ability to accurately predict
the properties of bulk water. For example, the TIP5P model improves on
the structural and transport properties of water relative to the TIP3P
and TIP4P models, yet this comes at a greater than 50\% increase in
computational cost.\cite{Mahoney00,Mahoney01} This cost is entirely
due to the additional distance and electrostatic calculations that
come from the extra point charges in the model description. Thus, the
main criteria for choosing a water model are:

\begin{enumerate}
\item capturing the liquid state properties and 
\item having the fewest number of points to insure efficient performance. 
\end{enumerate}
As researchers have begun to study larger systems, such as entire
viruses, the choice has shifted towards efficiency over accuracy in
order to make the calculations feasible.\cite{Freddolino06}

\section{Soft Sticky Dipole Model for Water}

One recently-developed model that largely succeeds in retaining the
accuracy of bulk properties while greatly reducing the computational
cost is the Soft Sticky Dipole (SSD) water
model.\cite{Liu96,Liu96b,Chandra99,Tan03} The SSD model was developed
as a modified form of the hard-sphere water model proposed by Bratko,
Blum, and Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point}
model which has an interaction site that is both a point dipole and a
Lennard-Jones core.  However, since the normal aligned and
anti-aligned geometries favored by point dipoles are poor mimics of
local structure in liquid water, a short ranged ``sticky'' potential
is also added.  The sticky potential directs the molecules to assume
the proper hydrogen bond orientation in the first solvation shell.

The interaction between two SSD water molecules \emph{i} and \emph{j}
is given by the potential
\begin{equation}
u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)\ +
u_{ij}^{sp}
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j),
\end{equation}
where the ${\bf r}_{ij}$ is the position vector between molecules
\emph{i} and \emph{j} with magnitude $r_{ij}$, and
${\bf \Omega}_i$ and ${\bf \Omega}_j$ represent the orientations of
the two molecules. The Lennard-Jones and dipole interactions are given
by the following familiar forms:
\begin{equation}
u_{ij}^{LJ}(r_{ij}) = 4\epsilon
\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]
\ ,
\end{equation}
and
\begin{equation}
u_{ij}^{dp} = \frac{|\mu_i||\mu_j|}{4 \pi \epsilon_0 r_{ij}^3} \left(
\hat{\bf u}_i \cdot \hat{\bf u}_j - 3(\hat{\bf u}_i\cdot\hat{\bf
r}_{ij})(\hat{\bf u}_j\cdot\hat{\bf r}_{ij}) \right)\ ,
\end{equation}
where $\hat{\bf u}_i$ and $\hat{\bf u}_j$ are the unit vectors along
the dipoles of molecules $i$ and $j$ respectively. $|\mu_i|$ and
$|\mu_j|$ are the strengths of the dipole moments, and $\hat{\bf
r}_{ij}$ is the unit vector pointing from molecule $j$ to molecule
$i$.

The sticky potential is somewhat less familiar:
\begin{equation}
u_{ij}^{sp}
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)
+ s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf
\Omega}_j)]\ .
\label{eq:stickyfunction}
\end{equation}
Here, $\nu_0$ is a strength parameter for the sticky potential, and
$s$ and $s^\prime$ are cubic switching functions which turn off the
sticky interaction beyond the first solvation shell. The $w$ function
can be thought of as an attractive potential with tetrahedral
geometry:
\begin{equation}
w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
\end{equation}
while the $w^\prime$ function counters the normal aligned and
anti-aligned structures favored by point dipoles:
\begin{equation}
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^\circ,
\end{equation}
It should be noted that $w$ is proportional to the sum of the $Y_3^2$
and $Y_3^{-2}$ spherical harmonics (a linear combination which
enhances the tetrahedral geometry for hydrogen bonded structures),
while $w^\prime$ is a purely empirical function.  A more detailed
description of the functional parts and variables in this potential
can be found in the original SSD
articles.\cite{Liu96,Liu96b,Chandra99,Tan03}

Since SSD is a single-point {\it dipolar} model, the force
calculations are simplified significantly relative to the standard
{\it charged} multi-point models. In the original Monte Carlo
simulations using this model, Liu and Ichiye reported that using SSD
decreased computer time by a factor of 6-7 compared to other
models.\cite{Liu96b} What is most impressive is that this savings did
not come at the expense of accurate depiction of the liquid state
properties.  Indeed, SSD maintains reasonable agreement with the Soper
data for the structural features of liquid water.\cite{Soper86,Liu96b}
Additionally, the dynamical properties exhibited by SSD agree with
experiment better than those of more computationally expensive models
(like TIP3P and SPC/E).\cite{Chandra99} The combination of speed and
accurate depiction of solvent properties makes SSD a very attractive
model for the simulation of large scale biochemical simulations.

It is important to note that the SSD model was originally developed
for use with the Ewald summation for handling long-range
electrostatics.\cite{Ewald21} In applying this water model in a
variety of molecular systems, it would be useful to know its
properties and behavior under the more computationally efficient
reaction field (RF) technique, the correction techniques discussed in
the previous chapter, or even a simple
cutoff.\cite{Onsager36,Fennell06} This chapter addresses these issues
by looking at the structural and transport behavior of SSD over a
variety of temperatures with the purpose of utilizing the RF
correction technique.  We then suggest modifications to the parameters
that result in more realistic bulk phase behavior.  It should be noted
that in a recent publication, some of the original investigators of
the SSD water model have suggested adjustments to the SSD water model
to address abnormal density behavior (also observed here), calling the
corrected model SSD1.\cite{Tan03} In the later sections of this
chapter, we compare our modified variants of SSD with both the
original SSD and SSD1 models and discuss how our changes improve the
depiction of water.

\section{Simulation Methods}{\label{sec:waterSimMethods}

Most of the calculations in this particular study were performed using
an internally developed simulation code that was one of the precursors
of the {\sc oopse} molecular dynamics (MD) package.\cite{Meineke05}
All of the capabilities of this code have been incorporated into {\sc
oopse}, and calculation results are consistent between the two
simulation packages. The later calculations involving the damped
shifted force ({\sc sf}) techniques were performed using {\sc oopse}.

In the primary simulations of this study, long-range dipole-dipole
corrections were accounted for by using either the reaction field
technique or a simple cubic switching function at the cutoff
radius. Interestingly, one of the early applications of a reaction
field was in Monte Carlo simulations of liquid water.\cite{Barker73}
In this method, the magnitude of the reaction field acting on dipole
$i$ is
\begin{equation}
\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}),
\label{eq: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, ${\bf\mu}_{j}$ is the dipole moment vector of particle $j$,
and $s(r_{ij})$ is a cubic switching function.\cite{Allen87} The
reaction field contribution to the total energy from particle $i$ is
given by $-\frac{1}{2}{\bf\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque
on dipole $i$ by ${\bf\mu}_{i}\times\mathcal{E}_{i}$.\cite{Allen87} An
applied reaction field will alter the bulk orientational properties of
simulated water, and there is particular sensitivity of these
properties on changes in the length of the cutoff
radius.\cite{vanderSpoel98} This behavior makes the reaction field a
less attractive method than the Ewald sum; however, for very large
systems, the computational benefit of reaction field is significant.
 
In contrast to the simulations with a reaction field, we have also
performed a companion set of simulations {\it without} a surrounding
dielectric (i.e. using a simple cubic switching function at the cutoff
radius). As a result, we have developed two reparametrizations of SSD
which can be used either with or without an active reaction field.

To determine the preferred densities of the models, we performed
simulations in the isobaric-isothermal ({\it NPT}) ensemble. All
dynamical properties for these models were then obtained from
microcanonical ({\it NVE}) simulations done at densities matching the
{\it NPT} density for a particular target temperature.  The constant
pressure simulations were implemented using an integral thermostat and
barostat as outlined by Hoover.\cite{Hoover85,Hoover86} All molecules
were treated as non-linear rigid bodies. Vibrational constraints are
not necessary in simulations of SSD, because there are no explicit
hydrogen atoms, and thus no molecular vibrational modes need to be
considered.

The symplectic splitting method proposed by Dullweber, Leimkuhler, and
McLachlan ({\sc dlm}, see section \ref{sec:IntroIntegrate}) was used
to carry out the integration of the equations of motion in place of
the more prevalent quaternion
method.\cite{Dullweber97,Evans77,Evans77b,Allen87} The reason behind
this decision was that, in {\it NVE} simulations, the energy drift
when using quaternions was substantially greater than when using the
{\sc dlm} method (Fig. \ref{fig:timeStepIntegration}). This steady
drift in the total energy has also been observed in other
studies.\cite{Kol97}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/timeStepIntegration.pdf}
\caption{Energy conservation using both quaternion-based integration 
and the {\sc dlm} method with increasing time step. The larger time
step plots are shifted from the true energy baseline (that of $\Delta
t$ = 0.1~fs) for clarity.}
\label{fig:timeStepIntegration}
\end{figure}

The {\sc dlm} method allows for Verlet style integration orientational
motion of rigid bodies via a sequence of rotation matrix
operations. Because these matrix operations are more costly than the
simpler arithmetic operations for quaternion propagation, typical SSD
particle simulations using {\sc dlm} are 5-10\% slower than
simulations using the quaternion method and an identical time
step. This additional expense is justified because of the ability to
use time steps that are more that twice as long and still achieve the
same energy conservation.

Figure \ref{fig:timeStepIntegration} shows the resulting energy drift
at various time steps for both {\sc dlm} and quaternion
integration. All of the 1000 SSD particle simulations started with the
same configuration, and the only difference was the method used to
handle orientational motion. At time steps of 0.1 and 0.5~fs, both
methods for propagating the orientational degrees of freedom conserve
energy fairly well, with the quaternion method showing a slight energy
drift over time in the 0.5~fs time step simulation. Time steps of 1 and
2~fs clearly demonstrate the benefits in energy conservation that come
with the {\sc dlm} method. 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 drifts in water simulations using {\sc dlm} integration were
unnoticeable for time steps up to 3~fs. We observed a slight energy
drift on the order of 0.012~kcal/mol per nanosecond with a time step
of 4~fs. As expected, this drift increases dramatically with increasing
time step. To insure accuracy in our {\it NVE} simulations, time steps
were set at 2~fs and were also kept at this value for {\it NPT}
simulations.

Proton-disordered ice crystals in both the I$_\textrm{h}$ and
I$_\textrm{c}$ lattices were generated as starting points for all
simulations. The I$_\textrm{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$_\textrm{h}$ ice, the simulation boxes were
orthorhombic in shape with an edge length ratio of approximately
1.00$\times$1.06$\times$1.23. We then allowed the particles to orient
freely about their fixed lattice positions with angular momenta values
randomly sampled at 400~K. 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$_\textrm{h}$ ice crystals that obey the
Bernal-Fowler rules.\cite{Bernal33,Rahman72} This method was also
utilized in the making of diamond lattice I$_\textrm{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{SSD Density Behavior}

Melting studies were performed on the randomized ice crystals using
the {\it NPT} ensemble. 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. It should be noted that the calculated melting temperature
($T_\textrm{m}$) will not be the true $T_\textrm{m}$ because of
super-heating due to the relatively short time scales in molecular
simulations. This behavior results in inflated $T_\textrm{m}$ values;
however, these values provide a reasonable initial estimate of
$T_\textrm{m}$. 

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 fluctuations in the volume autocorrelation function damped
out in all of the simulations in under 20~ps.

Our initial simulations focused on the original SSD water model, and
an average density versus temperature plot is shown in figure
\ref{fig:ssdDense}. Note that the density maximum when using a
reaction field appears between 255 and 265~K.  There were smaller
fluctuations in the density at 260~K than at either 255 or 265~K, so we
report this value as the location of the density maximum. Figure
\ref{fig:ssdDense} was constructed using ice I$_\textrm{h}$ crystals
for the initial configuration; though not pictured, the simulations
starting from ice I$_\textrm{c}$ crystal configurations showed similar
results, with a liquid-phase density maximum at the same temperature.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdDense.pdf}
\caption{ Density versus temperature for TIP3P, SPC/E, TIP4P, SSD,
SSD with a reaction field, and
experiment.\cite{Jorgensen98b,Baez94,CRC80}. Note that using a
reaction field lowers the density more than the already lowered SSD
densities. The lower than expected densities for the SSD model
prompted the original reparametrization of SSD to SSD1.\cite{Tan03}}
\label{fig:ssdDense}
\end{figure}

The density maximum for SSD compares quite favorably to other simple
water models. Figure \ref{fig:ssdDense} also shows calculated
densities of several other models and experiment obtained from other
sources.\cite{Jorgensen98b,Baez94,CRC80} Of the water models shown,
SSD has a temperature closest to the experimentally observed density
maximum. Of the {\it charge-based} models in figure
\ref{fig:ssdDense}, TIP4P has a density maximum behavior most like
that seen in SSD. Though not included in this plot, it is useful to
note that TIP5P has a density maximum nearly identical to the
experimentally measured temperature (see section
\ref{sec:t5peDensity}.

Liquid state densities in water have been observed to be dependent on
the cutoff radius ($R_\textrm{c}$), both with and without the use of a
reaction field.\cite{vanderSpoel98} In order to address the possible
effect of $R_\textrm{c}$, simulations were performed with a cutoff
radius of 12~\AA\, complementing the 9~\AA\ $R_\textrm{c}$ used in the
previous SSD simulations. All of the resulting densities overlapped
within error and showed no significant trend toward lower or higher
densities in simulations both with and without reaction field. 

The key feature to recognize in figure \ref{fig:ssdDense} is the
density scaling of SSD relative to other common models at any given
temperature. SSD assumes a lower density than any of the other listed
models at the same pressure. This behavior 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 reduced
densities is the presence of the reaction
field.\cite{vanderSpoel98,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{fig:ssdDense}. Without
the reaction field, the densities increase to more experimentally
reasonable values, especially around the freezing point of liquid
water. The shape of the curve is similar to the curve produced from
SSD simulations using reaction field, specifically the rapidly
decreasing densities at higher temperatures; however, a shift in the
density maximum location, down to 245~K, is observed. This is a more
accurate comparison to the other listed water models, in that no long
range corrections were applied in those
simulations.\cite{Baez94,Jorgensen98b} However, even without the
reaction field, the density around 300~K is still significantly lower
than experiment and comparable water models. This anomalous behavior
was what lead Tan {\it et al.} to recently reparametrize
SSD.\cite{Tan03} Throughout the remainder of the paper our
reparametrizations of SSD will be compared with their newer SSD1
model.

\section{SSD Transport Behavior}

Accurate dynamical properties of a water model are particularly
important when using the model to study permeation or transport across
biological membranes.  In order to probe transport in bulk water, {\it
NVE} simulations were performed at the average densities obtained from
the {\it NPT} simulations at an identical target
temperature. Simulations started with randomized velocities and
underwent 50~ps of temperature scaling and 50~ps of constant energy
equilibration before a 200~ps data collection run. Diffusion constants
were calculated via linear fits to the long-time behavior of the
mean-square displacement as a function of time.\cite{Allen87} The
averaged results from five sets of {\it NVE} simulations are displayed
in figure \ref{fig:ssdDiffuse}, alongside experimental, SPC/E, and TIP5P
results.\cite{Gillen72,Holz00,Baez94,Mahoney01} 

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdDiffuse.pdf} 
\caption{ Average self-diffusion constant as a function of temperature for 
SSD, SPC/E, and TIP5P compared with experimental
data.\cite{Baez94,Mahoney01,Gillen72,Holz00} Of the three water
models shown, SSD has the least deviation from the experimental
values. The rapidly increasing diffusion constants for TIP5P and SSD
correspond to significant decreases in density at the higher
temperatures.}
\label{fig:ssdDiffuse}
\end{figure} 

The observed values for the diffusion constant point out one of the
strengths of the SSD model. Of the three models shown, the SSD model
has the most accurate depiction of self-diffusion in both the
supercooled and liquid regimes.  SPC/E does a respectable job by
reproducing values similar to experiment around 290~K; however, it
deviates at both higher and lower temperatures, failing to predict the
correct thermal trend. TIP5P and SSD both start off low at colder
temperatures and tend to diffuse too rapidly at higher temperatures.
This behavior at higher temperatures is not particularly surprising
since the densities of both TIP5P and SSD are lower than experimental
water densities at higher temperatures.  When calculating the
diffusion coefficients for SSD at experimental densities (instead of
the densities from the {\it NPT} simulations), the resulting values
fall more in line with experiment at these temperatures.

\section{Structural Changes and Characterization}

By starting the simulations from the crystalline state, we can get an
estimation of the $T_\textrm{m}$ of the ice structure, and beyond the
melting point, we study the phase behavior of the liquid. The constant
pressure heat capacity ($C_\textrm{p}$) was monitored to locate
$T_\textrm{m}$ in each of the simulations. In the melting simulations
of the 1024 particle ice I$_\textrm{h}$ simulations, a large spike in
$C_\textrm{p}$ occurs at 245~K, indicating a first order phase
transition for the melting of these ice crystals (see figure
\ref{fig:ssdCp}. When the reaction field is turned off, the melting
transition occurs at 235~K.  These melting transitions are considerably
lower than the experimental value of 273~K, indicating that the solid
ice I$_\textrm{h}$ is not thermodynamically preferred relative to the
liquid state at these lower temperatures.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdCp.pdf}
\caption{Heat capacity versus temperature for the SSD model with an 
active reaction field. Note the large spike in $C_p$ around 245~K,
indicating a phase transition from the ordered crystal to disordered
liquid.}
\label{fig:ssdCp}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/fullContour.pdf} 
\caption{ Contour plots of 2D angular pair correlation functions for
512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
signify regions of enhanced density while light areas signify
depletion relative to the bulk density. White areas have pair
correlation values below 0.5 and black areas have values above 1.5.}
\label{fig:contour} 
\end{figure} 

\begin{figure}
\centering
\includegraphics[width=2.5in]{./figures/corrDiag.pdf} 
\caption{ An illustration of angles involved in the correlations 
displayed in figure \ref{fig:contour}.}
\label{fig:corrAngle}
\end{figure}

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

\begin{equation}
g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
\end{equation}
\begin{equation}
g_{\text{AB}}(r,\cos\omega) = 
\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
\end{equation}
where $\theta$ and $\omega$ refer to the angles shown in figure
\ref{fig:corrAngle}. By binning over both distance and the cosine of the
desired angle between the two dipoles, the $g(r)$ can be analyzed to
determine the common dipole arrangements that constitute the peaks and
troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
of figure \ref{fig:contour} show results from an ice I$_\textrm{c}$
simulation. The first peak in the $g(r)$ consists primarily of the
preferred hydrogen bonding arrangements as dictated by the tetrahedral
sticky potential - one peak for the hydrogen bond donor and the other
for the hydrogen bond acceptor.  Due to the high degree of
crystallinity of the sample, the second and third solvation shells
show a repeated peak arrangement which decays at distances around the
fourth solvation shell, near the imposed cutoff for the Lennard-Jones
and dipole-dipole interactions.  In the higher temperature simulation
shown in frames C and D, these long-range features deteriorate
rapidly. The first solvation shell still shows the strong effect of
the sticky-potential, although it covers a larger area, extending to
include a fraction of aligned dipole peaks within the first solvation
shell. The latter peaks lose due to thermal motion and as the
competing dipole force overcomes the sticky potential's tight
tetrahedral structuring of the crystal.

This complex interplay between dipole and sticky interactions was
remarked upon as a possible reason for the split second peak in the
oxygen-oxygen pair correlation function,
$g_\textrm{OO}(r)$.\cite{Liu96b} At low temperatures, the second
solvation shell peak appears to have two distinct components that
blend together to form one observable peak. At higher temperatures,
this split character alters to show the leading 4~\AA\ peak dominated
by equatorial anti-parallel dipole orientations. There is also a
tightly bunched group of axially arranged dipoles that most likely
consist of the smaller fraction of aligned dipole pairs. The trailing
component of the split peak at 5~\AA\ is dominated by aligned dipoles
that assume hydrogen bond arrangements similar to those seen in the
first solvation shell. This evidence indicates that the dipole pair
interaction begins to dominate outside of the range of the dipolar
repulsion term.  The energetically favorable dipole arrangements
populate the region immediately outside this repulsion region (around
4~\AA), while arrangements that seek to satisfy both the sticky and
dipole forces locate themselves just beyond this initial buildup
(around 5~\AA).

This analysis indicates that the split second peak is primarily the
product of the dipolar repulsion term of the sticky potential. In
fact, the inner peak can be pushed out and merged with the outer split
peak just by extending the switching function ($s^\prime(r_{ij})$)
from its normal 4~\AA\ cutoff to values of 4.5 or even 5~\AA. This
type of correction is not recommended for improving the liquid
structure, since the second solvation shell would still be shifted too
far out. In addition, this would have an even more detrimental effect
on the system densities, leading to a liquid with a more open
structure and a density considerably lower than the already low SSD
density.  A better correction would be to include the
quadrupole-quadrupole interactions for the water particles outside of
the first solvation shell, but this would remove the simplicity and
speed advantage of SSD.

\section{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 reparametrization, it is
important not to focus on just one property and neglect the others. In
this case, it would be ideal to correct the densities while
maintaining the accurate transport behavior.

The parameters available for tuning include the $\sigma$ and
$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the
strength of the sticky potential ($\nu_0$), and the cutoff distances
for the sticky attractive and dipole repulsive cubic switching
function cutoffs ($r_l$, $r_u$ and $r_l^\prime$, $r_u^\prime$
respectively). The results of the reparametrizations are shown in
table \ref{tab:ssdParams}. We are calling these reparametrizations the
Soft Sticky Dipole Reaction Field (SSD/RF - for use with a reaction
field) and Soft Sticky Dipole Extended (SSD/E - an attempt to improve
the liquid structure in simulations without a long-range correction).

\begin{table}
\caption{PARAMETERS FOR THE ORIGINAL AND ADJUSTED SSD MODELS}

\centering
\begin{tabular}{ llcccc } 
\toprule
\toprule
Parameters & & SSD & SSD1 & SSD/E & SSD/RF \\
\midrule
$\sigma$ & (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
$\epsilon$ & (kcal mol$^{-1}$) & 0.152 & 0.152 & 0.152 & 0.152\\
$\mu$ & (D) & 2.35 & 2.35 & 2.42 & 2.48\\
$\nu_0$ & (kcal mol$^{-1}$) & 3.7284 & 3.6613 & 3.90 & 3.90\\
$\omega^\circ$ & & 0.07715 & 0.07715 & 0.07715 & 0.07715\\
$r_l$ & (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
$r_u$ & (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
$r_l^\prime$ & (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\
$r_u^\prime$ & (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
\bottomrule
\end{tabular}
\label{tab:ssdParams}
\end{table}

\begin{figure}
\centering 
\includegraphics[width=4.5in]{./figures/newGofRCompare.pdf} 
\caption{ Plots showing the experimental $g(r)$ (Ref. \cite{Hura00}) 
with SSD/E and SSD1 without reaction field (top), as well as SSD/RF
and SSD1 with reaction field turned on (bottom). The changes in
parameters have lowered and broadened the first peak of SSD/E and
SSD/RF, resulting in a better fit to the first solvation shell.}
\label{fig:gofrCompare} 
\end{figure} 

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/dualPotentials.pdf} 
\caption{ Positive and negative isosurfaces of the sticky potential for
SSD and SSD1 (A) and SSD/E \& SSD/RF (B). Gold areas correspond to the
tetrahedral attractive component, and blue areas correspond to the
dipolar repulsive component.}
\label{fig:isosurface} 
\end{figure} 

In the original paper detailing the development of SSD, Liu and Ichiye
placed particular emphasis on an accurate description of the first
solvation shell. This resulted in a somewhat tall and narrow first
peak in $g(r)$ that integrated to give similar coordination numbers to
the experimental data obtained by Soper and
Phillips.\cite{Liu96b,Soper86} New experimental x-ray scattering data
from Hura {\it et al.} indicates a slightly lower and shifted first
peak in the $g_\textrm{OO}(r)$, so our adjustments to SSD were made
after taking into consideration the new experimental
findings.\cite{Hura00} Figure \ref{fig:gofrCompare} shows the
relocation of the first peak of the oxygen-oxygen $g(r)$ by comparing
the revised SSD model (SSD1), SSD/E, and SSD/RF to the new
experimental results. Both modified water models have shorter peaks
that match more closely to the experimental peak (as seen in the
insets of figure \ref{fig:gofrCompare}).  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{tab:ssdParams}, 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{fig: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 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(r)$, this modification improves the densities considerably by
allowing the persistence of full dipolar character below the previous
4~\AA\ cutoff.

While adjusting the location and shape of the first peak of $g(r)$
improves the densities, these changes alone are insufficient to bring
the system densities up to the values observed experimentally.  To
further increase the densities, the dipole moments were increased in
both of our adjusted models. Since SSD is a dipole based model, the
structure and transport are very sensitive to changes in the dipole
moment. The original SSD simply used the dipole moment calculated from
the TIP3P water model, which at 2.35~D is significantly greater than
the experimental gas phase value of 1.84~D. The larger dipole moment
is a more realistic value and improves the dielectric properties of
the fluid. Both theoretical and experimental measurements indicate a
liquid phase dipole moment ranging from 2.4~D to values as high as
3.11~D, providing a substantial range of reasonable values for a
dipole moment.\cite{Sprik91,Gubskaya02,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.

\subsection{Density Behavior}

In order to demonstrate the benefits of these reparametrizations, we
performed a series of {\it NPT} and {\it NVE} simulations to probe the
density and transport properties of the adapted models and compare the
results to the original SSD model. This comparison involved full {\it
NPT} melting sequences for both SSD/E and SSD/RF, as well as {\it NVE}
transport calculations at the calculated self-consistent
densities. Again, the results were obtained from five separate
simulations of 1024 particle systems, and the melting sequences were
started from different ice I$_\textrm{h}$ crystals constructed as
described previously. Each {\it 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 {\it NVE} simulations had the same
thermalization, equilibration, and data collection times as stated
previously.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdeDense.pdf} 
\caption{ Comparison of densities calculated with SSD/E to 
SSD1 without a reaction field, TIP3P, SPC/E, TIP5P, and
experiment.\cite{Jorgensen98b,Baez94,Mahoney00,CRC80} Both SSD1 and
SSD/E show good agreement with experiment when the long-range
correction is neglected.}
\label{fig:ssdeDense} 
\end{figure} 

Figure \ref{fig:ssdeDense} shows the density profiles for SSD/E, SSD1,
TIP3P, TIP4P, and SPC/E alongside the experimental results. The
calculated densities for both SSD/E and SSD1 have increased
significantly over the original SSD model (see figure
\ref{fig:ssdDense}) 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~g/cm$^3$ and 0.999~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~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{fig:gofrCompare}), but they have shifted the density maximum to
much lower temperatures. This comes about via an increase in the
liquid disorder through the weakening of the sticky potential and
strengthening of the dipolar character. However, this increasing
disorder in the SSD/E model has little effect on the melting
transition. By monitoring $C_p$ throughout these simulations, we
observed a melting transition for SSD/E at 235~K, the same as SSD and
SSD1.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdrfDense.pdf} 
\caption{ Comparison of densities calculated with SSD/RF to 
SSD1 with a reaction field, TIP3P, SPC/E, TIP5P, and
experiment.\cite{Jorgensen98b,Baez94,Mahoney00,CRC80} This plot
shows the benefit afforded by the reparametrization for use with a
reaction field correction - SSD/RF provides significantly more
accurate densities than SSD1 when performing room temperature
simulations.}
\label{fig:ssdrfDense} 
\end{figure} 

Including the reaction field long-range correction results is a more
interesting comparison.  A density profile including SSD/RF and SSD1
with an active reaction field is shown in figure \ref{fig: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{fig:ssdDense}). At 298~K, SSD/RF has a density of 0.997
g/cm$^3$, directly in line with experiment and considerably better
than the original SSD value of 0.941~g/cm$^3$ and the SSD1 value of
0.972~g/cm$^3$. These results further emphasize the importance of
reparametrization 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.

\subsection{Transport Behavior}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdeDiffuse.pdf} 
\caption{ The diffusion constants calculated from SSD/E and 
SSD1 (both without a reaction field) along with experimental
results.\cite{Gillen72,Holz00} The {\it NVE} calculations were
performed at the average densities from the {\it NPT} simulations for
the respective models. SSD/E is slightly more mobile than experiment
at all of the temperatures, but it is closer to experiment at
biologically relevant temperatures than SSD1 without a long-range
correction.}
\label{fig:ssdeDiffuse} 
\end{figure} 

The reparametrization 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 proper simulation of liquid water. In addition to
improving the densities, it is important that the diffusive behavior
of SSD be maintained or improved. Figure \ref{fig:ssdeDiffuse}
compares the temperature dependence of the diffusion constant of SSD/E
to SSD1 without an active reaction field at the densities calculated
from their respective {\it NPT} simulations at 1 atm. The diffusion
constant for SSD/E is consistently higher than experiment, while SSD1
remains lower than experiment until relatively high temperatures
(around 360~K). Both models follow the shape of the experimental curve
below 300~K but tend to diffuse too rapidly at higher temperatures, as
seen in SSD1 crossing above 360~K.  This increasing diffusion relative
to the experimental values is caused by the rapidly decreasing system
density with increasing temperature.  Both SSD1 and SSD/E show this
deviation in particle mobility, but this trend has different
implications on the diffusive behavior of the models.  While SSD1
shows more experimentally accurate diffusive behavior in the high
temperature regimes, SSD/E shows more accurate behavior in the
supercooled and biologically relevant temperature ranges.  Thus, the
changes made to improve the liquid structure may have had an adverse
affect on the density maximum, but they improve the transport behavior
of SSD/E relative to SSD1 under the most commonly simulated
conditions.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdrfDiffuse.pdf} 
\caption{ The diffusion constants calculated from SSD/RF and 
SSD1 (both with an active reaction field) along with experimental
results.\cite{Gillen72,Holz00} The {\it NVE} calculations were
performed at the average densities from the {\it NPT} simulations for
both of the models. SSD/RF captures the self-diffusion of water
throughout most of this temperature range. The increasing diffusion
constants at high temperatures for both models can be attributed to
lower calculated densities than those observed in experiment.}
\label{fig:ssdrfDiffuse} 
\end{figure} 

In figure \ref{fig:ssdrfDiffuse}, the diffusion constants for SSD/RF are
compared to SSD1 with an active reaction field. Note that SSD/RF
tracks the experimental results quantitatively, identical within error
throughout most of the temperature range shown and exhibiting only a
slight increasing trend at higher temperatures. SSD1 tends to diffuse
more slowly at low temperatures and deviates to diffuse too rapidly at
temperatures greater than 330~K.  As stated above, this deviation away
from the ideal trend is due to a rapid decrease in density at higher
temperatures. SSD/RF does not suffer from this problem as much as SSD1
because the calculated densities are closer to the experimental
values. These results again emphasize the importance of careful
reparametrization when using an altered long-range correction.

\subsection{Summary of Liquid State Properties}

\begin{table}
\caption{PROPERTIES OF THE SINGLE-POINT WATER MODELS COMPARED WITH
EXPERIMENTAL DATA AT AMBIENT CONDITIONS}
\footnotesize
\centering
\begin{tabular}{ llccccc } 
\toprule
\toprule
& & SSD1 & SSD/E & SSD1 (RF) & SSD/RF & Experiment [Ref.] \\
\midrule
$\rho$ & (g cm$^{-3}$) & 0.999(1) & 0.996(1) & 0.972(2) & 0.997(1) & 0.997 \cite{CRC80}\\
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 28.80(11) & 25.45(9) & 28.28(6) & 23.83(16) & 18.005 \cite{Wagner02}\\
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 1.78(7) & 2.51(18) & 2.00(17) & 2.32(6) & 2.299 \cite{Mills73}\\
$n_C$ & & 3.9 & 4.3 & 3.8 & 4.4 & 4.7 \cite{Hura00}\\ 
$n_H$ & & 3.7 & 3.6 & 3.7 & 3.7 & 3.5 \cite{Soper86}\\
$\tau_1$ & (ps) & 10.9(6) & 7.3(4) & 7.5(7) & 7.2(4) & 5.7 \cite{Eisenberg69}\\
$\tau_2$ & (ps) & 4.7(4) & 3.1(2) & 3.5(3) & 3.2(2) & 2.3 \cite{Krynicki66}\\
\bottomrule
\end{tabular}
\label{tab:liquidProperties}
\end{table}

Table \ref{tab:liquidProperties} gives a synopsis of the liquid state
properties of the water models compared in this study along with the
experimental values for liquid water at ambient conditions. The
coordination number ($n_C$) and number of hydrogen bonds per particle
($n_H$) were calculated by integrating the following relations:
\begin{equation}
n_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2g_{\textrm{OO}}(r)dr,
\end{equation}
\begin{equation}
n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2g_{\textrm{OH}}(r)dr,
\end{equation}
where $\rho$ is the number density of the specified pair interactions,
$a$ and $b$ are the radial locations of the minima following the first
peak in $g_\textrm{OO}(r)$ or $g_\textrm{OH}(r)$ respectively. The
number of hydrogen bonds stays relatively constant across all of the
models, but the coordination numbers of SSD/E and SSD/RF show an
improvement over SSD1.  This improvement is primarily due to extension
of the first solvation shell in the new parameter sets.  Because $n_H$
and $n_C$ are nearly identical in SSD1, it appears that all molecules
in the first solvation shell are involved in hydrogen bonds.  Since
$n_H$ and $n_C$ differ in the newly parameterized models, the
orientations in the first solvation shell are a bit more ``fluid''.
Therefore SSD1 over-structures the first solvation shell and our
reparametrizations have returned this shell to more realistic
liquid-like behavior.

The time constants for the orientational autocorrelation functions
are also displayed in Table \ref{tab:liquidProperties}. The dipolar
orientational time correlation functions ($C_{l}$) are described
by:
\begin{equation}
C_{l}(t) = \langle P_l[\hat{\mathbf{u}}_j(0)\cdot\hat{\mathbf{u}}_j(t)]\rangle,
\label{eq:reorientCorr}
\end{equation}
where $P_l$ are Legendre polynomials of order $l$ and
$\hat{\mathbf{u}}_j$ is the unit vector pointing along the molecular
dipole.\cite{Rahman71} Note that this is identical to equation
(\ref{eq:OrientCorr}) were $\alpha$ is equal to $z$. From these
correlation functions, the orientational relaxation time of the dipole
vector can be calculated from an exponential fit in the long-time
regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these time
constants were averaged over five detailed {\it NVE} simulations
performed at the ambient conditions for each of the respective
models. It should be noted that the commonly cited value of 1.9~ps for
$\tau_2$ was determined from the NMR data in Ref. \cite{Krynicki66} at
a temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong
temperature dependence of $\tau_2$, it is necessary to recalculate it
at 298~K to make proper comparisons. The value shown in Table
\ref{tab:liquidProperties} was calculated from the same NMR data in the 
fashion described in Ref. \cite{Krynicki66}. Similarly, $\tau_1$ was
recomputed for 298~K from the data in Ref. \cite{Eisenberg69}.
Again, SSD/E and SSD/RF show improved behavior over SSD1, both with
and without an active reaction field. Turning on the reaction field
leads to much improved time constants for SSD1; however, these results
also include a corresponding decrease in system density.
Orientational relaxation times published in the original SSD dynamics
paper are smaller than the values observed here, and this difference
can be attributed to the use of the Ewald sum.\cite{Chandra99}

\subsection{SSD/RF and Damped Electrostatics}\label{sec:ssdrfDamped}

In section \ref{sec:dampingMultipoles}, a method was described for
applying the damped {\sc sf} or {\sc sp} techniques to for systems
containing point multipoles. The SSD family of water models is the
perfect test case because of the dipole-dipole (and
charge-dipole/quadrupole) interactions that are present. The {\sc sf}
and {\sc sp} techniques were presented as a pairwise replacement for
the Ewald summation. It has been suggested that models parametrized
for the Ewald summation (like TIP5P-E) would be appropriate for use
with a reaction field and vice versa.\cite{Rick04} Therefore, we
decided to test the SSD/RF water model with this damped electrostatic
technique in place of the reaction field to see how the calculated
properties change.

\begin{table}
\caption{PROPERTIES OF SSD/RF WHEN USING DIFFERENT ELECTROSTATIC 
CORRECTION METHODS}
\footnotesize
\centering
\begin{tabular}{ llccc } 
\toprule
\toprule
& & Reaction Field & Damped Electrostatics & Experiment [Ref.] \\
& & $\epsilon = 80$ & $\alpha = 0.2125$~\AA$^{-1}$ & \\
\midrule
$\rho$ & (g cm$^{-3}$) & 0.997(1) & 1.004(4) & 0.997 \cite{CRC80}\\
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 23.8(2) & 27(1) & 18.005 \cite{Wagner02} \\
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.32(6) & 2.33(2) & 2.299 \cite{Mills73}\\
$n_C$ & & 4.4 & 4.2 & 4.7 \cite{Hura00}\\ 
$n_H$ & & 3.7 & 3.7 & 3.5 \cite{Soper86}\\
$\tau_1$ & (ps) & 7.2(4) & 5.86(8) & 5.7 \cite{Eisenberg69}\\
$\tau_2$ & (ps) & 3.2(2) & 2.45(7) & 2.3 \cite{Krynicki66}\\
\bottomrule
\end{tabular}
\label{tab:dampedSSDRF}
\end{table}

The properties shown in table \ref{tab:dampedSSDRF} compare
surprisingly well. The average density shows a modest increase when
using damped electrostatics in place of the reaction field. This comes
about because we neglect the pressure effect due to the surroundings
outside of the cuttoff, instead relying on screening effects to
neutralize electrostatic interactions at long distances. The $C_p$
also shows a slight increase, indicating greater fluctuation in the
enthalpy at constant pressure. The only other differences between the
damped and reaction field results are the dipole reorientational time
constants, $\tau_1$ and $\tau_2$. When using damped electrostatics,
the water molecules relax more quickly and exhibit behavior very
similar to the experimental values. These results indicate that not
only is it reasonable to use damped electrostatics with SSD/RF, it is
recommended if capturing realistic dynamics is of primary
importance. This is an encouraging result because the damping methods
are more generally applicable than reaction field. Using damping,
SSD/RF can be used effectively with mixed systems, such as dissolved
ions, dissolved organic molecules, or even proteins.

In addition to the properties tabulated in table
\ref{tab:dampedSSDRF}, we calculated the static dielectric constant
from a 5~ns simulation of SSD/RF using the damped electrostatics. The
resulting value of 82.6(6) compares very favorably with the
experimental value of 78.3.\cite{Malmberg56} This value is closer to
the experimental value than what was expected according to figure
\ref{fig:dielectricMap}, raising some questions as to the accuracy of
the visual contours in the figure. This highlights the qualitative
nature of contour plotting.

\section{Tetrahedrally Restructured Elongated Dipole (TRED) Water Model}\label{sec:tredWater}

The SSD/RF model works well with damped electrostatics, but because of
its point multipole character, there is no charge neutralization
correction at $R_\textrm{c}$. This has the effect of increasing the
density, since there is no consideration of the ``surroundings''. In
an attempt to optimize a water model for these conditions, we decided
to both simplify the parameters of the SSD type models and break the
point dipole into a charge pair so that it will gain some effect from
the shifting action in the {\sc sf} technique. This process resulted
in a two-point model that we are calling the Tetrahedrally
Restructured Elongated Dipole (TRED) water model.

\begin{figure}
\centering
\includegraphics[width=2.75in]{./figures/tredLayout.pdf}
\caption{Geometry of TRED water. In this two-point model, the origin 
is the center-of-mass of the water molecule with the same moments of
inertia as SSD/RF. The negatively charged site at the origin is also
both a Lennard-Jones and sticky interaction site.}
\label{fig:tredLayout}
\end{figure}
\begin{table}
\centering
\caption{PARAMETERS FOR THE TRED WATER MODEL}
\begin{tabular}{ llr }
\toprule
\toprule
$\sigma$ & (\AA) & 2.980 \\
$\epsilon$ & (kcal mol$^{-1}$) & 0.2045 \\
$q$ & ($e$) & 1.041 \\
$v_0$ & (kcal mol$^{-1}$) & 4.22 \\
$\omega^\circ$ & & 0.07715 \\
$r_l$ \& $r_l^\prime$ & (\AA) & 2.4 \\
$r_u$ \& $r_u^\prime$ & (\AA) & 4.0 \\
Q$_xx$ & (D \AA) & -1.682 \\
Q$_yy$ & (D \AA) & 1.762 \\
Q$_zz$ & (D \AA) & -0.080 \\
\bottomrule
\end{tabular}
\label{tab:tredParams}
\end{table}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/tredGofR.pdf}
\caption{The $g_\textrm{OO}(r)$ for TRED water. The first peak has a 
closer to the experimental plot; however, the second and third peaks
exhibit a more expanded structure similar to SSD/RF. The minimum
following the first peak is located at 3.55~\AA , which is further out
than both experiment and SSD/RF (3.42 and 3.3~\AA\ respectively),
leading to a larger coordination number. If all the curves were
integrated to the experimental minimum, the $n_C$ results would all be
similar, with TRED having an $n_C$ only 0.2 larger than experiment.}
\label{fig:tredGofR}
\end{figure}
The geometric layout of TRED water is shown in figure
\ref{fig:tredLayout} and the electrostatic, Lennard-Jones, and sticky
parameters are displayed in table \ref{tab:tredParams}. The negatively
charged site is located at the center of mass of the molecule and is
also home to the Lennard-Jones and sticky interactions. The charge
separation distance of 0.5~\AA\ along the dipolar ($z$) axis combined
with the charge magnitude ($q$) results in a dipole moment of
2.5~D. This value is simply the value used for SSD/RF (2.48~D) rounded
to the tenths place. We also unified the sticky parameters for the
switching functions on the repulsive and attractive interactions in
the interest of simplicity, and we left the quadrupole moment elements
and $\omega^\circ$ unaltered. It should be noted that additional logic
needs to be included into the electrostatic code when using TRED to
insure that the charges of each water do not interact with the other
water's quadrupole moment. Finally, the strength of the sticky
interaction ($v_0$) was modified to improve the shape of the first
peaks in $g_\textrm{OO}(r)$ and $g_\textrm{OH}(r)$, while the $\sigma$
and $\epsilon$ values were varied to adjust the location of the first
peak of $g_\textrm{OO}(r)$ (see figure \ref{fig:tredGofR}) and the
density. The $\sigma$ and $\epsilon$ optimization was carried out by
separating the Lennard-Jones potential into near entirely repulsive
($A$) and attractive ($C$) parts, such that
\begin{equation}
\sigma = \left(\frac{A}{C}\right)^\frac{1}{6},
\end{equation}
and
\begin{equation}
\epsilon = \frac{C^2}{4A}.
\end{equation}
The location of the peak is adjusted through $A$, while the
interaction strength is modified through $C$. All of the above changes
were made with the goal of capturing the experimental density and
translational diffusion constant at 298~K and 1~atm.

Since TRED is a two-point water model, we expect its computational
efficiency to fall some place in between the one-point and three-point
water models. In comparative simulations, TRED was approximately 85\%
slower than SSD/RF, while SPC/E was 225\% slower than SSD/RF. While
TRED loses some of the performance advantage of SSD, it is still
nearly twice as computationally efficient as SPC/E and TIP3P.

\begin{table}
\caption{PROPERTIES OF TRED COMPARED WITH SSD/RF AND EXPERIMENT}
\footnotesize
\centering
\begin{tabular}{ llccc } 
\toprule
\toprule
& & SSD/RF & TRED & Experiment [Ref.]\\
& & $\alpha = 0.2125$~\AA$^{-1}$ & $\alpha = 0.2125$~\AA$^{-1}$ & \\
\midrule
$\rho$ & (g cm$^{-3}$) & 1.004(4) & 0.995(5) & 0.997 \cite{CRC80}\\
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 27(1) & 23(1) & 18.005 \cite{Wagner02} \\
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.33(2) & 2.30(5) & 2.299 \cite{Mills73}\\
$n_C$ & & 4.2 & 5.3 & 4.7 \cite{Hura00}\\ 
$n_H$ & & 3.7 & 4.1 & 3.5 \cite{Soper86}\\
$\tau_1$ & (ps) & 5.86(8) & 6.0(1) & 5.7 \cite{Eisenberg69}\\
$\tau_2$ & (ps) & 2.45(7) & 2.49(5) & 2.3 \cite{Krynicki66}\\
$\epsilon_0$ & & 82.6(6) & 83(1) & 78.3 \cite{Malmberg56}\\
$\tau_D$ & (ps) & 9.1(2) & 10.6(3) & 8.2(4) \cite{Kindt96}\\
\bottomrule
\end{tabular}
\label{tab:tredProperties}
\end{table}
We calculated thermodynamic and dynamic properties in the same manner
described in section \ref{sec:ssdrfDamped} for SSD/RF, and the results
are presented in table \ref{tab:tredProperties}. These results show that
TRED improves upon the $\rho$ and $C_p$ properties of SSD/RF with
damped electrostatics while maintaining the excellent dynamics
behavior of both the translational self-diffusivity and the
reorientational time constants, $\tau_1$ and $\tau_2$. The structural
results show some differences between the two models. Despite the
improved peak shape for the first solvation shell (see figure
\ref{fig:tredGofR}), $n_C$ and $n_H$ counts increase because of the
further first minimum distance locations. This results in the
integration extending over a larger water volume. If we integrate to
the first minimum value of the experimental $g_\textrm{OO}(r)$
(3.42~\AA ) in both the SSD/RF and TRED plots, the $n_C$ values for
both are much closer to experiment (4.7 for SSD/RF and 4.9 for TRED).

\begin{figure}
\centering
\includegraphics[width=3in]{./figures/tredDielectric.pdf}
\caption{Contour map of the dielectric constant for TRED as a function
of damping parameter and cutoff radius. The dielectric behavior for TRED
is very similar to SSD/RF (see figure \ref{fig:dielectricMap}D), which 
is to be expected due to the similar dipole moment and sticky interaction
strength.}
\label{fig:tredDielectric}
\end{figure}
A comparison of the dielectric behavior was also included at the
bottom of table \ref{tab:tredProperties}. The static dielectric
constant results for SSD/RF and TRED are identical within error. This
is not surprising given the similar dipole moment, similar sticky
interaction strength, and identical applied damping
constant. Comparing the static dielectric constant contour map (figure
\ref{fig:tredDielectric}) with the dielectric map for SSD/RF (figure
\ref{fig:dielectricMap}D) highlights the similarities in how these
models respond to dielectric damping and how the dipolar and monopolar
electrostatic damping act in an equivalent fashion. Both these
dielectric maps span a larger range than the 3, 4, and 5 point-charge
water models; however, the SSD/RF range is greater than TRED,
indicating that multipoles are a little more sensitive to damping than
monopoles. 

The final dielectric comparison comes through the Debye relaxation
time ($\tau_D$) or the collective dipolar relaxation time when
assuming a Debye treatment for the dielectric
relaxation.\cite{Chandra99,Kindt96} This value is calculated through
equation (\ref{eq:reorientCorr}) applied to the total system dipole
moment. The values for both of the models are a slower than the
experimental relaxation; however, they compare compare very well to
experiment considering the Debye relaxation times calculated for the
original SSD (11.95~ps) and the SPC/E (6.95~ps) and TIP3P (6.1~ps)
values. The $\tau_D$ for TRED is about 1.5~ps slower than the $\tau_D$
for SSD/RF, most likely due to the slower decay of the charge-charge
interaction, even when screened by the same damping constant.

\section{Conclusions}

In the above sections, the density maximum and temperature dependence
of the self-diffusion constant were studied for the SSD water model,
both with and without the use of reaction field, via a series of {\it
NPT} and {\it NVE} simulations. The constant pressure simulations
showed a density maximum near 260~K. In most cases, the calculated
densities were significantly lower than the densities obtained from
other water models (and experiment). Analysis of self-diffusion showed
SSD to capture the transport properties of water well in both the
liquid and supercooled liquid regimes.

In order to correct the density behavior, we reparametrized the
original SSD model for use both with and without a reaction field
(SSD/RF and SSD/E), and made comparisons with SSD1, an alternate
density corrected version of SSD. Both models improve the liquid
structure, densities, and diffusive properties under their respective
simulation conditions, indicating the necessity of reparametrization
when changing the method of calculating long-range electrostatic
interactions.  

We also showed that SSD/RF performs well under the alternative damped
electrostatic conditions, validating the multipolar damping work in
the previous chapter. To improve the modeling of water when {\sc sf},
the TRED water model was developed. This model maintains improves upon
the thermodynamic properties of SSD/RF using damped electrostatics
while maintaining the exceptional depiction of water dynamics.

The simple water models investigated here are excellent choices for
representing water in large scale simulations of biochemical
systems. They are more computationally efficient than the common
charge based water models, and often exhibit more realistic bulk phase
fluid properties. These models are one of the few cases in which
maximizing efficiency does not result in a loss in realistic
representation.
Revision:	3023
Committed:	Tue Sep 26 03:07:59 2006 UTC (17 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	61005 byte(s)
Log Message:	Refinements. Just a little proof checking left...