| 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 |
| 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, an ensemble |
| 321 |
< |
average was accumulated from five separate simulation progressions, |
| 322 |
< |
each starting from a different ice crystal. |
| 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 |
