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 |