139 |
|
\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) |
140 |
|
+ s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf |
141 |
|
\Omega}_j)]\ . |
142 |
+ |
\label{stickyfunction} |
143 |
|
\end{equation} |
144 |
|
Here, $\nu_0$ is a strength parameter for the sticky potential, and |
145 |
|
$s$ and $s^\prime$ are cubic switching functions which turn off the |
152 |
|
while the $w^\prime$ function counters the normal aligned and |
153 |
|
anti-aligned structures favored by point dipoles: |
154 |
|
\begin{equation} |
155 |
< |
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0, |
155 |
> |
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^\circ, |
156 |
|
\end{equation} |
157 |
|
It should be noted that $w$ is proportional to the sum of the $Y_3^2$ |
158 |
|
and $Y_3^{-2}$ spherical harmonics (a linear combination which |
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})\ , |
218 |
> |
\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij}), |
219 |
|
\label{rfequation} |
220 |
|
\end{equation} |
221 |
|
where $\mathcal{R}$ is the cavity defined by the cutoff radius |
222 |
|
($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the |
223 |
|
system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole |
224 |
< |
moment vector of particle $j$ and $f(r_{ij})$ is a cubic switching |
224 |
> |
moment vector of particle $j$, and $f(r_{ij})$ is a cubic switching |
225 |
|
function.\cite{AllenTildesley} The reaction field contribution to the |
226 |
|
total energy by particle $i$ is given by $-\frac{1}{2}{\bf |
227 |
|
\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf |
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 |
251 |
|
need to be considered. |
252 |
|
|
253 |
|
Integration of the equations of motion was carried out using 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 |
259 |
< |
energy has also been observed by Kol {\it et al.}\cite{Laird97} |
254 |
> |
symplectic splitting method proposed by Dullweber, Leimkuhler, and |
255 |
> |
McLachlan (DLM).\cite{Dullweber1997} Our reason for selecting this |
256 |
> |
integrator centers on poor energy conservation of rigid body dynamics |
257 |
> |
using traditional quaternion integration.\cite{Evans77,Evans77b} In |
258 |
> |
typical microcanonical ensemble simulations, the energy drift when |
259 |
> |
using quaternions was substantially greater than when using the DLM |
260 |
> |
method (fig. \ref{timestep}). This steady drift in the total energy |
261 |
> |
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 |
267 |
|
rotation matrix into quaternions for storage purposes makes trajectory |
268 |
|
data quite compact. |
269 |
|
|
270 |
< |
The symplectic splitting method allows for Verlet style integration of |
271 |
< |
both translational and orientational motion of rigid bodies. In this |
270 |
> |
The DML method allows for Verlet style integration of both |
271 |
> |
translational and orientational motion of rigid bodies. In this |
272 |
|
integration method, the orientational propagation involves a sequence |
273 |
|
of matrix evaluations to update the rotation |
274 |
|
matrix.\cite{Dullweber1997} These matrix rotations are more costly |
275 |
|
than the simpler arithmetic quaternion propagation. With the same time |
276 |
|
step, a 1000 SSD particle simulation shows an average 7\% increase in |
277 |
< |
computation time using the symplectic step method in place of |
278 |
< |
quaternions. The additional expense per step is justified when one |
279 |
< |
considers the ability to use time steps that are nearly twice as large |
280 |
< |
under symplectic splitting than would be usable under quaternion |
281 |
< |
dynamics. The energy conservation of the two methods using a number |
282 |
< |
of different time steps is illustrated in figure |
277 |
> |
computation time using the DML method in place of quaternions. The |
278 |
> |
additional expense per step is justified when one considers the |
279 |
> |
ability to use time steps that are nearly twice as large under DML |
280 |
> |
than would be usable under quaternion dynamics. The energy |
281 |
> |
conservation of the two methods using a number of different time steps |
282 |
> |
is illustrated in figure |
283 |
|
\ref{timestep}. |
284 |
|
|
285 |
|
\begin{figure} |
287 |
|
\epsfxsize=6in |
288 |
|
\epsfbox{timeStep.epsi} |
289 |
|
\caption{Energy conservation using both quaternion based integration and |
290 |
< |
the symplectic step method proposed by Dullweber \emph{et al.} with |
291 |
< |
increasing time step. The larger time step plots are shifted from the |
292 |
< |
true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.} |
290 |
> |
the symplectic splitting method proposed by Dullweber \emph{et al.} |
291 |
> |
with increasing time step. The larger time step plots are shifted from |
292 |
> |
the true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.} |
293 |
|
\label{timestep} |
294 |
|
\end{center} |
295 |
|
\end{figure} |
296 |
|
|
297 |
|
In figure \ref{timestep}, the resulting energy drift at various time |
298 |
< |
steps for both the symplectic step and quaternion integration schemes |
299 |
< |
is compared. All of the 1000 SSD particle simulations started with |
300 |
< |
the same configuration, and the only difference was the method used to |
301 |
< |
handle orientational motion. At time steps of 0.1 and 0.5 fs, both |
302 |
< |
methods for propagating the orientational degrees of freedom conserve |
303 |
< |
energy fairly well, with the quaternion method showing a slight energy |
304 |
< |
drift over time in the 0.5 fs time step simulation. At time steps of 1 |
305 |
< |
and 2 fs, the energy conservation benefits of the symplectic step |
306 |
< |
method are clearly demonstrated. Thus, while maintaining the same |
307 |
< |
degree of energy conservation, one can take considerably longer time |
308 |
< |
steps, leading to an overall reduction in computation time. |
298 |
> |
steps for both the DML and quaternion integration schemes is compared. |
299 |
> |
All of the 1000 SSD particle simulations started with the same |
300 |
> |
configuration, and the only difference was the method used to handle |
301 |
> |
orientational motion. At time steps of 0.1 and 0.5 fs, both methods |
302 |
> |
for propagating the orientational degrees of freedom conserve energy |
303 |
> |
fairly well, with the quaternion method showing a slight energy drift |
304 |
> |
over time in the 0.5 fs time step simulation. At time steps of 1 and 2 |
305 |
> |
fs, the energy conservation benefits of the DML method are clearly |
306 |
> |
demonstrated. Thus, while maintaining the same degree of energy |
307 |
> |
conservation, one can take considerably longer time steps, leading to |
308 |
> |
an overall reduction in computation time. |
309 |
|
|
310 |
< |
Energy drift in the symplectic step simulations was unnoticeable for |
311 |
< |
time steps up to 3 fs. A slight energy drift on the |
312 |
< |
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
313 |
< |
4 fs, and as expected, this drift increases dramatically |
314 |
< |
with increasing time step. To insure accuracy in our microcanonical |
315 |
< |
simulations, time steps were set at 2 fs and kept at this value for |
316 |
< |
constant pressure simulations as well. |
310 |
> |
Energy drift in the simulations using DML integration was unnoticeable |
311 |
> |
for time steps up to 3 fs. A slight energy drift on the order of 0.012 |
312 |
> |
kcal/mol per nanosecond was observed at a time step of 4 fs, and as |
313 |
> |
expected, this drift increases dramatically with increasing time |
314 |
> |
step. To insure accuracy in our microcanonical simulations, time steps |
315 |
> |
were set at 2 fs and kept at this value for constant pressure |
316 |
> |
simulations as well. |
317 |
|
|
318 |
|
Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices |
319 |
|
were generated as starting points for all simulations. The $I_h$ |
430 |
|
simulations.\cite{Clancy94,Jorgensen98b} However, even without the |
431 |
|
reaction field, the density around 300 K is still significantly lower |
432 |
|
than experiment and comparable water models. This anomalous behavior |
433 |
< |
was what lead Ichiye {\it et al.} to recently reparameterize |
433 |
> |
was what lead Tan {\it et al.} to recently reparameterize |
434 |
|
SSD.\cite{Ichiye03} Throughout the remainder of the paper our |
435 |
|
reparamaterizations of SSD will be compared with the newer SSD1 model. |
436 |
|
|
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 |
593 |
|
important properties. In this case, it would be ideal to correct the |
594 |
|
densities while maintaining the accurate transport behavior. |
595 |
|
|
596 |
< |
The parameters available for tuning include the $\sigma$ and $\epsilon$ |
597 |
< |
Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky |
598 |
< |
attractive and dipole repulsive terms with their respective |
599 |
< |
cutoffs. To alter the attractive and repulsive terms of the sticky |
600 |
< |
potential independently, it is necessary to separate the terms as |
601 |
< |
follows: |
602 |
< |
\begin{equation} |
601 |
< |
u_{ij}^{sp} |
602 |
< |
({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = |
603 |
< |
\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)], |
604 |
< |
\end{equation} |
605 |
< |
where $\nu_0$ scales the strength of the tetrahedral attraction and |
606 |
< |
$\nu_0^\prime$ scales the dipole repulsion term independently. The |
607 |
< |
separation was performed for purposes of the reparameterization, but |
608 |
< |
the final parameters were adjusted so that it is not necessary to |
609 |
< |
separate the terms when implementing the adjusted water |
610 |
< |
potentials. The results of the reparameterizations are shown in table |
611 |
< |
\ref{params}. Note that the tetrahedral attractive and dipolar |
612 |
< |
repulsive terms do not share the same lower cutoff ($r_l$) in the |
613 |
< |
newly parameterized potentials. We are calling these |
614 |
< |
reparameterizations the Soft Sticky Dipole / Reaction Field |
596 |
> |
The parameters available for tuning include the $\sigma$ and |
597 |
> |
$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the |
598 |
> |
strength of the sticky potential ($\nu_0$), and the sticky attractive |
599 |
> |
and dipole repulsive cubic switching function cutoffs ($r_l$, $r_u$ |
600 |
> |
and $r_l^\prime$, $r_u^\prime$ respectively). The results of the |
601 |
> |
reparameterizations are shown in table \ref{params}. We are calling |
602 |
> |
these reparameterizations the Soft Sticky Dipole / Reaction Field |
603 |
|
(SSD/RF - for use with a reaction field) and Soft Sticky Dipole |
604 |
< |
Enhanced (SSD/E - an attempt to improve the liquid structure in |
604 |
> |
Extended (SSD/E - an attempt to improve the liquid structure in |
605 |
|
simulations without a long-range correction). |
606 |
|
|
607 |
|
\begin{table} |
616 |
|
\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ |
617 |
|
\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\ |
618 |
|
\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
619 |
+ |
\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\ |
620 |
|
\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ |
621 |
|
\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\ |
633 |
– |
\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
622 |
|
\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ |
623 |
|
\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ |
624 |
|
\end{tabular} |
642 |
|
\begin{figure} |
643 |
|
\begin{center} |
644 |
|
\epsfxsize=6in |
645 |
< |
\epsfbox{dualsticky.ps} |
645 |
> |
\epsfbox{dualsticky_bw.eps} |
646 |
|
\caption{Isosurfaces of the sticky potential for SSD1 (left) and SSD/E \& |
647 |
|
SSD/RF (right). Light areas correspond to the tetrahedral attractive |
648 |
|
component, and darker areas correspond to the dipolar repulsive |
794 |
|
\begin{center} |
795 |
|
\epsfxsize=6in |
796 |
|
\epsfbox{ssdeDiffuse.epsi} |
797 |
< |
\caption{Plots of the diffusion constants calculated from SSD/E and SSD1, |
798 |
< |
both without a reaction field, along with experimental results |
799 |
< |
[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations were |
800 |
< |
performed at the average densities observed in the 1 atm NPT |
801 |
< |
simulations for the respective models. SSD/E is slightly more fluid |
802 |
< |
than experiment at all of the temperatures, but it is closer than SSD1 |
803 |
< |
without a long-range correction.} |
797 |
> |
\caption{The diffusion constants calculated from SSD/E and SSD1, |
798 |
> |
both without a reaction field, along with experimental results |
799 |
> |
[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations |
800 |
> |
were performed at the average densities observed in the 1 atm NPT |
801 |
> |
simulations for the respective models. SSD/E is slightly more mobile |
802 |
> |
than experiment at all of the temperatures, but it is closer to |
803 |
> |
experiment at biologically relevant temperatures than SSD1 without a |
804 |
> |
long-range correction.} |
805 |
|
\label{ssdediffuse} |
806 |
|
\end{center} |
807 |
|
\end{figure} |
812 |
|
the densities, it is important that the excellent diffusive behavior |
813 |
|
of SSD be maintained or improved. Figure \ref{ssdediffuse} compares |
814 |
|
the temperature dependence of the diffusion constant of SSD/E to SSD1 |
815 |
< |
without an active reaction field, both at the densities calculated at |
816 |
< |
1 atm and at the experimentally calculated densities for super-cooled |
817 |
< |
and liquid water. The diffusion constant for SSD/E is consistently |
818 |
< |
higher than experiment, while SSD1 remains lower than experiment until |
819 |
< |
relatively high temperatures (greater than 330 K). Both models follow |
820 |
< |
the shape of the experimental curve well below 300 K but tend to |
821 |
< |
diffuse too rapidly at higher temperatures, something that is |
833 |
< |
especially apparent with SSD1. This increasing diffusion relative to |
815 |
> |
without an active reaction field at the densities calculated from |
816 |
> |
their respective NPT simulations at 1 atm. The diffusion constant for |
817 |
> |
SSD/E is consistently higher than experiment, while SSD1 remains lower |
818 |
> |
than experiment until relatively high temperatures (around 360 |
819 |
> |
K). Both models follow the shape of the experimental curve well below |
820 |
> |
300 K but tend to diffuse too rapidly at higher temperatures, as seen |
821 |
> |
in SSD1's crossing above 360 K. This increasing diffusion relative to |
822 |
|
the experimental values is caused by the rapidly decreasing system |
823 |
< |
density with increasing temperature. The densities of SSD1 decay more |
824 |
< |
rapidly with temperature than do those of SSD/E, leading to more |
825 |
< |
visible deviation from the experimental diffusion trend. Thus, the |
823 |
> |
density with increasing temperature. Both SSD1 and SSD/E show this |
824 |
> |
deviation in particle mobility, but this trend has different |
825 |
> |
implications on the diffusive behavior of the models. While SSD1 |
826 |
> |
shows more experimentally accurate diffusive behavior in the high |
827 |
> |
temperature regimes, SSD/E shows more accurate behavior in the |
828 |
> |
supercooled and biologically relevant temperature ranges. Thus, the |
829 |
|
changes made to improve the liquid structure may have had an adverse |
830 |
|
affect on the density maximum, but they improve the transport behavior |
831 |
< |
of SSD/E relative to SSD1. |
831 |
> |
of SSD/E relative to SSD1 under the most commonly simulated |
832 |
> |
conditions. |
833 |
|
|
834 |
|
\begin{figure} |
835 |
|
\begin{center} |
836 |
|
\epsfxsize=6in |
837 |
|
\epsfbox{ssdrfDiffuse.epsi} |
838 |
< |
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1, |
838 |
> |
\caption{The diffusion constants calculated from SSD/RF and SSD1, |
839 |
|
both with an active reaction field, along with experimental results |
840 |
< |
[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations |
840 |
> |
[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations |
841 |
|
were performed at the average densities observed in the 1 atm NPT |
842 |
|
simulations for both of the models. Note how accurately SSD/RF |
843 |
|
simulates the diffusion of water throughout this temperature |
844 |
|
range. The more rapidly increasing diffusion constants at high |
845 |
< |
temperatures for both models is attributed to the significantly lower |
846 |
< |
densities than observed in experiment.} |
845 |
> |
temperatures for both models is attributed to lower calculated |
846 |
> |
densities than those observed in experiment.} |
847 |
|
\label{ssdrfdiffuse} |
848 |
|
\end{center} |
849 |
|
\end{figure} |
851 |
|
In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are |
852 |
|
compared to SSD1 with an active reaction field. Note that SSD/RF |
853 |
|
tracks the experimental results quantitatively, identical within error |
854 |
< |
throughout the temperature range shown and with only a slight |
855 |
< |
increasing trend at higher temperatures. SSD1 tends to diffuse more |
856 |
< |
slowly at low temperatures and deviates to diffuse too rapidly at |
854 |
> |
throughout most of the temperature range shown and exhibiting only a |
855 |
> |
slight increasing trend at higher temperatures. SSD1 tends to diffuse |
856 |
> |
more slowly at low temperatures and deviates to diffuse too rapidly at |
857 |
|
temperatures greater than 330 K. As stated above, this deviation away |
858 |
|
from the ideal trend is due to a rapid decrease in density at higher |
859 |
|
temperatures. SSD/RF does not suffer from this problem as much as SSD1 |
861 |
|
values. These results again emphasize the importance of careful |
862 |
|
reparameterization when using an altered long-range correction. |
863 |
|
|
864 |
+ |
\begin{table} |
865 |
+ |
\begin{center} |
866 |
+ |
\caption{Calculated and experimental properties of the single point waters and liquid water at 298 K and 1 atm. (a) Ref. [\citen{Mills73}]. (b) Calculated by integrating the data in ref. \citen{Head-Gordon00_1}. (c) Calculated by integrating the data in ref. \citen{Soper86}. (d) Calculated for 298 K from data in ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.} |
867 |
+ |
\begin{tabular}{ l c c c c c } |
868 |
+ |
\hline \\[-3mm] |
869 |
+ |
\ \ \ \ \ \ & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ |
870 |
+ |
\ & \ SSD/RF \ \ \ & \ Expt. \\ |
871 |
+ |
\hline \\[-3mm] |
872 |
+ |
\ \ \ $\rho$ (g/cm$^3$) & 0.999 $\pm$0.001 & 0.996 $\pm$0.001 & 0.972 $\pm$0.002 & 0.997 $\pm$0.001 & 0.997 \\ |
873 |
+ |
\ \ \ $C_p$ (cal/mol K) & 28.80 $\pm$0.11 & 25.45 $\pm$0.09 & 28.28 $\pm$0.06 & 23.83 $\pm$0.16 & 17.98 \\ |
874 |
+ |
\ \ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78 $\pm$0.07 & 2.51 $\pm$0.18 & 2.00 $\pm$0.17 & 2.32 $\pm$0.06 & 2.299$^\text{a}$ \\ |
875 |
+ |
\ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\ |
876 |
+ |
\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.5$^\text{c}$ \\ |
877 |
+ |
\ \ \ $\tau_1$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 5.7$^\text{d}$ \\ |
878 |
+ |
\ \ \ $\tau_2$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2 $\pm$0.2 & 2.3$^\text{e}$ \\ |
879 |
+ |
\end{tabular} |
880 |
+ |
\label{liquidproperties} |
881 |
+ |
\end{center} |
882 |
+ |
\end{table} |
883 |
+ |
|
884 |
+ |
Table \ref{liquidproperties} gives a synopsis of the liquid state |
885 |
+ |
properties of the water models compared in this study along with the |
886 |
+ |
experimental values for liquid water at ambient conditions. The |
887 |
+ |
coordination number ($N_C$) and hydrogen bonds per particle ($N_H$) |
888 |
+ |
were calculated by integrating the following relations: |
889 |
+ |
\begin{equation} |
890 |
+ |
N_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr, |
891 |
+ |
\end{equation} |
892 |
+ |
\begin{equation} |
893 |
+ |
N_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr, |
894 |
+ |
\end{equation} |
895 |
+ |
where $\rho$ is the number density of the specified pair interactions, |
896 |
+ |
$a$ and $b$ are the radial locations of the minima following the first |
897 |
+ |
solvation shell peak in g$_\text{OO}(r)$ or g$_\text{OH}(r)$ |
898 |
+ |
respectively. The number of hydrogen bonds stays relatively constant |
899 |
+ |
across all of the models, but the coordination numbers of SSD/E and |
900 |
+ |
SSD/RF show an improvement over SSD1. This improvement is primarily |
901 |
+ |
due to the widening of the first solvation shell peak, allowing the |
902 |
+ |
first minima to push outward. Comparing the coordination number with |
903 |
+ |
the number of hydrogen bonds can lead to more insight into the |
904 |
+ |
structural character of the liquid. Because of the near identical |
905 |
+ |
values for SSD1, it appears to be a little too exclusive, in that all |
906 |
+ |
molecules in the first solvation shell are involved in forming ideal |
907 |
+ |
hydrogen bonds. The differing numbers for the newly parameterized |
908 |
+ |
models indicate the allowance of more fluid configurations in addition |
909 |
+ |
to the formation of an acceptable number of ideal hydrogen bonds. |
910 |
+ |
|
911 |
+ |
The time constants for the self orientational autocorrelation function |
912 |
+ |
are also displayed in Table \ref{liquidproperties}. The dipolar |
913 |
+ |
orientational time correlation function ($\Gamma_{l}$) is described |
914 |
+ |
by: |
915 |
+ |
\begin{equation} |
916 |
+ |
\Gamma_{l}(t) = \langle P_l[\mathbf{u}_j(0)\cdot\mathbf{u}_j(t)]\rangle, |
917 |
+ |
\end{equation} |
918 |
+ |
where $P_l$ is a Legendre polynomial of order $l$ and $\mathbf{u}_j$ |
919 |
+ |
is the unit vector of the particle dipole.\cite{Rahman71} From these |
920 |
+ |
correlation functions, the orientational relaxation time of the dipole |
921 |
+ |
vector can be calculated from an exponential fit in the long-time |
922 |
+ |
regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these |
923 |
+ |
time constants were averaged from five detailed NVE simulations |
924 |
+ |
performed at the STP density for each of the respective models. It |
925 |
+ |
should be noted that the commonly cited value for $\tau_2$ of 1.9 ps |
926 |
+ |
was determined from the NMR data in reference \citen{Krynicki66} at a |
927 |
+ |
temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong |
928 |
+ |
temperature dependence of $\tau_2$, it is necessary to recalculate it |
929 |
+ |
at 298 K to make proper comparisons. The value shown in Table |
930 |
+ |
\ref{liquidproperties} was calculated from the same NMR data in the |
931 |
+ |
fashion described in reference \citen{Krynicki66}. Similarly, $\tau_1$ |
932 |
+ |
was recomputed for 298 K from the data in ref \citen{Eisenberg69}. |
933 |
+ |
Again, SSD/E and SSD/RF show improved behavior over SSD1, both with |
934 |
+ |
and without an active reaction field. Turning on the reaction field |
935 |
+ |
leads to much improved time constants for SSD1; however, these results |
936 |
+ |
also include a corresponding decrease in system density. Numbers |
937 |
+ |
published from the original SSD dynamics studies are shorter than the |
938 |
+ |
values observed here, and this difference can be attributed to the use |
939 |
+ |
of the Ewald sum technique versus a reaction field.\cite{Ichiye99} |
940 |
+ |
|
941 |
|
\subsection{Additional Observations} |
942 |
|
|
943 |
|
\begin{figure} |
944 |
|
\begin{center} |
945 |
|
\epsfxsize=6in |
946 |
< |
\epsfbox{povIce.ps} |
946 |
> |
\epsfbox{icei_bw.eps} |
947 |
|
\caption{A water lattice built from the crystal structure assumed by |
948 |
|
SSD/E when undergoing an extremely restricted temperature NPT |
949 |
|
simulation. This form of ice is referred to as ice-{\it i} to |
1023 |
|
significantly lower than the densities obtained from other water |
1024 |
|
models (and experiment). Analysis of self-diffusion showed SSD to |
1025 |
|
capture the transport properties of water well in both the liquid and |
1026 |
< |
super-cooled liquid regimes. |
1026 |
> |
supercooled liquid regimes. |
1027 |
|
|
1028 |
|
In order to correct the density behavior, the original SSD model was |
1029 |
|
reparameterized for use both with and without a reaction field (SSD/RF |
1039 |
|
The existence of a novel low-density ice structure that is preferred |
1040 |
|
by the SSD family of water models is somewhat troubling, since liquid |
1041 |
|
simulations on this family of water models at room temperature are |
1042 |
< |
effectively simulations of super-cooled or metastable liquids. One |
1043 |
< |
way to de-stabilize this unphysical ice structure would be to make the |
1042 |
> |
effectively simulations of supercooled or metastable liquids. One |
1043 |
> |
way to destabilize this unphysical ice structure would be to make the |
1044 |
|
range of angles preferred by the attractive part of the sticky |
1045 |
|
potential much narrower. This would require extensive |
1046 |
|
reparameterization to maintain the same level of agreement with the |