| 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 |
| 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} 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} |
| 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 |
|
|
| 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 |
| 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 at the densities calculated from the |
| 816 |
< |
NPT simulations at 1 atm. The diffusion constant for SSD/E is |
| 817 |
< |
consistently higher than experiment, while SSD1 remains lower than |
| 818 |
< |
experiment until relatively high temperatures (around 360 K). Both |
| 819 |
< |
models follow the shape of the experimental curve well below 300 K but |
| 820 |
< |
tend to diffuse too rapidly at higher temperatures, as seen in SSD1's |
| 821 |
< |
crossing above 360 K. This increasing diffusion relative to the |
| 822 |
< |
experimental values is caused by the rapidly decreasing system density |
| 823 |
< |
with increasing temperature. Both SSD1 and SSD/E show this deviation |
| 824 |
< |
in diffusive behavior, but this trend has different implications on |
| 825 |
< |
the diffusive behavior of the models. While SSD1 shows more |
| 826 |
< |
experimentally accurate diffusive behavior in the high temperature |
| 827 |
< |
regimes, SSD/E shows more accurate behavior in the supercooled and |
| 828 |
< |
biologically relevant temperature ranges. Thus, the changes made to |
| 829 |
< |
improve the liquid structure may have had an adverse affect on the |
| 830 |
< |
density maximum, but they improve the transport behavior of SSD/E |
| 831 |
< |
relative to SSD1 under the most commonly simulated conditions. |
| 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. 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 under the most commonly simulated |
| 832 |
> |
conditions. |
| 833 |
|
|
| 834 |
|
\begin{figure} |
| 835 |
|
\begin{center} |
| 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) Ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.} |
| 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) \ \ |
| 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.4$^\text{c}$ \\ |
| 877 |
< |
\ \ \ $\tau_1^\mu$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 4.76$^\text{d}$ \\ |
| 878 |
< |
\ \ \ $\tau_2^\mu$ (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}$ \\ |
| 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} |
| 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 and hydrogen bonds per particle were calculated by |
| 888 |
< |
integrating the following relation: |
| 887 |
> |
coordination number ($N_C$) and hydrogen bonds per particle ($N_H$) |
| 888 |
> |
were calculated by integrating the following relations: |
| 889 |
|
\begin{equation} |
| 890 |
< |
4\pi\rho\int_{0}^{a}r^2\text{g}(r)dr, |
| 890 |
> |
N_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr, |
| 891 |
|
\end{equation} |
| 892 |
< |
where $\rho$ is the number density of pair interactions, $a$ is the |
| 893 |
< |
radial location of the minima following the first solvation shell |
| 894 |
< |
peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for |
| 895 |
< |
calculation of the coordination number or hydrogen bonds per particle |
| 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 |
| 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^\mu$).\cite{Rothschild84} Calculation of these |
| 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{Rahman73} Because of the strong |
| 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}. Again, SSD/E and |
| 932 |
< |
SSD/RF show improved behavior over SSD1, both with and without an |
| 933 |
< |
active reaction field. Turning on the reaction field leads to much |
| 934 |
< |
improved time constants for SSD1; however, these results also include |
| 935 |
< |
a corresponding decrease in system density. Numbers published from the |
| 936 |
< |
original SSD dynamics studies appear closer to the experimental |
| 937 |
< |
values, and this difference can be attributed to the use of the Ewald |
| 938 |
< |
sum technique versus a reaction field.\cite{Ichiye99} |
| 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 |
