| 221 |
|
\caption{Energy conservation using both quaternion-based integration |
| 222 |
|
and the {\sc dlm} method with increasing time step. The larger time |
| 223 |
|
step plots are shifted from the true energy baseline (that of $\Delta |
| 224 |
< |
t$ = 0.1fs) for clarity.} |
| 224 |
> |
t$ = 0.1~fs) for clarity.} |
| 225 |
|
\label{fig:timeStepIntegration} |
| 226 |
|
\end{figure} |
| 227 |
|
|
| 239 |
|
at various time steps for both {\sc dlm} and quaternion |
| 240 |
|
integration. All of the 1000 SSD particle simulations started with the |
| 241 |
|
same configuration, and the only difference was the method used to |
| 242 |
< |
handle orientational motion. At time steps of 0.1 and 0.5fs, both |
| 242 |
> |
handle orientational motion. At time steps of 0.1 and 0.5~fs, both |
| 243 |
|
methods for propagating the orientational degrees of freedom conserve |
| 244 |
|
energy fairly well, with the quaternion method showing a slight energy |
| 245 |
< |
drift over time in the 0.5fs time step simulation. Time steps of 1 and |
| 246 |
< |
2fs clearly demonstrate the benefits in energy conservation that come |
| 245 |
> |
drift over time in the 0.5~fs time step simulation. Time steps of 1 and |
| 246 |
> |
2~fs clearly demonstrate the benefits in energy conservation that come |
| 247 |
|
with the {\sc dlm} method. Thus, while maintaining the same degree of |
| 248 |
|
energy conservation, one can take considerably longer time steps, |
| 249 |
|
leading to an overall reduction in computation time. |
| 250 |
|
|
| 251 |
|
Energy drifts in water simulations using {\sc dlm} integration were |
| 252 |
< |
unnoticeable for time steps up to 3fs. We observed a slight energy |
| 253 |
< |
drift on the order of 0.012~kcal/mol per nanosecond with a time step |
| 254 |
< |
of 4fs. As expected, this drift increases dramatically with increasing |
| 252 |
> |
unnoticeable for time steps up to 3~fs. We observed a slight energy |
| 253 |
> |
drift on the order of 0.012~kcal/mol per nanosecond with a time step |
| 254 |
> |
of 4~fs. As expected, this drift increases dramatically with increasing |
| 255 |
|
time step. To insure accuracy in our {\it NVE} simulations, time steps |
| 256 |
< |
were set at 2fs and were also kept at this value for {\it NPT} |
| 256 |
> |
were set at 2~fs and were also kept at this value for {\it NPT} |
| 257 |
|
simulations. |
| 258 |
|
|
| 259 |
|
Proton-disordered ice crystals in both the I$_\textrm{h}$ and |
| 265 |
|
orthorhombic in shape with an edge length ratio of approximately |
| 266 |
|
1.00$\times$1.06$\times$1.23. We then allowed the particles to orient |
| 267 |
|
freely about their fixed lattice positions with angular momenta values |
| 268 |
< |
randomly sampled at 400K. The rotational temperature was then scaled |
| 269 |
< |
down in stages to slowly cool the crystals to 25K. The particles were |
| 268 |
> |
randomly sampled at 400~K. The rotational temperature was then scaled |
| 269 |
> |
down in stages to slowly cool the crystals to 25~K. The particles were |
| 270 |
|
then allowed to translate with fixed orientations at a constant |
| 271 |
< |
pressure of 1atm for 50ps at 25K. Finally, all constraints were |
| 272 |
< |
removed and the ice crystals were allowed to equilibrate for 50ps at |
| 273 |
< |
25K and a constant pressure of 1atm. This procedure resulted in |
| 271 |
> |
pressure of 1atm for 50~ps at 25~K. Finally, all constraints were |
| 272 |
> |
removed and the ice crystals were allowed to equilibrate for 50~ps at |
| 273 |
> |
25~K and a constant pressure of 1atm. This procedure resulted in |
| 274 |
|
structurally stable I$_\textrm{h}$ ice crystals that obey the |
| 275 |
|
Bernal-Fowler rules.\cite{Bernal33,Rahman72} This method was also |
| 276 |
|
utilized in the making of diamond lattice I$_\textrm{c}$ ice crystals, |
| 294 |
|
|
| 295 |
|
An ensemble average from five separate melting simulations was |
| 296 |
|
acquired, each starting from different ice crystals generated as |
| 297 |
< |
described previously. All simulations were equilibrated for 100ps |
| 298 |
< |
prior to a 200ps data collection run at each temperature setting. The |
| 299 |
< |
temperature range of study spanned from 25 to 400K, with a maximum |
| 300 |
< |
degree increment of 25K. For regions of interest along this stepwise |
| 301 |
< |
progression, the temperature increment was decreased from 25K to 10 |
| 302 |
< |
and 5K. The above equilibration and production times were sufficient |
| 297 |
> |
described previously. All simulations were equilibrated for 100~ps |
| 298 |
> |
prior to a 200~ps data collection run at each temperature setting. The |
| 299 |
> |
temperature range of study spanned from 25 to 400~K, with a maximum |
| 300 |
> |
degree increment of 25~K. For regions of interest along this stepwise |
| 301 |
> |
progression, the temperature increment was decreased from 25~K to 10 |
| 302 |
> |
and 5~K. The above equilibration and production times were sufficient |
| 303 |
|
in that the fluctuations in the volume autocorrelation function damped |
| 304 |
< |
out in all of the simulations in under 20ps. |
| 304 |
> |
out in all of the simulations in under 20~ps. |
| 305 |
|
|
| 306 |
|
Our initial simulations focused on the original SSD water model, and |
| 307 |
|
an average density versus temperature plot is shown in figure |
| 308 |
|
\ref{fig:ssdDense}. Note that the density maximum when using a |
| 309 |
< |
reaction field appears between 255 and 265K. There were smaller |
| 310 |
< |
fluctuations in the density at 260K than at either 255 or 265K, so we |
| 309 |
> |
reaction field appears between 255 and 265~K. There were smaller |
| 310 |
> |
fluctuations in the density at 260~K than at either 255 or 265~K, so we |
| 311 |
|
report this value as the location of the density maximum. Figure |
| 312 |
|
\ref{fig:ssdDense} was constructed using ice I$_\textrm{h}$ crystals |
| 313 |
|
for the initial configuration; though not pictured, the simulations |
| 342 |
|
the cutoff radius ($R_\textrm{c}$), both with and without the use of a |
| 343 |
|
reaction field.\cite{vanderSpoel98} In order to address the possible |
| 344 |
|
effect of $R_\textrm{c}$, simulations were performed with a cutoff |
| 345 |
< |
radius of 12\AA\, complementing the 9\AA\ $R_\textrm{c}$ used in the |
| 345 |
> |
radius of 12~\AA\, complementing the 9~\AA\ $R_\textrm{c}$ used in the |
| 346 |
|
previous SSD simulations. All of the resulting densities overlapped |
| 347 |
|
within error and showed no significant trend toward lower or higher |
| 348 |
|
densities in simulations both with and without reaction field. |
| 351 |
|
density scaling of SSD relative to other common models at any given |
| 352 |
|
temperature. SSD assumes a lower density than any of the other listed |
| 353 |
|
models at the same pressure, behavior which is especially apparent at |
| 354 |
< |
temperatures greater than 300K. Lower than expected densities have |
| 354 |
> |
temperatures greater than 300~K. Lower than expected densities have |
| 355 |
|
been observed for other systems using a reaction field for long-range |
| 356 |
|
electrostatic interactions, so the most likely reason for the reduced |
| 357 |
|
densities is the presence of the reaction |
| 364 |
|
water. The shape of the curve is similar to the curve produced from |
| 365 |
|
SSD simulations using reaction field, specifically the rapidly |
| 366 |
|
decreasing densities at higher temperatures; however, a shift in the |
| 367 |
< |
density maximum location, down to 245K, is observed. This is a more |
| 367 |
> |
density maximum location, down to 245~K, is observed. This is a more |
| 368 |
|
accurate comparison to the other listed water models, in that no long |
| 369 |
|
range corrections were applied in those |
| 370 |
|
simulations.\cite{Baez94,Jorgensen98b} However, even without the |
| 371 |
< |
reaction field, the density around 300K is still significantly lower |
| 371 |
> |
reaction field, the density around 300~K is still significantly lower |
| 372 |
|
than experiment and comparable water models. This anomalous behavior |
| 373 |
|
was what lead Tan {\it et al.} to recently reparametrize |
| 374 |
|
SSD.\cite{Tan03} Throughout the remainder of the paper our |
| 383 |
|
NVE} simulations were performed at the average densities obtained from |
| 384 |
|
the {\it NPT} simulations at an identical target |
| 385 |
|
temperature. Simulations started with randomized velocities and |
| 386 |
< |
underwent 50ps of temperature scaling and 50ps of constant energy |
| 387 |
< |
equilibration before a 200ps data collection run. Diffusion constants |
| 386 |
> |
underwent 50~ps of temperature scaling and 50~ps of constant energy |
| 387 |
> |
equilibration before a 200~ps data collection run. Diffusion constants |
| 388 |
|
were calculated via linear fits to the long-time behavior of the |
| 389 |
|
mean-square displacement as a function of time.\cite{Allen87} The |
| 390 |
|
averaged results from five sets of {\it NVE} simulations are displayed |
| 408 |
|
strengths of the SSD model. Of the three models shown, the SSD model |
| 409 |
|
has the most accurate depiction of self-diffusion in both the |
| 410 |
|
supercooled and liquid regimes. SPC/E does a respectable job by |
| 411 |
< |
reproducing values similar to experiment around 290K; however, it |
| 411 |
> |
reproducing values similar to experiment around 290~K; however, it |
| 412 |
|
deviates at both higher and lower temperatures, failing to predict the |
| 413 |
|
correct thermal trend. TIP5P and SSD both start off low at colder |
| 414 |
|
temperatures and tend to diffuse too rapidly at higher temperatures. |
| 427 |
|
pressure heat capacity ($C_\textrm{p}$) was monitored to locate |
| 428 |
|
$T_\textrm{m}$ in each of the simulations. In the melting simulations |
| 429 |
|
of the 1024 particle ice I$_\textrm{h}$ simulations, a large spike in |
| 430 |
< |
$C_\textrm{p}$ occurs at 245K, indicating a first order phase |
| 430 |
> |
$C_\textrm{p}$ occurs at 245~K, indicating a first order phase |
| 431 |
|
transition for the melting of these ice crystals (see figure |
| 432 |
|
\ref{fig:ssdCp}. When the reaction field is turned off, the melting |
| 433 |
< |
transition occurs at 235K. These melting transitions are considerably |
| 434 |
< |
lower than the experimental value of 273K, indicating that the solid |
| 433 |
> |
transition occurs at 235~K. These melting transitions are considerably |
| 434 |
> |
lower than the experimental value of 273~K, indicating that the solid |
| 435 |
|
ice I$_\textrm{h}$ is not thermodynamically preferred relative to the |
| 436 |
|
liquid state at these lower temperatures. |
| 437 |
|
\begin{figure} |
| 438 |
|
\centering |
| 439 |
|
\includegraphics[width=\linewidth]{./figures/ssdCp.pdf} |
| 440 |
|
\caption{Heat capacity versus temperature for the SSD model with an |
| 441 |
< |
active reaction field. Note the large spike in $C_p$ around 245K, |
| 441 |
> |
active reaction field. Note the large spike in $C_p$ around 245~K, |
| 442 |
|
indicating a phase transition from the ordered crystal to disordered |
| 443 |
|
liquid.} |
| 444 |
|
\label{fig:ssdCp} |
| 448 |
|
\centering |
| 449 |
|
\includegraphics[width=\linewidth]{./figures/fullContour.pdf} |
| 450 |
|
\caption{ Contour plots of 2D angular pair correlation functions for |
| 451 |
< |
512 SSD molecules at 100K (A \& B) and 300K (C \& D). Dark areas |
| 451 |
> |
512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas |
| 452 |
|
signify regions of enhanced density while light areas signify |
| 453 |
|
depletion relative to the bulk density. White areas have pair |
| 454 |
|
correlation values below 0.5 and black areas have values above 1.5.} |
| 502 |
|
$g_\textrm{OO}(r)$.\cite{Liu96b} At low temperatures, the second |
| 503 |
|
solvation shell peak appears to have two distinct components that |
| 504 |
|
blend together to form one observable peak. At higher temperatures, |
| 505 |
< |
this split character alters to show the leading 4\AA\ peak dominated |
| 505 |
> |
this split character alters to show the leading 4~\AA\ peak dominated |
| 506 |
|
by equatorial anti-parallel dipole orientations. There is also a |
| 507 |
|
tightly bunched group of axially arranged dipoles that most likely |
| 508 |
|
consist of the smaller fraction of aligned dipole pairs. The trailing |
| 509 |
< |
component of the split peak at 5\AA\ is dominated by aligned dipoles |
| 509 |
> |
component of the split peak at 5~\AA\ is dominated by aligned dipoles |
| 510 |
|
that assume hydrogen bond arrangements similar to those seen in the |
| 511 |
|
first solvation shell. This evidence indicates that the dipole pair |
| 512 |
|
interaction begins to dominate outside of the range of the dipolar |
| 513 |
|
repulsion term. The energetically favorable dipole arrangements |
| 514 |
|
populate the region immediately outside this repulsion region (around |
| 515 |
< |
4\AA), while arrangements that seek to satisfy both the sticky and |
| 515 |
> |
4~\AA), while arrangements that seek to satisfy both the sticky and |
| 516 |
|
dipole forces locate themselves just beyond this initial buildup |
| 517 |
< |
(around 5\AA). |
| 517 |
> |
(around 5~\AA). |
| 518 |
|
|
| 519 |
|
This analysis indicates that the split second peak is primarily the |
| 520 |
|
product of the dipolar repulsion term of the sticky potential. In |
| 521 |
|
fact, the inner peak can be pushed out and merged with the outer split |
| 522 |
|
peak just by extending the switching function ($s^\prime(r_{ij})$) |
| 523 |
< |
from its normal 4\AA\ cutoff to values of 4.5 or even 5\AA. This |
| 523 |
> |
from its normal 4~\AA\ cutoff to values of 4.5 or even 5~\AA. This |
| 524 |
|
type of correction is not recommended for improving the liquid |
| 525 |
|
structure, since the second solvation shell would still be shifted too |
| 526 |
|
far out. In addition, this would have an even more detrimental effect |
| 557 |
|
\caption{PARAMETERS FOR THE ORIGINAL AND ADJUSTED SSD MODELS} |
| 558 |
|
|
| 559 |
|
\centering |
| 560 |
< |
\begin{tabular}{ lcccc } |
| 560 |
> |
\begin{tabular}{ llcccc } |
| 561 |
|
\toprule |
| 562 |
|
\toprule |
| 563 |
< |
Parameters & SSD & SSD1 & SSD/E & SSD/RF \\ |
| 563 |
> |
Parameters & & SSD & SSD1 & SSD/E & SSD/RF \\ |
| 564 |
|
\midrule |
| 565 |
< |
$\sigma$ (\AA) & 3.051 & 3.016 & 3.035 & 3.019\\ |
| 566 |
< |
$\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ |
| 567 |
< |
$\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\ |
| 568 |
< |
$\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
| 569 |
< |
$\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\ |
| 570 |
< |
$r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ |
| 571 |
< |
$r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\ |
| 572 |
< |
$r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ |
| 573 |
< |
$r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ |
| 565 |
> |
$\sigma$ & (\AA) & 3.051 & 3.016 & 3.035 & 3.019\\ |
| 566 |
> |
$\epsilon$ & (kcal mol$^{-1}$) & 0.152 & 0.152 & 0.152 & 0.152\\ |
| 567 |
> |
$\mu$ & (D) & 2.35 & 2.35 & 2.42 & 2.48\\ |
| 568 |
> |
$\nu_0$ & (kcal mol$^{-1}$) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
| 569 |
> |
$\omega^\circ$ & & 0.07715 & 0.07715 & 0.07715 & 0.07715\\ |
| 570 |
> |
$r_l$ & (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ |
| 571 |
> |
$r_u$ & (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\ |
| 572 |
> |
$r_l^\prime$ & (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ |
| 573 |
> |
$r_u^\prime$ & (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ |
| 574 |
|
\bottomrule |
| 575 |
|
\end{tabular} |
| 576 |
|
\label{tab:ssdParams} |
| 632 |
|
bonds''. Aside from improving the shape of the first peak in the |
| 633 |
|
$g(r)$, this modification improves the densities considerably by |
| 634 |
|
allowing the persistence of full dipolar character below the previous |
| 635 |
< |
4\AA\ cutoff. |
| 635 |
> |
4~\AA\ cutoff. |
| 636 |
|
|
| 637 |
|
While adjusting the location and shape of the first peak of $g(r)$ |
| 638 |
|
improves the densities, these changes alone are insufficient to bring |
| 664 |
|
simulations of 1024 particle systems, and the melting sequences were |
| 665 |
|
started from different ice I$_\textrm{h}$ crystals constructed as |
| 666 |
|
described previously. Each {\it NPT} simulation was equilibrated for |
| 667 |
< |
100ps before a 200ps data collection run at each temperature step, |
| 667 |
> |
100~ps before a 200~ps data collection run at each temperature step, |
| 668 |
|
and the final configuration from the previous temperature simulation |
| 669 |
|
was used as a starting point. All {\it NVE} simulations had the same |
| 670 |
|
thermalization, equilibration, and data collection times as stated |
| 686 |
|
calculated densities for both SSD/E and SSD1 have increased |
| 687 |
|
significantly over the original SSD model (see figure |
| 688 |
|
\ref{fig:ssdDense}) and are in better agreement with the experimental |
| 689 |
< |
values. At 298 K, the densities of SSD/E and SSD1 without a long-range |
| 690 |
< |
correction are 0.996 g/cm$^3$ and 0.999 g/cm$^3$ respectively. These |
| 691 |
< |
both compare well with the experimental value of 0.997 g/cm$^3$, and |
| 692 |
< |
they are considerably better than the SSD value of 0.967 g/cm$^3$. The |
| 689 |
> |
values. At 298~K, the densities of SSD/E and SSD1 without a long-range |
| 690 |
> |
correction are 0.996~g/cm$^3$ and 0.999~g/cm$^3$ respectively. These |
| 691 |
> |
both compare well with the experimental value of 0.997~g/cm$^3$, and |
| 692 |
> |
they are considerably better than the SSD value of 0.967~g/cm$^3$. The |
| 693 |
|
changes to the dipole moment and sticky switching functions have |
| 694 |
|
improved the structuring of the liquid (as seen in figure |
| 695 |
|
\ref{fig:gofrCompare}), but they have shifted the density maximum to |
| 698 |
|
strengthening of the dipolar character. However, this increasing |
| 699 |
|
disorder in the SSD/E model has little effect on the melting |
| 700 |
|
transition. By monitoring $C_p$ throughout these simulations, we |
| 701 |
< |
observed a melting transition for SSD/E at 235K, the same as SSD and |
| 701 |
> |
observed a melting transition for SSD/E at 235~K, the same as SSD and |
| 702 |
|
SSD1. |
| 703 |
|
|
| 704 |
|
\begin{figure} |
| 719 |
|
with an active reaction field is shown in figure \ref{fig:ssdrfDense}. |
| 720 |
|
As observed in the simulations without a reaction field, the densities |
| 721 |
|
of SSD/RF and SSD1 show a dramatic increase over normal SSD (see |
| 722 |
< |
figure \ref{fig:ssdDense}). At 298 K, SSD/RF has a density of 0.997 |
| 722 |
> |
figure \ref{fig:ssdDense}). At 298~K, SSD/RF has a density of 0.997 |
| 723 |
|
g/cm$^3$, directly in line with experiment and considerably better |
| 724 |
< |
than the original SSD value of 0.941 g/cm$^3$ and the SSD1 value of |
| 725 |
< |
0.972 g/cm$^3$. These results further emphasize the importance of |
| 724 |
> |
than the original SSD value of 0.941~g/cm$^3$ and the SSD1 value of |
| 725 |
> |
0.972~g/cm$^3$. These results further emphasize the importance of |
| 726 |
|
reparametrization in order to model the density properly under |
| 727 |
|
different simulation conditions. Again, these changes have only a |
| 728 |
< |
minor effect on the melting point, which observed at 245K for SSD/RF, |
| 729 |
< |
is identical to SSD and only 5K lower than SSD1 with a reaction |
| 728 |
> |
minor effect on the melting point, which observed at 245~K for SSD/RF, |
| 729 |
> |
is identical to SSD and only 5~K lower than SSD1 with a reaction |
| 730 |
|
field. Additionally, the difference in density maxima is not as |
| 731 |
< |
extreme, with SSD/RF showing a density maximum at 255K, fairly close |
| 732 |
< |
to the density maxima of 260K and 265K, shown by SSD and SSD1 |
| 731 |
> |
extreme, with SSD/RF showing a density maximum at 255~K, fairly close |
| 732 |
> |
to the density maxima of 260~K and 265~K, shown by SSD and SSD1 |
| 733 |
|
respectively. |
| 734 |
|
|
| 735 |
|
\subsection{Transport Behavior} |
| 758 |
|
from their respective {\it NPT} simulations at 1 atm. The diffusion |
| 759 |
|
constant for SSD/E is consistently higher than experiment, while SSD1 |
| 760 |
|
remains lower than experiment until relatively high temperatures |
| 761 |
< |
(around 360K). Both models follow the shape of the experimental curve |
| 762 |
< |
below 300K but tend to diffuse too rapidly at higher temperatures, as |
| 763 |
< |
seen in SSD1 crossing above 360K. This increasing diffusion relative |
| 761 |
> |
(around 360~K). Both models follow the shape of the experimental curve |
| 762 |
> |
below 300~K but tend to diffuse too rapidly at higher temperatures, as |
| 763 |
> |
seen in SSD1 crossing above 360~K. This increasing diffusion relative |
| 764 |
|
to the experimental values is caused by the rapidly decreasing system |
| 765 |
|
density with increasing temperature. Both SSD1 and SSD/E show this |
| 766 |
|
deviation in particle mobility, but this trend has different |
| 793 |
|
throughout most of the temperature range shown and exhibiting only a |
| 794 |
|
slight increasing trend at higher temperatures. SSD1 tends to diffuse |
| 795 |
|
more slowly at low temperatures and deviates to diffuse too rapidly at |
| 796 |
< |
temperatures greater than 330K. As stated above, this deviation away |
| 796 |
> |
temperatures greater than 330~K. As stated above, this deviation away |
| 797 |
|
from the ideal trend is due to a rapid decrease in density at higher |
| 798 |
|
temperatures. SSD/RF does not suffer from this problem as much as SSD1 |
| 799 |
|
because the calculated densities are closer to the experimental |
| 866 |
|
regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these time |
| 867 |
|
constants were averaged over five detailed {\it NVE} simulations |
| 868 |
|
performed at the ambient conditions for each of the respective |
| 869 |
< |
models. It should be noted that the commonly cited value of 1.9 ps for |
| 869 |
> |
models. It should be noted that the commonly cited value of 1.9~ps for |
| 870 |
|
$\tau_2$ was determined from the NMR data in Ref. \cite{Krynicki66} at |
| 871 |
|
a temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong |
| 872 |
|
temperature dependence of $\tau_2$, it is necessary to recalculate it |
| 873 |
< |
at 298K to make proper comparisons. The value shown in Table |
| 873 |
> |
at 298~K to make proper comparisons. The value shown in Table |
| 874 |
|
\ref{tab:liquidProperties} was calculated from the same NMR data in the |
| 875 |
|
fashion described in Ref. \cite{Krynicki66}. Similarly, $\tau_1$ was |
| 876 |
< |
recomputed for 298K from the data in Ref. \cite{Eisenberg69}. |
| 876 |
> |
recomputed for 298~K from the data in Ref. \cite{Eisenberg69}. |
| 877 |
|
Again, SSD/E and SSD/RF show improved behavior over SSD1, both with |
| 878 |
|
and without an active reaction field. Turning on the reaction field |
| 879 |
|
leads to much improved time constants for SSD1; however, these results |
| 882 |
|
paper are smaller than the values observed here, and this difference |
| 883 |
|
can be attributed to the use of the Ewald sum.\cite{Chandra99} |
| 884 |
|
|
| 885 |
< |
\subsection{SSD/RF and Damped Electrostatics} |
| 885 |
> |
\subsection{SSD/RF and Damped Electrostatics}\label{sec:ssdrfDamped} |
| 886 |
|
|
| 887 |
|
In section \ref{sec:dampingMultipoles}, a method was described for |
| 888 |
|
applying the damped {\sc sf} or {\sc sp} techniques to for systems |
| 906 |
|
\toprule |
| 907 |
|
\toprule |
| 908 |
|
& & Reaction Field & Damped Electrostatics & Experiment [Ref.] \\ |
| 909 |
< |
& & $\epsilon = 80$ & $\alpha = 0.2125$\AA$^{-1}$ & \\ |
| 909 |
> |
& & $\epsilon = 80$ & $\alpha = 0.2125$~\AA$^{-1}$ & \\ |
| 910 |
|
\midrule |
| 911 |
|
$\rho$ & (g cm$^{-3}$) & 0.997(1) & 1.004(4) & 0.997 \cite{CRC80}\\ |
| 912 |
|
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 23.8(2) & 27(1) & 18.005 \cite{Wagner02} \\ |
| 913 |
|
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.32(6) & 2.33(2) & 2.299 \cite{Mills73}\\ |
| 914 |
< |
$n_C$ & & 4.4 & 4.4 & 4.7 \cite{Hura00}\\ |
| 914 |
> |
$n_C$ & & 4.4 & 4.2 & 4.7 \cite{Hura00}\\ |
| 915 |
|
$n_H$ & & 3.7 & 3.7 & 3.5 \cite{Soper86}\\ |
| 916 |
|
$\tau_1$ & (ps) & 7.2(4) & 5.86(8) & 5.7 \cite{Eisenberg69}\\ |
| 917 |
|
$\tau_2$ & (ps) & 3.2(2) & 2.45(7) & 2.3 \cite{Krynicki66}\\ |
| 937 |
|
encouraging result because of the more varied applicability of damping |
| 938 |
|
over the reaction field technique. Rather than be limited to |
| 939 |
|
homogeneous systems, SSD/RF can be used effectively with mixed |
| 940 |
< |
systems, such as dissolved ions, small organic molecules, or even |
| 940 |
> |
systems, such as dissolved ions, dissolved organic molecules, or even |
| 941 |
|
proteins. |
| 942 |
|
|
| 943 |
|
In addition to the properties tabulated in table |
| 944 |
|
\ref{tab:dampedSSDRF}, we calculated the static dielectric constant |
| 945 |
< |
from a 5ns simulation of SSD/RF using the damped electrostatics. The |
| 945 |
> |
from a 5~ns simulation of SSD/RF using the damped electrostatics. The |
| 946 |
|
resulting value of 82.6(6) compares very favorably with the |
| 947 |
|
experimental value of 78.3.\cite{Malmberg56} This value is closer to |
| 948 |
|
the experimental value than what was expected according to figure |
| 952 |
|
|
| 953 |
|
\section{Tetrahedrally Restructured Elongated Dipole (TRED) Water Model} |
| 954 |
|
|
| 955 |
< |
The SSD/RF model works well with damped electrostatics, but because of its point multipole character, there is no charge neutralization correction at $R_\textrm{c}$. This has the effect of increasing the density, since there is no consideration of the ``surroundings''. |
| 955 |
> |
The SSD/RF model works well with damped electrostatics, but because of |
| 956 |
> |
its point multipole character, there is no charge neutralization |
| 957 |
> |
correction at $R_\textrm{c}$. This has the effect of increasing the |
| 958 |
> |
density, since there is no consideration of the ``surroundings''. In |
| 959 |
> |
an attempt to optimize a water model for these conditions, we decided |
| 960 |
> |
to both simplify the parameters of the SSD type models and break the |
| 961 |
> |
point dipole into a charge pair so that it will gain some effect from |
| 962 |
> |
the shifting action in the {\sc sf} technique. This process resulted |
| 963 |
> |
in a two-point model that we are calling the Tetrahedrally |
| 964 |
> |
Restructured Elongated Dipole (TRED) water model. |
| 965 |
|
|
| 966 |
+ |
\begin{figure} |
| 967 |
+ |
\centering |
| 968 |
+ |
\includegraphics[width=2.75in]{./figures/tredLayout.pdf} |
| 969 |
+ |
\caption{Geometry of TRED water. In this two-point model, the origin |
| 970 |
+ |
is the center-of-mass of the water molecule with the same moments of |
| 971 |
+ |
inertia as SSD/RF. The negatively charged site at the origin is also |
| 972 |
+ |
both a Lennard-Jones and sticky interaction site.} |
| 973 |
+ |
\label{fig:tredLayout} |
| 974 |
+ |
\end{figure} |
| 975 |
|
\begin{table} |
| 976 |
+ |
\centering |
| 977 |
+ |
\caption{PARAMETERS FOR THE TRED WATER MODEL} |
| 978 |
+ |
\begin{tabular}{ llr } |
| 979 |
+ |
\toprule |
| 980 |
+ |
\toprule |
| 981 |
+ |
$\sigma$ & (\AA) & 2.980 \\ |
| 982 |
+ |
$\epsilon$ & (kcal mol$^{-1}$) & 0.2045 \\ |
| 983 |
+ |
$q$ & ($e$) & 1.041 \\ |
| 984 |
+ |
$v_0$ & (kcal mol$^{-1}$) & 4.22 \\ |
| 985 |
+ |
$\omega^\circ$ & & 0.07715 \\ |
| 986 |
+ |
$r_l$ \& $r_l^\prime$ & (\AA) & 2.4 \\ |
| 987 |
+ |
$r_u$ \& $r_u^\prime$ & (\AA) & 4.0 \\ |
| 988 |
+ |
Q$_xx$ & (D \AA) & -1.682 \\ |
| 989 |
+ |
Q$_yy$ & (D \AA) & 1.762 \\ |
| 990 |
+ |
Q$_zz$ & (D \AA) & -0.080 \\ |
| 991 |
+ |
\bottomrule |
| 992 |
+ |
\end{tabular} |
| 993 |
+ |
\label{tab:tredParams} |
| 994 |
+ |
\end{table} |
| 995 |
+ |
\begin{figure} |
| 996 |
+ |
\centering |
| 997 |
+ |
\includegraphics[width=\linewidth]{./figures/tredGofR.pdf} |
| 998 |
+ |
\caption{The $g_\textrm{OO}(r)$ for TRED water. The first peak has a |
| 999 |
+ |
closer to the experimental plot; however, the second and third peaks |
| 1000 |
+ |
exhibit a more expanded structure similar to SSD/RF. The minimum |
| 1001 |
+ |
following the first peak is located at 3.55~\AA , which is further out |
| 1002 |
+ |
than both experiment and SSD/RF (3.42 and 3.3~\AA\ respectively), |
| 1003 |
+ |
leading to a larger coordination number. If all the curves were |
| 1004 |
+ |
integrated to the experimental minimum, the $n_C$ results would all be |
| 1005 |
+ |
similar, with TRED having an $n_C$ only 0.2 larger than experiment.} |
| 1006 |
+ |
\label{fig:tredGofR} |
| 1007 |
+ |
\end{figure} |
| 1008 |
+ |
The geometric layout of TRED water is shown in figure |
| 1009 |
+ |
\ref{fig:tredLayout} and the electrostatic, Lennard-Jones, and sticky |
| 1010 |
+ |
parameters are displayed in table \ref{tab:tredParams}. The negatively |
| 1011 |
+ |
charged site is located at the center of mass of the molecule and is |
| 1012 |
+ |
also home to the Lennard-Jones and sticky interactions. The charge |
| 1013 |
+ |
separation distance of 0.5~\AA\ along the dipolar ($z$) axis combined |
| 1014 |
+ |
with the charge magnitude ($q$) results in a dipole moment of |
| 1015 |
+ |
2.5~D. This value is simply the value used for SSD/RF (2.48~D) rounded |
| 1016 |
+ |
to the tenths place. We also unified the sticky parameters for the |
| 1017 |
+ |
switching functions on the repulsive and attractive interactions in |
| 1018 |
+ |
the interest of simplicity, and we left the quadrupole moment elements |
| 1019 |
+ |
and $\omega^\circ$ unaltered. Finally, the strength of the sticky |
| 1020 |
+ |
interaction ($v_0$) was modified to improve the shape of the first |
| 1021 |
+ |
peaks in $g_\textrm{OO}(r)$ and $g_\textrm{OH}(r)$, while the $\sigma$ |
| 1022 |
+ |
and $\epsilon$ values were varied to adjust the location of the first |
| 1023 |
+ |
peak of $g_\textrm{OO}(r)$ (see figure \ref{fig:tredGofR}) and the |
| 1024 |
+ |
density. The $\sigma$ and $epsilon$ optimization was carried out by |
| 1025 |
+ |
separating the Lennard-Jones potential into near entirely repulsive |
| 1026 |
+ |
($A$) and attractive ($C$) parts, such that |
| 1027 |
+ |
\begin{equation} |
| 1028 |
+ |
\sigma = \left(\frac{A}{C}\right)^\frac{1}{6}, |
| 1029 |
+ |
\end{equation} |
| 1030 |
+ |
and |
| 1031 |
+ |
\begin{equation} |
| 1032 |
+ |
\epsilon = \frac{C^2}{4A}. |
| 1033 |
+ |
\end{equation} |
| 1034 |
+ |
The location of the peak is adjusted through $A$, while the |
| 1035 |
+ |
interaction strength is modified through $C$. All of the above changes |
| 1036 |
+ |
were made with the goal of capturing the experimental density and |
| 1037 |
+ |
translational diffusion constant at 298~K and 1~atm. |
| 1038 |
+ |
|
| 1039 |
+ |
Being that TRED is a two-point water model, we expect its |
| 1040 |
+ |
computational efficiency to fall some place in between the single and |
| 1041 |
+ |
three-point water models. In comparative simulations, TRED was |
| 1042 |
+ |
approximately 85\% slower than SSD/RF, while SPC/E was 225\% slower |
| 1043 |
+ |
than SSD/RF. While TRED loses some of the performance advantage of |
| 1044 |
+ |
SSD, it is still nearly twice as computationally efficient as SPC/E |
| 1045 |
+ |
and TIP3P. |
| 1046 |
+ |
|
| 1047 |
+ |
\begin{table} |
| 1048 |
|
\caption{PROPERTIES OF TRED COMPARED WITH SSD/RF AND EXPERIMENT} |
| 1049 |
|
\footnotesize |
| 1050 |
|
\centering |
| 1052 |
|
\toprule |
| 1053 |
|
\toprule |
| 1054 |
|
& & SSD/RF & TRED & Experiment [Ref.]\\ |
| 1055 |
< |
& & $\alpha = 0.2125$\AA$^{-1}$ & $\alpha = 0.2125$\AA$^{-1}$ & \\ |
| 1055 |
> |
& & $\alpha = 0.2125$~\AA$^{-1}$ & $\alpha = 0.2125$~\AA$^{-1}$ & \\ |
| 1056 |
|
\midrule |
| 1057 |
|
$\rho$ & (g cm$^{-3}$) & 1.004(4) & 0.995(5) & 0.997 \cite{CRC80}\\ |
| 1058 |
|
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 27(1) & 23(1) & 18.005 \cite{Wagner02} \\ |
| 1059 |
|
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.33(2) & 2.30(5) & 2.299 \cite{Mills73}\\ |
| 1060 |
< |
$n_C$ & & 4.4 & 5.3 & 4.7 \cite{Hura00}\\ |
| 1060 |
> |
$n_C$ & & 4.2 & 5.3 & 4.7 \cite{Hura00}\\ |
| 1061 |
|
$n_H$ & & 3.7 & 4.1 & 3.5 \cite{Soper86}\\ |
| 1062 |
|
$\tau_1$ & (ps) & 5.86(8) & 6.0(1) & 5.7 \cite{Eisenberg69}\\ |
| 1063 |
|
$\tau_2$ & (ps) & 2.45(7) & 2.49(5) & 2.3 \cite{Krynicki66}\\ |
| 1065 |
|
$\tau_D$ & (ps) & 9.1(2) & 10.6(3) & 8.2(4) \cite{Kindt96}\\ |
| 1066 |
|
\bottomrule |
| 1067 |
|
\end{tabular} |
| 1068 |
< |
\label{tab:tredProps} |
| 1068 |
> |
\label{tab:tredProperties} |
| 1069 |
|
\end{table} |
| 1070 |
+ |
We calculated thermodynamic and dynamic properties in the same manner |
| 1071 |
+ |
described in section \ref{sec:ssdrfDamped} for SSD/RF, and the results |
| 1072 |
+ |
are presented in table \ref{tab:tredProperties}. These results show that |
| 1073 |
+ |
TRED improves upon the $\rho$ and $C_p$ properties of SSD/RF with |
| 1074 |
+ |
damped electrostatics while maintaining the excellent dynamics |
| 1075 |
+ |
behavior of both the translational self-diffusivity and the |
| 1076 |
+ |
reorientational time constants, $\tau_1$ and $\tau_2$. The structural |
| 1077 |
+ |
results show some differences between the two models. Despite the |
| 1078 |
+ |
improved peak shape for the first solvation shell (see figure |
| 1079 |
+ |
\ref{fig:tredGofR}), $n_C$ and $n_H$ counts increase because of the |
| 1080 |
+ |
further first minimum distance locations. This results in the |
| 1081 |
+ |
integration extending over a larger water volume. If we integrate to |
| 1082 |
+ |
the first minimum value of the experimental $g_\textrm{OO}(r)$ (3.42 |
| 1083 |
+ |
\AA ) in both the SSD/RF and TRED plots, the $n_C$ values for both are |
| 1084 |
+ |
much closer to experiment (4.7 for SSD/RF and 4.9 for TRED). |
| 1085 |
|
|
| 1086 |
+ |
\begin{figure} |
| 1087 |
+ |
\centering |
| 1088 |
+ |
\includegraphics[width=3in]{./figures/tredDielectric.pdf} |
| 1089 |
+ |
\caption{Contour map of the dielectric constant for TRED as a function |
| 1090 |
+ |
of damping parameter and cutoff radius. The dielectric behavior for TRED |
| 1091 |
+ |
is very similar to SSD/RF (see figure \ref{fig:dielectricMap}D), which |
| 1092 |
+ |
is to be expected due to the similar dipole moment and sticky interaction |
| 1093 |
+ |
strength.} |
| 1094 |
+ |
\label{fig:tredDielectric} |
| 1095 |
+ |
\end{figure} |
| 1096 |
+ |
A comparison of the dielectric behavior was also included at the |
| 1097 |
+ |
bottom of table \ref{tab:tredProperties}. The static dielectric |
| 1098 |
+ |
constant results for SSD/RF and TRED are identical within error. This |
| 1099 |
+ |
is not surprising given the similar dipole moment, similar sticky |
| 1100 |
+ |
interaction strength, and identical applied damping constant. Comparing the static dielectric constant contour map (figure \ref{fig:tredDielectric}) with the dielectric map for SSD/RF |
| 1101 |
+ |
|
| 1102 |
|
\section{Conclusions} |
| 1103 |
|
|
| 1104 |
|
In the above sections, the density maximum and temperature dependence |
| 1105 |
|
of the self-diffusion constant were studied for the SSD water model, |
| 1106 |
|
both with and without the use of reaction field, via a series of {\it |
| 1107 |
|
NPT} and {\it NVE} simulations. The constant pressure simulations |
| 1108 |
< |
showed a density maximum near 260K. In most cases, the calculated |
| 1108 |
> |
showed a density maximum near 260~K. In most cases, the calculated |
| 1109 |
|
densities were significantly lower than the densities obtained from |
| 1110 |
|
other water models (and experiment). Analysis of self-diffusion showed |
| 1111 |
|
SSD to capture the transport properties of water well in both the |
| 1120 |
|
when changing the method of calculating long-range electrostatic |
| 1121 |
|
interactions. |
| 1122 |
|
|
| 1123 |
< |
These simple water models are excellent choices for representing |
| 1124 |
< |
explicit water in large scale simulations of biochemical systems. They |
| 1125 |
< |
are more computationally efficient than the common charge based water |
| 1126 |
< |
models, and, in many cases, exhibit more realistic bulk phase fluid |
| 1127 |
< |
properties. These models are one of the few cases in which maximizing |
| 1128 |
< |
efficiency does not result in a loss in realistic representation. |
| 1123 |
> |
The simple water models investigated here are excellent choices for |
| 1124 |
> |
representing explicit water in large scale simulations of biochemical |
| 1125 |
> |
systems. They are more computationally efficient than the common |
| 1126 |
> |
charge based water models, and, in many cases, exhibit more realistic |
| 1127 |
> |
bulk phase fluid properties. These models are one of the few cases in |
| 1128 |
> |
which maximizing efficiency does not result in a loss in realistic |
| 1129 |
> |
representation. |
