| 19 |
|
\label{fig:waterModels} |
| 20 |
|
\end{figure} |
| 21 |
|
|
| 22 |
< |
As discussed in the previous chapter, water it typically modeled with |
| 22 |
> |
As discussed in the previous chapter, water is typically modeled with |
| 23 |
|
fixed geometries of point charges shielded by the repulsive part of a |
| 24 |
|
Lennard-Jones interaction. Some of the common water models are shown |
| 25 |
|
in figure \ref{fig:waterModels}. The various models all have their |
| 31 |
|
computational cost.\cite{Mahoney00,Mahoney01} This cost is entirely |
| 32 |
|
due to the additional distance and electrostatic calculations that |
| 33 |
|
come from the extra point charges in the model description. Thus, the |
| 34 |
< |
criteria for choosing a water model are capturing the liquid state |
| 35 |
< |
properties and having the fewest number of points to insure efficient |
| 36 |
< |
performance. As researchers have begun to study larger systems, such |
| 37 |
< |
as entire viruses, the choice readily shifts towards efficiency over |
| 38 |
< |
accuracy in order to make the calculations |
| 39 |
< |
feasible.\cite{Freddolino06} |
| 34 |
> |
main criteria for choosing a water model are: |
| 35 |
|
|
| 36 |
+ |
\begin{enumerate} |
| 37 |
+ |
\item capturing the liquid state properties and |
| 38 |
+ |
\item having the fewest number of points to insure efficient performance. |
| 39 |
+ |
\end{enumerate} |
| 40 |
+ |
As researchers have begun to study larger systems, such as entire |
| 41 |
+ |
viruses, the choice has shifted towards efficiency over accuracy in |
| 42 |
+ |
order to make the calculations feasible.\cite{Freddolino06} |
| 43 |
+ |
|
| 44 |
|
\section{Soft Sticky Dipole Model for Water} |
| 45 |
|
|
| 46 |
< |
One recently developed model that largely succeeds in retaining the |
| 46 |
> |
One recently-developed model that largely succeeds in retaining the |
| 47 |
|
accuracy of bulk properties while greatly reducing the computational |
| 48 |
|
cost is the Soft Sticky Dipole (SSD) water |
| 49 |
|
model.\cite{Liu96,Liu96b,Chandra99,Tan03} The SSD model was developed |
| 138 |
|
properties and behavior under the more computationally efficient |
| 139 |
|
reaction field (RF) technique, the correction techniques discussed in |
| 140 |
|
the previous chapter, or even a simple |
| 141 |
< |
cutoff.\cite{Onsager36,Fennell06} This study addresses these issues by |
| 142 |
< |
looking at the structural and transport behavior of SSD over a variety |
| 143 |
< |
of temperatures with the purpose of utilizing the RF correction |
| 144 |
< |
technique. We then suggest modifications to the parameters that |
| 145 |
< |
result in more realistic bulk phase behavior. It should be noted that |
| 146 |
< |
in a recent publication, some of the original investigators of the SSD |
| 147 |
< |
water model have suggested adjustments to the SSD water model to |
| 148 |
< |
address abnormal density behavior (also observed here), calling the |
| 141 |
> |
cutoff.\cite{Onsager36,Fennell06} This chapter addresses these issues |
| 142 |
> |
by looking at the structural and transport behavior of SSD over a |
| 143 |
> |
variety of temperatures with the purpose of utilizing the RF |
| 144 |
> |
correction technique. We then suggest modifications to the parameters |
| 145 |
> |
that result in more realistic bulk phase behavior. It should be noted |
| 146 |
> |
that in a recent publication, some of the original investigators of |
| 147 |
> |
the SSD water model have suggested adjustments to the SSD water model |
| 148 |
> |
to address abnormal density behavior (also observed here), calling the |
| 149 |
|
corrected model SSD1.\cite{Tan03} In the later sections of this |
| 150 |
|
chapter, we compare our modified variants of SSD with both the |
| 151 |
|
original SSD and SSD1 models and discuss how our changes improve the |
| 154 |
|
\section{Simulation Methods}{\label{sec:waterSimMethods} |
| 155 |
|
|
| 156 |
|
Most of the calculations in this particular study were performed using |
| 157 |
< |
a internally developed simulation code that was one of the precursors |
| 157 |
> |
an internally developed simulation code that was one of the precursors |
| 158 |
|
of the {\sc oopse} molecular dynamics (MD) package.\cite{Meineke05} |
| 159 |
< |
All of the capabilities of this code have been efficiently |
| 160 |
< |
incorporated into {\sc oopse}, and calculation results are consistent |
| 161 |
< |
between the two simulation packages. The later calculations involving |
| 162 |
< |
the damped shifted force ({\sc sf}) techniques were performed using |
| 160 |
< |
{\sc oopse}. |
| 159 |
> |
All of the capabilities of this code have been incorporated into {\sc |
| 160 |
> |
oopse}, and calculation results are consistent between the two |
| 161 |
> |
simulation packages. The later calculations involving the damped |
| 162 |
> |
shifted force ({\sc sf}) techniques were performed using {\sc oopse}. |
| 163 |
|
|
| 164 |
|
In the primary simulations of this study, long-range dipole-dipole |
| 165 |
< |
interaction corrections were accounted for by using either the |
| 166 |
< |
reaction field technique or a simple cubic switching function at the |
| 167 |
< |
cutoff radius. Interestingly, one of the early applications of a |
| 168 |
< |
reaction field was in Monte Carlo simulations of liquid |
| 169 |
< |
water.\cite{Barker73} In this method, the magnitude of the reaction |
| 170 |
< |
field acting on dipole $i$ is |
| 165 |
> |
corrections were accounted for by using either the reaction field |
| 166 |
> |
technique or a simple cubic switching function at the cutoff |
| 167 |
> |
radius. Interestingly, one of the early applications of a reaction |
| 168 |
> |
field was in Monte Carlo simulations of liquid water.\cite{Barker73} |
| 169 |
> |
In this method, the magnitude of the reaction field acting on dipole |
| 170 |
> |
$i$ is |
| 171 |
|
\begin{equation} |
| 172 |
|
\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1} |
| 173 |
|
\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}), |
| 177 |
|
($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the |
| 178 |
|
system, ${\bf\mu}_{j}$ is the dipole moment vector of particle $j$, |
| 179 |
|
and $s(r_{ij})$ is a cubic switching function.\cite{Allen87} The |
| 180 |
< |
reaction field contribution to the total energy by particle $i$ is |
| 180 |
> |
reaction field contribution to the total energy from particle $i$ is |
| 181 |
|
given by $-\frac{1}{2}{\bf\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque |
| 182 |
|
on dipole $i$ by ${\bf\mu}_{i}\times\mathcal{E}_{i}$.\cite{Allen87} An |
| 183 |
|
applied reaction field will alter the bulk orientational properties of |
| 184 |
|
simulated water, and there is particular sensitivity of these |
| 185 |
|
properties on changes in the length of the cutoff |
| 186 |
< |
radius.\cite{vanderSpoel98} This variable behavior makes reaction |
| 187 |
< |
field a less attractive method than the Ewald sum; however, for very |
| 188 |
< |
large systems, the computational benefit of reaction field is is |
| 187 |
< |
significant. |
| 186 |
> |
radius.\cite{vanderSpoel98} This behavior makes the reaction field a |
| 187 |
> |
less attractive method than the Ewald sum; however, for very large |
| 188 |
> |
systems, the computational benefit of reaction field is significant. |
| 189 |
|
|
| 190 |
|
In contrast to the simulations with a reaction field, we have also |
| 191 |
|
performed a companion set of simulations {\it without} a surrounding |
| 269 |
|
randomly sampled at 400~K. The rotational temperature was then scaled |
| 270 |
|
down in stages to slowly cool the crystals to 25~K. The particles were |
| 271 |
|
then allowed to translate with fixed orientations at a constant |
| 272 |
< |
pressure of 1atm for 50~ps at 25~K. Finally, all constraints were |
| 272 |
> |
pressure of 1~atm for 50~ps at 25~K. Finally, all constraints were |
| 273 |
|
removed and the ice crystals were allowed to equilibrate for 50~ps at |
| 274 |
< |
25~K and a constant pressure of 1atm. This procedure resulted in |
| 274 |
> |
25~K and a constant pressure of 1~atm. This procedure resulted in |
| 275 |
|
structurally stable I$_\textrm{h}$ ice crystals that obey the |
| 276 |
|
Bernal-Fowler rules.\cite{Bernal33,Rahman72} This method was also |
| 277 |
|
utilized in the making of diamond lattice I$_\textrm{c}$ ice crystals, |
| 330 |
|
The density maximum for SSD compares quite favorably to other simple |
| 331 |
|
water models. Figure \ref{fig:ssdDense} also shows calculated |
| 332 |
|
densities of several other models and experiment obtained from other |
| 333 |
< |
sources.\cite{Jorgensen98b,Baez94,CRC80} Of the listed simple water |
| 334 |
< |
models, SSD has a temperature closest to the experimentally observed |
| 335 |
< |
density maximum. Of the {\it charge-based} models in figure |
| 333 |
> |
sources.\cite{Jorgensen98b,Baez94,CRC80} Of the water models shown, |
| 334 |
> |
SSD has a temperature closest to the experimentally observed density |
| 335 |
> |
maximum. Of the {\it charge-based} models in figure |
| 336 |
|
\ref{fig:ssdDense}, TIP4P has a density maximum behavior most like |
| 337 |
|
that seen in SSD. Though not included in this plot, it is useful to |
| 338 |
|
note that TIP5P has a density maximum nearly identical to the |
| 351 |
|
The key feature to recognize in figure \ref{fig:ssdDense} is the |
| 352 |
|
density scaling of SSD relative to other common models at any given |
| 353 |
|
temperature. SSD assumes a lower density than any of the other listed |
| 354 |
< |
models at the same pressure, behavior which is especially apparent at |
| 354 |
> |
models at the same pressure. This behavior is especially apparent at |
| 355 |
|
temperatures greater than 300~K. Lower than expected densities have |
| 356 |
|
been observed for other systems using a reaction field for long-range |
| 357 |
|
electrostatic interactions, so the most likely reason for the reduced |
| 460 |
|
\centering |
| 461 |
|
\includegraphics[width=2.5in]{./figures/corrDiag.pdf} |
| 462 |
|
\caption{ An illustration of angles involved in the correlations |
| 463 |
< |
observed in figure \ref{fig:contour}.} |
| 463 |
> |
displayed in figure \ref{fig:contour}.} |
| 464 |
|
\label{fig:corrAngle} |
| 465 |
|
\end{figure} |
| 466 |
|
|
| 715 |
|
\label{fig:ssdrfDense} |
| 716 |
|
\end{figure} |
| 717 |
|
|
| 718 |
< |
Including the reaction field long-range correction results in a more |
| 718 |
> |
Including the reaction field long-range correction results is a more |
| 719 |
|
interesting comparison. A density profile including SSD/RF and SSD1 |
| 720 |
|
with an active reaction field is shown in figure \ref{fig:ssdrfDense}. |
| 721 |
|
As observed in the simulations without a reaction field, the densities |
| 932 |
|
enthalpy at constant pressure. The only other differences between the |
| 933 |
|
damped and reaction field results are the dipole reorientational time |
| 934 |
|
constants, $\tau_1$ and $\tau_2$. When using damped electrostatics, |
| 935 |
< |
the water molecules relax more quickly and are almost identical to the |
| 936 |
< |
experimental values. These results indicate that not only is it |
| 937 |
< |
reasonable to use damped electrostatics with SSD/RF, it is recommended |
| 938 |
< |
if capturing realistic dynamics is of primary importance. This is an |
| 939 |
< |
encouraging result because of the more varied applicability of damping |
| 940 |
< |
over the reaction field technique. Rather than be limited to |
| 941 |
< |
homogeneous systems, SSD/RF can be used effectively with mixed |
| 942 |
< |
systems, such as dissolved ions, dissolved organic molecules, or even |
| 942 |
< |
proteins. |
| 935 |
> |
the water molecules relax more quickly and exhibit behavior very |
| 936 |
> |
similar to the experimental values. These results indicate that not |
| 937 |
> |
only is it reasonable to use damped electrostatics with SSD/RF, it is |
| 938 |
> |
recommended if capturing realistic dynamics is of primary |
| 939 |
> |
importance. This is an encouraging result because the damping methods |
| 940 |
> |
are more generally applicable than reaction field. Using damping, |
| 941 |
> |
SSD/RF can be used effectively with mixed systems, such as dissolved |
| 942 |
> |
ions, dissolved organic molecules, or even proteins. |
| 943 |
|
|
| 944 |
|
In addition to the properties tabulated in table |
| 945 |
|
\ref{tab:dampedSSDRF}, we calculated the static dielectric constant |
| 1025 |
|
peaks in $g_\textrm{OO}(r)$ and $g_\textrm{OH}(r)$, while the $\sigma$ |
| 1026 |
|
and $\epsilon$ values were varied to adjust the location of the first |
| 1027 |
|
peak of $g_\textrm{OO}(r)$ (see figure \ref{fig:tredGofR}) and the |
| 1028 |
< |
density. The $\sigma$ and $epsilon$ optimization was carried out by |
| 1028 |
> |
density. The $\sigma$ and $\epsilon$ optimization was carried out by |
| 1029 |
|
separating the Lennard-Jones potential into near entirely repulsive |
| 1030 |
|
($A$) and attractive ($C$) parts, such that |
| 1031 |
|
\begin{equation} |
| 1040 |
|
were made with the goal of capturing the experimental density and |
| 1041 |
|
translational diffusion constant at 298~K and 1~atm. |
| 1042 |
|
|
| 1043 |
< |
Being that TRED is a two-point water model, we expect its |
| 1044 |
< |
computational efficiency to fall some place in between the single and |
| 1045 |
< |
three-point water models. In comparative simulations, TRED was |
| 1046 |
< |
approximately 85\% slower than SSD/RF, while SPC/E was 225\% slower |
| 1047 |
< |
than SSD/RF. While TRED loses some of the performance advantage of |
| 1048 |
< |
SSD, it is still nearly twice as computationally efficient as SPC/E |
| 1049 |
< |
and TIP3P. |
| 1043 |
> |
Since TRED is a two-point water model, we expect its computational |
| 1044 |
> |
efficiency to fall some place in between the one-point and three-point |
| 1045 |
> |
water models. In comparative simulations, TRED was approximately 85\% |
| 1046 |
> |
slower than SSD/RF, while SPC/E was 225\% slower than SSD/RF. While |
| 1047 |
> |
TRED loses some of the performance advantage of SSD, it is still |
| 1048 |
> |
nearly twice as computationally efficient as SPC/E and TIP3P. |
| 1049 |
|
|
| 1050 |
|
\begin{table} |
| 1051 |
|
\caption{PROPERTIES OF TRED COMPARED WITH SSD/RF AND EXPERIMENT} |
| 1153 |
|
while maintaining the exceptional depiction of water dynamics. |
| 1154 |
|
|
| 1155 |
|
The simple water models investigated here are excellent choices for |
| 1156 |
< |
representing explicit water in large scale simulations of biochemical |
| 1156 |
> |
representing water in large scale simulations of biochemical |
| 1157 |
|
systems. They are more computationally efficient than the common |
| 1158 |
< |
charge based water models, and, in many cases, exhibit more realistic |
| 1159 |
< |
bulk phase fluid properties. These models are one of the few cases in |
| 1160 |
< |
which maximizing efficiency does not result in a loss in realistic |
| 1161 |
< |
liquid water representation. |
| 1158 |
> |
charge based water models, and often exhibit more realistic bulk phase |
| 1159 |
> |
fluid properties. These models are one of the few cases in which |
| 1160 |
> |
maximizing efficiency does not result in a loss in realistic |
| 1161 |
> |
representation. |
