| 36 |
|
\begin{abstract} |
| 37 |
|
NVE and NPT molecular dynamics simulations were performed in order to |
| 38 |
|
investigate the density maximum and temperature dependent transport |
| 39 |
< |
for the SSD water model, both with and without the use of reaction |
| 40 |
< |
field. The constant pressure simulations of the melting of both $I_h$ |
| 41 |
< |
and $I_c$ ice showed a density maximum near 260 K. In most cases, the |
| 42 |
< |
calculated densities were significantly lower than the densities |
| 43 |
< |
calculated in simulations of other water models. Analysis of particle |
| 44 |
< |
diffusion showed SSD to capture the transport properties of |
| 39 |
> |
for SSD and related water models, both with and without the use of |
| 40 |
> |
reaction field. The constant pressure simulations of the melting of |
| 41 |
> |
both $I_h$ and $I_c$ ice showed a density maximum near 260 K. In most |
| 42 |
> |
cases, the calculated densities were significantly lower than the |
| 43 |
> |
densities calculated in simulations of other water models. Analysis of |
| 44 |
> |
particle diffusion showed SSD to capture the transport properties of |
| 45 |
|
experimental very well in both the normal and super-cooled liquid |
| 46 |
|
regimes. In order to correct the density behavior, SSD was |
| 47 |
|
reparameterized for use both with and without a long-range interaction |
| 146 |
|
simulations using this model, Ichiye \emph{et al.} reported a |
| 147 |
|
calculation speed up of up to an order of magnitude over other |
| 148 |
|
comparable models while maintaining the structural behavior of |
| 149 |
< |
water.\cite{Ichiye96} In the original molecular dynamics studies of |
| 150 |
< |
SSD, it was shown that it actually improves upon the prediction of |
| 151 |
< |
water's dynamical properties 3 and 4-point models.\cite{Ichiye99} This |
| 149 |
> |
water.\cite{Ichiye96} In the original molecular dynamics studies, it |
| 150 |
> |
was shown that SSD improves on the prediction of many of water's |
| 151 |
> |
dynamical properties over TIP3P and SPC/E.\cite{Ichiye99} This |
| 152 |
|
attractive combination of speed and accurate depiction of solvent |
| 153 |
|
properties makes SSD a model of interest for the simulation of large |
| 154 |
|
scale biological systems, such as membrane phase behavior, a specific |
| 205 |
|
to the use of reaction field, simulations were also performed without |
| 206 |
|
a surrounding dielectric and suggestions are proposed on how to make |
| 207 |
|
SSD more compatible with a reaction field. |
| 208 |
< |
|
| 208 |
> |
|
| 209 |
|
Simulations were performed in both the isobaric-isothermal and |
| 210 |
|
microcanonical ensembles. The constant pressure simulations were |
| 211 |
|
implemented using an integral thermostat and barostat as outlined by |
| 212 |
< |
Hoover.\cite{Hoover85,Hoover86} For the constant pressure |
| 213 |
< |
simulations, the \emph{Q} parameter for the was set to 5.0 amu |
| 214 |
< |
\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at |
| 215 |
< |
100 ps. |
| 212 |
> |
Hoover.\cite{Hoover85,Hoover86} All particles were treated as |
| 213 |
> |
non-linear rigid bodies. Vibrational constraints are not necessary in |
| 214 |
> |
simulations of SSD, because there are no explicit hydrogen atoms, and |
| 215 |
> |
thus no molecular vibrational modes need to be considered. |
| 216 |
|
|
| 217 |
|
Integration of the equations of motion was carried out using the |
| 218 |
|
symplectic splitting method proposed by Dullweber \emph{et |
| 220 |
|
deals with poor energy conservation of rigid body systems using |
| 221 |
|
quaternions. While quaternions work well for orientational motion in |
| 222 |
|
alternate ensembles, the microcanonical ensemble has a constant energy |
| 223 |
< |
requirement that is actually quite sensitive to errors in the |
| 224 |
< |
equations of motion. The original implementation of this code utilized |
| 225 |
< |
quaternions for rotational motion propagation; however, a detailed |
| 226 |
< |
investigation showed that they resulted in a steady drift in the total |
| 227 |
< |
energy, something that has been observed by others.\cite{Laird97} |
| 223 |
> |
requirement that is quite sensitive to errors in the equations of |
| 224 |
> |
motion. The original implementation of this code utilized quaternions |
| 225 |
> |
for rotational motion propagation; however, a detailed investigation |
| 226 |
> |
showed that they resulted in a steady drift in the total energy, |
| 227 |
> |
something that has been observed by others.\cite{Laird97} |
| 228 |
|
|
| 229 |
|
The key difference in the integration method proposed by Dullweber |
| 230 |
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 244 |
|
method, the orientational propagation involves a sequence of matrix |
| 245 |
|
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 246 |
|
matrix rotations end up being more costly computationally than the |
| 247 |
< |
simpler arithmetic quaternion propagation. On average, a 1000 SSD |
| 248 |
< |
particle simulation shows a 7\% increase in computation time using the |
| 249 |
< |
symplectic step method in place of quaternions. This cost is more than |
| 250 |
< |
justified when comparing the energy conservation of the two methods as |
| 251 |
< |
illustrated in figure \ref{timestep}. |
| 247 |
> |
simpler arithmetic quaternion propagation. With the same time step, a |
| 248 |
> |
1000 SSD particle simulation shows an average 7\% increase in |
| 249 |
> |
computation time using the symplectic step method in place of |
| 250 |
> |
quaternions. This cost is more than justified when comparing the |
| 251 |
> |
energy conservation of the two methods as illustrated in figure |
| 252 |
> |
\ref{timestep}. |
| 253 |
|
|
| 254 |
|
\begin{figure} |
| 255 |
|
\includegraphics[width=61mm, angle=-90]{timeStep.epsi} |
| 308 |
|
\section{Results and discussion} |
| 309 |
|
|
| 310 |
|
Melting studies were performed on the randomized ice crystals using |
| 311 |
< |
constant pressure and temperature dynamics. This involved an initial |
| 312 |
< |
randomization of velocities about the starting temperature of 25 K for |
| 313 |
< |
varying amounts of time. The systems were all equilibrated for 100 ps |
| 314 |
< |
prior to a 200 ps data collection run at each temperature setting, |
| 315 |
< |
ranging from 25 to 400 K, with a maximum degree increment of 25 K. For |
| 316 |
< |
regions of interest along this stepwise progression, the temperature |
| 317 |
< |
increment was decreased from 25 K to 10 and then 5 K. The above |
| 318 |
< |
equilibration and production times were sufficient in that the system |
| 319 |
< |
volume fluctuations dampened out in all but the very cold simulations |
| 320 |
< |
(below 225 K). In order to further improve statistics, five separate |
| 321 |
< |
simulation progressions were performed, and the averaged results from |
| 322 |
< |
the $I_h$ melting simulations are shown in figure \ref{dense1}. |
| 323 |
< |
|
| 324 |
< |
\begin{figure} |
| 324 |
< |
\includegraphics[width=65mm, angle=-90]{1hdense.epsi} |
| 325 |
< |
\caption{Average density of SSD water at increasing temperatures |
| 326 |
< |
starting from ice $I_h$ lattice.} |
| 327 |
< |
\label{dense1} |
| 328 |
< |
\end{figure} |
| 311 |
> |
constant pressure and temperature dynamics. By performing melting |
| 312 |
> |
simulations, the melting transition can be determined by monitoring |
| 313 |
> |
the heat capacity, in addition to determining the density maximum, |
| 314 |
> |
provided that the density maximum occurs in the liquid and not the |
| 315 |
> |
supercooled regimes. An ensemble average from five separate melting |
| 316 |
> |
simulations was acquired, each starting from different ice crystals |
| 317 |
> |
generated as described previously. All simulations were equilibrated |
| 318 |
> |
for 100 ps prior to a 200 ps data collection run at each temperature |
| 319 |
> |
setting, ranging from 25 to 400 K, with a maximum degree increment of |
| 320 |
> |
25 K. For regions of interest along this stepwise progression, the |
| 321 |
> |
temperature increment was decreased from 25 K to 10 and then 5 K. The |
| 322 |
> |
above equilibration and production times were sufficient in that the |
| 323 |
> |
system volume fluctuations dampened out in all but the very cold |
| 324 |
> |
simulations (below 225 K). |
| 325 |
|
|
| 326 |
|
\subsection{Density Behavior} |
| 327 |
|
In the initial average density versus temperature plot, the density |
| 328 |
< |
maximum clearly appears between 255 and 265 K. The calculated |
| 329 |
< |
densities within this range were nearly indistinguishable, as can be |
| 330 |
< |
seen in the zoom of this region of interest, shown in figure |
| 328 |
> |
maximum appears between 255 and 265 K. The calculated densities within |
| 329 |
> |
this range were nearly indistinguishable, as can be seen in the zoom |
| 330 |
> |
of this region of interest, shown in figure |
| 331 |
|
\ref{dense1}. The greater certainty of the average value at 260 K makes |
| 332 |
|
a good argument for the actual density maximum residing at this |
| 333 |
|
midpoint value. Figure \ref{dense1} was constructed using ice $I_h$ |
| 343 |
|
\begin{figure} |
| 344 |
|
\includegraphics[width=65mm,angle=-90]{dense2.eps} |
| 345 |
|
\caption{Density versus temperature for TIP4P\cite{Jorgensen98b}, |
| 346 |
< |
TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction |
| 347 |
< |
Field, SSD, and Experiment\cite{CRC80}. } |
| 346 |
> |
TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction |
| 347 |
> |
Field, SSD, and Experiment\cite{CRC80}. The arrows indicate the |
| 348 |
> |
change in densities observed when turning off the reaction field. The |
| 349 |
> |
the lower than expected densities for the SSD model were what |
| 350 |
> |
prompted the original reparameterization.\cite{Ichiye03}} |
| 351 |
|
\label{dense2} |
| 352 |
|
\end{figure} |
| 353 |
|
|
| 575 |
|
|
| 576 |
|
\begin{table} |
| 577 |
|
\caption{Parameters for the original and adjusted models} |
| 578 |
< |
\begin{tabular}{ l c c c } |
| 578 |
> |
\begin{tabular}{ l c c c c } |
| 579 |
|
\hline \\[-3mm] |
| 580 |
< |
\ Parameters & \ \ \ SSD$^\dagger$\ \ \ \ & \ SSD/E\ \ & \ SSD/RF\ \ \\ |
| 580 |
> |
\ \ \ Parameters\ \ \ & \ \ \ SSD$^\dagger$ \ \ \ & \ SSD1$^\ddagger$\ \ & \ SSD/E\ \ & \ SSD/RF \\ |
| 581 |
|
\hline \\[-3mm] |
| 582 |
< |
\ \ \ $\sigma$ (\AA) & 3.051 & 3.035 & 3.019\\ |
| 583 |
< |
\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\ |
| 584 |
< |
\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\ |
| 585 |
< |
\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\ |
| 586 |
< |
\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\ |
| 587 |
< |
\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\ |
| 588 |
< |
\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\ |
| 589 |
< |
\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\ |
| 590 |
< |
\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\ |
| 582 |
> |
\ \ \ $\sigma$ (\AA) & 3.051 & 3.016 & 3.035 & 3.019\\ |
| 583 |
> |
\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ |
| 584 |
> |
\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\ |
| 585 |
> |
\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
| 586 |
> |
\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ |
| 587 |
> |
\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\ |
| 588 |
> |
\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
| 589 |
> |
\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ |
| 590 |
> |
\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ |
| 591 |
|
\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96} |
| 592 |
+ |
\\$^\ddagger$ ref. \onlinecite{Ichiye03} |
| 593 |
|
\end{tabular} |
| 594 |
|
\label{params} |
| 595 |
|
\end{table} |
| 596 |
|
|
| 597 |
|
\begin{figure} |
| 598 |
< |
\includegraphics[width=85mm]{gofrCompare.epsi} |
| 598 |
> |
\includegraphics[width=85mm]{GofRCompare.epsi} |
| 599 |
|
\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E |
| 600 |
< |
and SSD without reaction field (top), as well as SSD/RF and SSD with |
| 600 |
> |
and SSD1 without reaction field (top), as well as SSD/RF and SSD1 with |
| 601 |
|
reaction field turned on (bottom). The insets show the respective |
| 602 |
|
first peaks in detail. Solid Line - experiment, dashed line - SSD/E |
| 603 |
< |
and SSD/RF, and dotted line - SSD (with and without reaction field).} |
| 603 |
> |
and SSD/RF, and dotted line - SSD1 (with and without reaction field).} |
| 604 |
|
\label{grcompare} |
| 605 |
|
\end{figure} |
| 606 |
|
|
| 607 |
|
\begin{figure} |
| 608 |
|
\includegraphics[width=85mm]{dualsticky.ps} |
| 609 |
< |
\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \& |
| 609 |
> |
\caption{Isosurfaces of the sticky potential for SSD1 (left) and SSD/E \& |
| 610 |
|
SSD/RF (right). Light areas correspond to the tetrahedral attractive |
| 611 |
|
part, and the darker areas correspond to the dipolar repulsive part.} |
| 612 |
|
\label{isosurface} |
| 683 |
|
stated earlier in this paper. |
| 684 |
|
|
| 685 |
|
\begin{figure} |
| 686 |
< |
\includegraphics[width=85mm]{ssdecompare.epsi} |
| 686 |
> |
\includegraphics[width=62mm, angle=-90]{ssdeDense.epsi} |
| 687 |
|
\caption{Comparison of densities calculated with SSD/E to SSD without a |
| 688 |
< |
reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b}, |
| 689 |
< |
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot |
| 690 |
< |
includes error bars, and the calculated results from the other |
| 691 |
< |
references were removed for clarity.} |
| 688 |
> |
reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00}, |
| 689 |
> |
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The window shows a |
| 690 |
> |
expansion around 300 K with error bars included to clarify this region |
| 691 |
> |
of interest. Note that both SSD1 and SSD/E show good agreement with |
| 692 |
> |
experiment when the long-range correction is neglected.} |
| 693 |
|
\label{ssdedense} |
| 694 |
|
\end{figure} |
| 695 |
|
|
| 715 |
|
justify the modifications taken. |
| 716 |
|
|
| 717 |
|
\begin{figure} |
| 718 |
< |
\includegraphics[width=85mm]{ssdrfcompare.epsi} |
| 718 |
> |
\includegraphics[width=62mm, angle=-90]{ssdrfDense.epsi} |
| 719 |
|
\caption{Comparison of densities calculated with SSD/RF to SSD with a |
| 720 |
< |
reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b}, |
| 721 |
< |
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot |
| 722 |
< |
includes error bars, and the calculated results from the other |
| 723 |
< |
references were removed for clarity.} |
| 720 |
> |
reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00}, |
| 721 |
> |
SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The inset shows the |
| 722 |
> |
necessity of reparameterization when utilizing a reaction field |
| 723 |
> |
long-ranged correction - SSD/RF provides significantly more accurate |
| 724 |
> |
densities than SSD1 when performing room temperature simulations.} |
| 725 |
|
\label{ssdrfdense} |
| 726 |
|
\end{figure} |
| 727 |
|
|
| 759 |
|
lower densities with increasing temperature as rapidly. |
| 760 |
|
|
| 761 |
|
\begin{figure} |
| 762 |
< |
\includegraphics[width=85mm]{ssdediffuse.epsi} |
| 763 |
< |
\caption{Plots of the diffusion constants calculated from SSD/E and SSD, |
| 764 |
< |
both without a reaction field along with experimental results from |
| 765 |
< |
Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The |
| 766 |
< |
upper plot is at densities calculated from the NPT simulations at a |
| 767 |
< |
pressure of 1 atm, while the lower plot is at the experimentally |
| 768 |
< |
calculated densities.} |
| 769 |
< |
\label{ssdediffuse} |
| 762 |
> |
\includegraphics[width=65mm, angle=-90]{ssdrfDiffuse.epsi} |
| 763 |
> |
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1, |
| 764 |
> |
both with an active reaction field, along with experimental results |
| 765 |
> |
from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The |
| 766 |
> |
NVE calculations were performed at the average densities observed in |
| 767 |
> |
the 1 atm NPT simulations for both of the models. Note how accurately |
| 768 |
> |
SSD/RF simulates the diffusion of water throughout this temperature |
| 769 |
> |
range. The more rapidly increasing diffusion constants at high |
| 770 |
> |
temperatures for both models is attributed to the significantly lower |
| 771 |
> |
densities than observed in experiment.} |
| 772 |
> |
\label{ssdrfdiffuse} |
| 773 |
|
\end{figure} |
| 774 |
|
|
| 775 |
|
\begin{figure} |
| 776 |
< |
\includegraphics[width=85mm]{ssdrfdiffuse.epsi} |
| 777 |
< |
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD, |
| 778 |
< |
both with an active reaction field along with experimental results |
| 776 |
> |
\includegraphics[width=65mm, angle=-90]{ssdeDiffuse.epsi} |
| 777 |
> |
\caption{Plots of the diffusion constants calculated from SSD/E and SSD1, |
| 778 |
> |
both without a reaction field, along with experimental results are |
| 779 |
|
from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The |
| 780 |
< |
upper plot is at densities calculated from the NPT simulations at a |
| 781 |
< |
pressure of 1 atm, while the lower plot is at the experimentally |
| 782 |
< |
calculated densities.} |
| 783 |
< |
\label{ssdrfdiffuse} |
| 780 |
> |
NVE calculations were performed at the average densities observed in |
| 781 |
> |
the 1 atm NPT simulations for the respective models. SSD/E is |
| 782 |
> |
slightly more fluid than experiment at all of the temperatures, but |
| 783 |
> |
it is closer than SSD1 without a long-range correction.} |
| 784 |
> |
\label{ssdediffuse} |
| 785 |
|
\end{figure} |
| 786 |
|
|
| 787 |
|
In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are |
| 908 |
|
simulations of biochemical systems. |
| 909 |
|
|
| 910 |
|
\section{Acknowledgments} |
| 911 |
< |
The authors would like to thank the National Science Foundation for |
| 912 |
< |
funding under grant CHE-0134881. Computation time was provided by the |
| 913 |
< |
Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR |
| 914 |
< |
00 79647. |
| 911 |
> |
Support for this project was provided by the National Science |
| 912 |
> |
Foundation under grant CHE-0134881. Computation time was provided by |
| 913 |
> |
the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant |
| 914 |
> |
DMR 00 79647. |
| 915 |
|
|
| 916 |
|
\bibliographystyle{jcp} |
| 917 |
|
|
