| 209 |
|
|
| 210 |
|
Long-range dipole-dipole interactions were accounted for in this study |
| 211 |
|
by using either the reaction field method or by resorting to a simple |
| 212 |
< |
cubic switching function at a cutoff radius. Under the first method, |
| 213 |
< |
the magnitude of the reaction field acting on dipole $i$ is |
| 212 |
> |
cubic switching function at a cutoff radius. The reaction field |
| 213 |
> |
method was actually first used in Monte Carlo simulations of liquid |
| 214 |
> |
water.\cite{Barker73} Under this method, the magnitude of the reaction |
| 215 |
> |
field acting on dipole $i$ is |
| 216 |
|
\begin{equation} |
| 217 |
|
\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1} |
| 218 |
|
\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij})\ , |
| 235 |
|
|
| 236 |
|
We have also performed a companion set of simulations {\it without} a |
| 237 |
|
surrounding dielectric (i.e. using a simple cubic switching function |
| 238 |
< |
at the cutoff radius) and as a result we have two reparamaterizations |
| 239 |
< |
of SSD which could be used either with or without the Reaction Field |
| 238 |
> |
at the cutoff radius), and as a result we have two reparamaterizations |
| 239 |
> |
of SSD which could be used either with or without the reaction field |
| 240 |
|
turned on. |
| 241 |
|
|
| 242 |
|
Simulations to obtain the preferred density were performed in the |
| 254 |
|
symplectic splitting method proposed by Dullweber {\it et |
| 255 |
|
al.}\cite{Dullweber1997} Our reason for selecting this integrator |
| 256 |
|
centers on poor energy conservation of rigid body dynamics using |
| 257 |
< |
traditional quaternion integration.\cite{Evans77,Evans77b} While quaternions |
| 258 |
< |
may work well for orientational motion under NVT or NPT integrators, |
| 259 |
< |
our limits on energy drift in the microcanonical ensemble were quite |
| 260 |
< |
strict, and the drift under quaternions was substantially greater than |
| 261 |
< |
in the symplectic splitting method. This steady drift in the total |
| 260 |
< |
energy has also been observed by Kol {\it et al.}\cite{Laird97} |
| 257 |
> |
traditional quaternion integration.\cite{Evans77,Evans77b} In typical |
| 258 |
> |
microcanonical ensemble simulations, the energy drift when using |
| 259 |
> |
quaternions was substantially greater than when using the symplectic |
| 260 |
> |
splitting method (fig. \ref{timestep}). This steady drift in the |
| 261 |
> |
total energy has also been observed by Kol {\it et al.}\cite{Laird97} |
| 262 |
|
|
| 263 |
|
The key difference in the integration method proposed by Dullweber |
| 264 |
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 448 |
|
mean-square displacement as a function of time. The averaged results |
| 449 |
|
from five sets of NVE simulations are displayed in figure |
| 450 |
|
\ref{diffuse}, alongside experimental, SPC/E, and TIP5P |
| 451 |
< |
results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} |
| 451 |
> |
results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} |
| 452 |
|
|
| 453 |
|
\begin{figure} |
| 454 |
|
\begin{center} |
| 456 |
|
\epsfbox{betterDiffuse.epsi} |
| 457 |
|
\caption{Average self-diffusion constant as a function of temperature for |
| 458 |
|
SSD, SPC/E [Ref. \citen{Clancy94}], TIP5P [Ref. \citen{Jorgensen01}], |
| 459 |
< |
and Experimental data [Refs. \citen{Gillen72} and \citen{Mills73}]. Of |
| 459 |
> |
and Experimental data [Refs. \citen{Gillen72} and \citen{Holz00}]. Of |
| 460 |
|
the three water models shown, SSD has the least deviation from the |
| 461 |
|
experimental values. The rapidly increasing diffusion constants for |
| 462 |
|
TIP5P and SSD correspond to significant decrease in density at the |
| 892 |
|
radial location of the minima following the first solvation shell |
| 893 |
|
peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for |
| 894 |
|
calculation of the coordination number or hydrogen bonds per particle |
| 895 |
< |
respectively. |
| 895 |
> |
respectively. The number of hydrogen bonds stays relatively constant |
| 896 |
> |
across all of the models, but the coordination numbers of SSD/E and |
| 897 |
> |
SSD/RF show an improvement over SSD1. This improvement is primarily |
| 898 |
> |
due to the widening of the first solvation shell peak, allowing the |
| 899 |
> |
first minima to push outward. Comparing the coordination number with |
| 900 |
> |
the number of hydrogen bonds can lead to more insight into the |
| 901 |
> |
structural character of the liquid. Because of the near identical |
| 902 |
> |
values for SSD1, it appears to be a little too exclusive, in that all |
| 903 |
> |
molecules in the first solvation shell are involved in forming ideal |
| 904 |
> |
hydrogen bonds. The differing numbers for the newly parameterized |
| 905 |
> |
models indicate the allowance of more fluid configurations in addition |
| 906 |
> |
to the formation of an acceptable number of ideal hydrogen bonds. |
| 907 |
|
|
| 908 |
|
The time constants for the self orientational autocorrelation function |
| 909 |
|
are also displayed in Table \ref{liquidproperties}. The dipolar |
| 918 |
|
vector can be calculated from an exponential fit in the long-time |
| 919 |
|
regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculation of these |
| 920 |
|
time constants were averaged from five detailed NVE simulations |
| 921 |
< |
performed at the STP density for each of the respective models. |
| 921 |
> |
performed at the STP density for each of the respective models. It |
| 922 |
> |
should be noted that the commonly cited value for $\tau_2$ of 1.9 ps |
| 923 |
> |
was determined from the NMR data in reference \citen{Krynicki66} at a |
| 924 |
> |
temperature near 34$^\circ$C.\cite{Rahman73} Because of the strong |
| 925 |
> |
temperature dependence of $\tau_2$, it is necessary to recalculate it |
| 926 |
> |
at 298 K to make proper comparisons. The value shown in Table |
| 927 |
> |
\ref{liquidproperties} was calculated from the same NMR data in the |
| 928 |
> |
fashion described in reference \citen{Krynicki66}. Again, SSD/E and |
| 929 |
> |
SSD/RF show improved behavior over SSD1, both with and without an |
| 930 |
> |
active reaction field. Turning on the reaction field leads to much |
| 931 |
> |
improved time constants for SSD1; however, these results also include |
| 932 |
> |
a corresponding decrease in system density. Numbers published from the |
| 933 |
> |
original SSD dynamics studies appear closer to the experimental |
| 934 |
> |
values, and this difference can be attributed to the use of the Ewald |
| 935 |
> |
sum technique versus a reaction field.\cite{Ichiye99} |
| 936 |
|
|
| 937 |
|
\subsection{Additional Observations} |
| 938 |
|
|
