| 5 |
|
molecules of interest. In some cases (such as in the simulation of |
| 6 |
|
phospholipid bilayers), the majority of the calculations that are |
| 7 |
|
performed involve interactions with or between solvent molecules. |
| 8 |
< |
Thus, the motion and behavior of molecules in biochemical simulations |
| 9 |
< |
are highly dependent on the properties of the water model that is |
| 10 |
< |
chosen as the solvent. |
| 8 |
> |
Thus, the motion and behavior of solute molecules in biochemical |
| 9 |
> |
simulations are highly dependent on the properties of the water model |
| 10 |
> |
that is chosen as the solvent. |
| 11 |
|
\begin{figure} |
| 12 |
|
\includegraphics[width=\linewidth]{./figures/waterModels.pdf} |
| 13 |
|
\caption{Partial-charge geometries for the TIP3P, TIP4P, TIP5P, and |
| 53 |
|
Lennard-Jones core. However, since the normal aligned and |
| 54 |
|
anti-aligned geometries favored by point dipoles are poor mimics of |
| 55 |
|
local structure in liquid water, a short ranged ``sticky'' potential |
| 56 |
< |
is also added. The sticky potential directs the molecules to assume |
| 56 |
> |
is also added. This sticky potential directs the molecules to assume |
| 57 |
|
the proper hydrogen bond orientation in the first solvation shell. |
| 58 |
|
|
| 59 |
|
The interaction between two SSD water molecules \emph{i} and \emph{j} |
| 136 |
|
electrostatics.\cite{Ewald21} In applying this water model in a |
| 137 |
|
variety of molecular systems, it would be useful to know its |
| 138 |
|
properties and behavior under the more computationally efficient |
| 139 |
< |
reaction field (RF) technique, the correction techniques discussed in |
| 140 |
< |
the previous chapter, or even a simple |
| 141 |
< |
cutoff.\cite{Onsager36,Fennell06} This chapter addresses these issues |
| 142 |
< |
by looking at the structural and transport behavior of SSD over a |
| 143 |
< |
variety of temperatures with the purpose of utilizing the RF |
| 139 |
> |
reaction field (RF) technique (and later, the {\sc sf} method), the |
| 140 |
> |
correction techniques discussed in the previous chapter, or even a |
| 141 |
> |
simple cutoff.\cite{Onsager36,Fennell06} This chapter addresses these |
| 142 |
> |
issues by looking at the structural and transport behavior of SSD over |
| 143 |
> |
a variety of temperatures with the purpose of utilizing the RF |
| 144 |
|
correction technique. We then suggest modifications to the parameters |
| 145 |
|
that result in more realistic bulk phase behavior. It should be noted |
| 146 |
|
that in a recent publication, some of the original investigators of |
| 147 |
< |
the SSD water model have suggested adjustments to the SSD water model |
| 148 |
< |
to address abnormal density behavior (also observed here), calling the |
| 149 |
< |
corrected model SSD1.\cite{Tan03} In the later sections of this |
| 150 |
< |
chapter, we compare our modified variants of SSD with both the |
| 151 |
< |
original SSD and SSD1 models and discuss how our changes improve the |
| 152 |
< |
depiction of water. |
| 147 |
> |
the SSD water model have suggested adjustments to address abnormal |
| 148 |
> |
density behavior (also observed here), calling the corrected model |
| 149 |
> |
SSD1.\cite{Tan03} In the later sections of this chapter, we compare |
| 150 |
> |
our modified variants of SSD with both the original SSD and SSD1 |
| 151 |
> |
models and discuss how our changes improve the depiction of water. |
| 152 |
|
|
| 153 |
|
\section{Simulation Methods}{\label{sec:waterSimMethods} |
| 154 |
|
|
| 184 |
|
properties on changes in the length of the cutoff |
| 185 |
|
radius.\cite{vanderSpoel98} This behavior makes the reaction field a |
| 186 |
|
less attractive method than the Ewald sum; however, for very large |
| 187 |
< |
systems, the computational benefit of reaction field is significant. |
| 187 |
> |
systems, the computational benefit of a reaction field is significant. |
| 188 |
|
|
| 189 |
|
In contrast to the simulations with a reaction field, we have also |
| 190 |
|
performed a companion set of simulations {\it without} a surrounding |
| 193 |
|
which can be used either with or without an active reaction field. |
| 194 |
|
|
| 195 |
|
To determine the preferred densities of the models, we performed |
| 196 |
< |
simulations in the isobaric-isothermal ({\it NPT}) ensemble. All |
| 196 |
> |
simulations in the isobaric-isothermal ($NPT$) ensemble. All |
| 197 |
|
dynamical properties for these models were then obtained from |
| 198 |
< |
microcanonical ({\it NVE}) simulations done at densities matching the |
| 199 |
< |
{\it NPT} density for a particular target temperature. The constant |
| 198 |
> |
microcanonical ($NVE$) simulations done at densities matching the |
| 199 |
> |
$NPT$ density for a particular target temperature. The constant |
| 200 |
|
pressure simulations were implemented using an integral thermostat and |
| 201 |
|
barostat as outlined by Hoover.\cite{Hoover85,Hoover86} All molecules |
| 202 |
|
were treated as non-linear rigid bodies. Vibrational constraints are |
| 209 |
|
to carry out the integration of the equations of motion in place of |
| 210 |
|
the more prevalent quaternion |
| 211 |
|
method.\cite{Dullweber97,Evans77,Evans77b,Allen87} The reason behind |
| 212 |
< |
this decision was that, in {\it NVE} simulations, the energy drift |
| 212 |
> |
this decision was that, in $NVE$ simulations, the energy drift |
| 213 |
|
when using quaternions was substantially greater than when using the |
| 214 |
< |
{\sc dlm} method (Fig. \ref{fig:timeStepIntegration}). This steady |
| 214 |
> |
{\sc dlm} method (figure \ref{fig:timeStepIntegration}). This steady |
| 215 |
|
drift in the total energy has also been observed in other |
| 216 |
|
studies.\cite{Kol97} |
| 217 |
|
|
| 233 |
|
simulations using the quaternion method and an identical time |
| 234 |
|
step. This additional expense is justified because of the ability to |
| 235 |
|
use time steps that are more that twice as long and still achieve the |
| 236 |
< |
same energy conservation. |
| 236 |
> |
same degree of energy conservation. |
| 237 |
|
|
| 238 |
|
Figure \ref{fig:timeStepIntegration} shows the resulting energy drift |
| 239 |
|
at various time steps for both {\sc dlm} and quaternion |
| 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 2~fs and were also kept at this value for {\it NPT} |
| 255 |
> |
time step. To insure accuracy in our $NVE$ simulations, time steps |
| 256 |
> |
were set at 2~fs and were also kept at this value for $NPT$ |
| 257 |
|
simulations. |
| 258 |
|
|
| 259 |
|
Proton-disordered ice crystals in both the I$_\textrm{h}$ and |
| 261 |
|
simulations. The I$_\textrm{h}$ crystals were formed by first |
| 262 |
|
arranging the centers of mass of the SSD particles into a |
| 263 |
|
``hexagonal'' ice lattice of 1024 particles. Because of the crystal |
| 264 |
< |
structure of I$_\textrm{h}$ ice, the simulation boxes were |
| 264 |
> |
structure of ice I$_\textrm{h}$, the simulation boxes were |
| 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 |
| 271 |
|
pressure of 1~atm 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 1~atm. This procedure resulted in |
| 274 |
< |
structurally stable I$_\textrm{h}$ ice crystals that obey the |
| 274 |
> |
structurally stable ice I$_\textrm{h}$ 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, |
| 276 |
> |
utilized in the making of diamond lattice ice I$_\textrm{c}$ crystals, |
| 277 |
|
with each cubic simulation box consisting of either 512 or 1000 |
| 278 |
|
particles. Only isotropic volume fluctuations were performed under |
| 279 |
|
constant pressure, so the ratio of edge lengths remained constant |
| 282 |
|
\section{SSD Density Behavior} |
| 283 |
|
|
| 284 |
|
Melting studies were performed on the randomized ice crystals using |
| 285 |
< |
the {\it NPT} ensemble. During melting simulations, the melting |
| 285 |
> |
the $NPT$ ensemble. During melting simulations, the melting |
| 286 |
|
transition and the density maximum can both be observed, provided that |
| 287 |
|
the density maximum occurs in the liquid and not the supercooled |
| 288 |
|
regime. It should be noted that the calculated melting temperature |
| 336 |
|
that seen in SSD. Though not included in this plot, it is useful to |
| 337 |
|
note that TIP5P has a density maximum nearly identical to the |
| 338 |
|
experimentally measured temperature (see section |
| 339 |
< |
\ref{sec:t5peDensity}. |
| 339 |
> |
\ref{sec:t5peDensity}). |
| 340 |
|
|
| 341 |
|
Liquid state densities in water have been observed to be dependent on |
| 342 |
|
the cutoff radius ($R_\textrm{c}$), both with and without the use of a |
| 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. |
| 348 |
> |
densities in simulations both with and without a reaction field. |
| 349 |
|
|
| 350 |
|
The key feature to recognize in figure \ref{fig:ssdDense} is the |
| 351 |
|
density scaling of SSD relative to other common models at any given |
| 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 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 |
| 367 |
> |
density maximum location, down to 245~K, is also observed. This is a |
| 368 |
> |
more accurate comparison to the other listed water models, in that no |
| 369 |
> |
long-range corrections were applied in those |
| 370 |
|
simulations.\cite{Baez94,Jorgensen98b} However, even without the |
| 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 |
| 374 |
> |
SSD.\cite{Tan03} Throughout the remainder of this chapter our |
| 375 |
|
reparametrizations of SSD will be compared with their newer SSD1 |
| 376 |
|
model. |
| 377 |
|
|
| 379 |
|
|
| 380 |
|
Accurate dynamical properties of a water model are particularly |
| 381 |
|
important when using the model to study permeation or transport across |
| 382 |
< |
biological membranes. In order to probe transport in bulk water, {\it |
| 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 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 |
| 391 |
< |
in figure \ref{fig:ssdDiffuse}, alongside experimental, SPC/E, and TIP5P |
| 392 |
< |
results.\cite{Gillen72,Holz00,Baez94,Mahoney01} |
| 382 |
> |
biological membranes. In order to probe transport in bulk water, |
| 383 |
> |
$NVE$ simulations were performed at the average densities obtained |
| 384 |
> |
from the $NPT$ simulations at the target temperature. Simulations |
| 385 |
> |
started with randomized velocities and underwent 50~ps of temperature |
| 386 |
> |
scaling and 50~ps of constant energy equilibration before a 200~ps |
| 387 |
> |
data collection run. Diffusion constants were calculated via linear |
| 388 |
> |
fits to the long-time behavior of the mean-square displacement as a |
| 389 |
> |
function of time.\cite{Allen87} The average results from the five sets |
| 390 |
> |
of $NVE$ simulations are displayed in figure \ref{fig:ssdDiffuse}, |
| 391 |
> |
alongside experimental, SPC/E, and TIP5P |
| 392 |
> |
results.\cite{Gillen72,Holz00,Baez94,Mahoney01} |
| 393 |
|
|
| 394 |
|
\begin{figure} |
| 395 |
|
\centering |
| 414 |
|
temperatures and tend to diffuse too rapidly at higher temperatures. |
| 415 |
|
This behavior at higher temperatures is not particularly surprising |
| 416 |
|
since the densities of both TIP5P and SSD are lower than experimental |
| 417 |
< |
water densities at higher temperatures. When calculating the |
| 417 |
> |
water densities at these higher temperatures. When calculating the |
| 418 |
|
diffusion coefficients for SSD at experimental densities (instead of |
| 419 |
< |
the densities from the {\it NPT} simulations), the resulting values |
| 420 |
< |
fall more in line with experiment at these temperatures. |
| 419 |
> |
the densities from the $NPT$ simulations), the resulting values |
| 420 |
> |
fall more in line with experiment at the respective temperatures. |
| 421 |
|
|
| 422 |
|
\section{Structural Changes and Characterization} |
| 423 |
|
|
| 424 |
|
By starting the simulations from the crystalline state, we can get an |
| 425 |
|
estimation of the $T_\textrm{m}$ of the ice structure, and beyond the |
| 426 |
< |
melting point, we study the phase behavior of the liquid. The constant |
| 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 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 |
| 426 |
> |
melting point, we can study the phase behavior of the liquid. The |
| 427 |
> |
constant pressure heat capacity ($C_\textrm{p}$) was monitored to |
| 428 |
> |
locate $T_\textrm{m}$ in each of the simulations. In the melting |
| 429 |
> |
simulations of the 1024 particle ice I$_\textrm{h}$ crystals, a large |
| 430 |
> |
spike in $C_\textrm{p}$ occurs at 245~K, indicating a first order |
| 431 |
> |
phase 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 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 |
| 466 |
|
Additional analysis of the melting process was performed using |
| 467 |
|
two-dimensional structure and dipole angle correlations. Expressions |
| 468 |
|
for these correlations are as follows: |
| 470 |
– |
|
| 469 |
|
\begin{equation} |
| 470 |
< |
g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ , |
| 470 |
> |
g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\left\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\right\rangle\ , |
| 471 |
|
\end{equation} |
| 472 |
|
\begin{equation} |
| 473 |
|
g_{\text{AB}}(r,\cos\omega) = |
| 474 |
< |
\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ , |
| 474 |
> |
\frac{V}{N_\text{A}N_\text{B}}\left\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\right\rangle\ , |
| 475 |
|
\end{equation} |
| 476 |
|
where $\theta$ and $\omega$ refer to the angles shown in figure |
| 477 |
|
\ref{fig:corrAngle}. By binning over both distance and the cosine of the |
| 491 |
|
rapidly. The first solvation shell still shows the strong effect of |
| 492 |
|
the sticky-potential, although it covers a larger area, extending to |
| 493 |
|
include a fraction of aligned dipole peaks within the first solvation |
| 494 |
< |
shell. The latter peaks lose due to thermal motion and as the |
| 494 |
> |
shell. The latter peaks are lost due to thermal motion and as the |
| 495 |
|
competing dipole force overcomes the sticky potential's tight |
| 496 |
|
tetrahedral structuring of the crystal. |
| 497 |
|
|
| 504 |
|
this split character alters to show the leading 4~\AA\ peak dominated |
| 505 |
|
by equatorial anti-parallel dipole orientations. There is also a |
| 506 |
|
tightly bunched group of axially arranged dipoles that most likely |
| 507 |
< |
consist of the smaller fraction of aligned dipole pairs. The trailing |
| 507 |
> |
consist of a smaller fraction of aligned dipole pairs. The trailing |
| 508 |
|
component of the split peak at 5~\AA\ is dominated by aligned dipoles |
| 509 |
|
that assume hydrogen bond arrangements similar to those seen in the |
| 510 |
|
first solvation shell. This evidence indicates that the dipole pair |
| 578 |
|
\begin{figure} |
| 579 |
|
\centering |
| 580 |
|
\includegraphics[width=4.5in]{./figures/newGofRCompare.pdf} |
| 581 |
< |
\caption{ Plots showing the experimental $g(r)$ (Ref. \cite{Hura00}) |
| 581 |
> |
\caption{ Plots showing the experimental $g(r)$ (reference \cite{Hura00}) |
| 582 |
|
with SSD/E and SSD1 without reaction field (top), as well as SSD/RF |
| 583 |
|
and SSD1 with reaction field turned on (bottom). The changes in |
| 584 |
|
parameters have lowered and broadened the first peak of SSD/E and |
| 635 |
|
|
| 636 |
|
While adjusting the location and shape of the first peak of $g(r)$ |
| 637 |
|
improves the densities, these changes alone are insufficient to bring |
| 638 |
< |
the system densities up to the values observed experimentally. To |
| 639 |
< |
further increase the densities, the dipole moments were increased in |
| 640 |
< |
both of our adjusted models. Since SSD is a dipole based model, the |
| 641 |
< |
structure and transport are very sensitive to changes in the dipole |
| 642 |
< |
moment. The original SSD simply used the dipole moment calculated from |
| 643 |
< |
the TIP3P water model, which at 2.35~D is significantly greater than |
| 644 |
< |
the experimental gas phase value of 1.84~D. The larger dipole moment |
| 645 |
< |
is a more realistic value and improves the dielectric properties of |
| 646 |
< |
the fluid. Both theoretical and experimental measurements indicate a |
| 647 |
< |
liquid phase dipole moment ranging from 2.4~D to values as high as |
| 648 |
< |
3.11~D, providing a substantial range of reasonable values for a |
| 649 |
< |
dipole moment.\cite{Sprik91,Gubskaya02,Badyal00,Barriol64} Moderately |
| 638 |
> |
the system densities up to the experimental values. To further |
| 639 |
> |
increase the densities, the dipole moments were increased in both of |
| 640 |
> |
our adjusted models. Since SSD is a dipole based model, the structure |
| 641 |
> |
and transport are very sensitive to changes in the dipole moment. The |
| 642 |
> |
original SSD used the dipole moment of the TIP3P water model, which at |
| 643 |
> |
2.35~D is significantly greater than the experimental gas phase value |
| 644 |
> |
of 1.84~D. The larger dipole moment is a more realistic value and |
| 645 |
> |
improves the dielectric properties of the fluid. Both theoretical and |
| 646 |
> |
experimental measurements indicate a liquid phase dipole moment |
| 647 |
> |
ranging from 2.4~D to values as high as 3.11~D, providing a |
| 648 |
> |
substantial range of reasonable values for a dipole |
| 649 |
> |
moment.\cite{Sprik91,Gubskaya02,Badyal00,Barriol64} Moderately |
| 650 |
|
increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF, |
| 651 |
|
respectively, leads to significant changes in the density and |
| 652 |
|
transport of the water models. |
| 654 |
|
\subsection{Density Behavior} |
| 655 |
|
|
| 656 |
|
In order to demonstrate the benefits of these reparametrizations, we |
| 657 |
< |
performed a series of {\it NPT} and {\it NVE} simulations to probe the |
| 658 |
< |
density and transport properties of the adapted models and compare the |
| 657 |
> |
performed a series of $NPT$ and $NVE$ simulations to probe the density |
| 658 |
> |
and transport properties of the adapted models and compared the |
| 659 |
|
results to the original SSD model. This comparison involved full {\it |
| 660 |
< |
NPT} melting sequences for both SSD/E and SSD/RF, as well as {\it NVE} |
| 660 |
> |
NPT} melting sequences for both SSD/E and SSD/RF, as well as $NVE$ |
| 661 |
|
transport calculations at the calculated self-consistent |
| 662 |
|
densities. Again, the results were obtained from five separate |
| 663 |
|
simulations of 1024 particle systems, and the melting sequences were |
| 664 |
|
started from different ice I$_\textrm{h}$ crystals constructed as |
| 665 |
< |
described previously. Each {\it NPT} simulation was equilibrated for |
| 665 |
> |
described previously. Each $NPT$ simulation was equilibrated for |
| 666 |
|
100~ps before a 200~ps data collection run at each temperature step, |
| 667 |
|
and the final configuration from the previous temperature simulation |
| 668 |
< |
was used as a starting point. All {\it NVE} simulations had the same |
| 668 |
> |
was used as a starting point. All $NVE$ simulations had the same |
| 669 |
|
thermalization, equilibration, and data collection times as stated |
| 670 |
|
previously. |
| 671 |
|
|
| 738 |
|
\includegraphics[width=\linewidth]{./figures/ssdeDiffuse.pdf} |
| 739 |
|
\caption{ The diffusion constants calculated from SSD/E and |
| 740 |
|
SSD1 (both without a reaction field) along with experimental |
| 741 |
< |
results.\cite{Gillen72,Holz00} The {\it NVE} calculations were |
| 742 |
< |
performed at the average densities from the {\it NPT} simulations for |
| 741 |
> |
results.\cite{Gillen72,Holz00} The $NVE$ calculations were |
| 742 |
> |
performed at the average densities from the $NPT$ simulations for |
| 743 |
|
the respective models. SSD/E is slightly more mobile than experiment |
| 744 |
|
at all of the temperatures, but it is closer to experiment at |
| 745 |
|
biologically relevant temperatures than SSD1 without a long-range |
| 754 |
|
of SSD be maintained or improved. Figure \ref{fig:ssdeDiffuse} |
| 755 |
|
compares the temperature dependence of the diffusion constant of SSD/E |
| 756 |
|
to SSD1 without an active reaction field at the densities calculated |
| 757 |
< |
from their respective {\it NPT} simulations at 1 atm. The diffusion |
| 757 |
> |
from their respective $NPT$ simulations at 1 atm. The diffusion |
| 758 |
|
constant for SSD/E is consistently higher than experiment, while SSD1 |
| 759 |
|
remains lower than experiment until relatively high temperatures |
| 760 |
|
(around 360~K). Both models follow the shape of the experimental curve |
| 777 |
|
\includegraphics[width=\linewidth]{./figures/ssdrfDiffuse.pdf} |
| 778 |
|
\caption{ The diffusion constants calculated from SSD/RF and |
| 779 |
|
SSD1 (both with an active reaction field) along with experimental |
| 780 |
< |
results.\cite{Gillen72,Holz00} The {\it NVE} calculations were |
| 781 |
< |
performed at the average densities from the {\it NPT} simulations for |
| 780 |
> |
results.\cite{Gillen72,Holz00} The $NVE$ calculations were |
| 781 |
> |
performed at the average densities from the $NPT$ simulations for |
| 782 |
|
both of the models. SSD/RF captures the self-diffusion of water |
| 783 |
|
throughout most of this temperature range. The increasing diffusion |
| 784 |
|
constants at high temperatures for both models can be attributed to |
| 835 |
|
n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2g_{\textrm{OH}}(r)dr, |
| 836 |
|
\end{equation} |
| 837 |
|
where $\rho$ is the number density of the specified pair interactions, |
| 838 |
< |
$a$ and $b$ are the radial locations of the minima following the first |
| 839 |
< |
peak in $g_\textrm{OO}(r)$ or $g_\textrm{OH}(r)$ respectively. The |
| 840 |
< |
number of hydrogen bonds stays relatively constant across all of the |
| 841 |
< |
models, but the coordination numbers of SSD/E and SSD/RF show an |
| 842 |
< |
improvement over SSD1. This improvement is primarily due to extension |
| 843 |
< |
of the first solvation shell in the new parameter sets. Because $n_H$ |
| 844 |
< |
and $n_C$ are nearly identical in SSD1, it appears that all molecules |
| 845 |
< |
in the first solvation shell are involved in hydrogen bonds. Since |
| 846 |
< |
$n_H$ and $n_C$ differ in the newly parameterized models, the |
| 847 |
< |
orientations in the first solvation shell are a bit more ``fluid''. |
| 848 |
< |
Therefore SSD1 over-structures the first solvation shell and our |
| 849 |
< |
reparametrizations have returned this shell to more realistic |
| 850 |
< |
liquid-like behavior. |
| 838 |
> |
and $a$ and $b$ are the radial locations of the minima following the |
| 839 |
> |
first peak in $g_\textrm{OO}(r)$ or $g_\textrm{OH}(r)$ |
| 840 |
> |
respectively. The number of hydrogen bonds stays relatively constant |
| 841 |
> |
across all of the models, but the coordination numbers of SSD/E and |
| 842 |
> |
SSD/RF show an improvement over SSD1. This improvement is primarily |
| 843 |
> |
due to extension of the first solvation shell in the new parameter |
| 844 |
> |
sets. Because $n_H$ and $n_C$ are nearly identical in SSD1, it |
| 845 |
> |
appears that all molecules in the first solvation shell are involved |
| 846 |
> |
in hydrogen bonds. Since $n_H$ and $n_C$ differ in the newly |
| 847 |
> |
parameterized models, the orientations in the first solvation shell |
| 848 |
> |
are a bit more ``fluid''. Therefore SSD1 over-structures the first |
| 849 |
> |
solvation shell and our reparametrizations have returned this shell to |
| 850 |
> |
more realistic liquid-like behavior. |
| 851 |
|
|
| 852 |
|
The time constants for the orientational autocorrelation functions |
| 853 |
|
are also displayed in Table \ref{tab:liquidProperties}. The dipolar |
| 864 |
|
correlation functions, the orientational relaxation time of the dipole |
| 865 |
|
vector can be calculated from an exponential fit in the long-time |
| 866 |
|
regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these time |
| 867 |
< |
constants were averaged over five detailed {\it NVE} simulations |
| 867 |
> |
constants were averaged over five detailed $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 |
| 870 |
< |
$\tau_2$ was determined from the NMR data in Ref. \cite{Krynicki66} at |
| 870 |
> |
$\tau_2$ was determined from the NMR data in reference \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 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 298~K from the data in Ref. \cite{Eisenberg69}. |
| 875 |
> |
fashion described in reference \cite{Krynicki66}. Similarly, $\tau_1$ was |
| 876 |
> |
recomputed for 298~K from the data in reference \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 |
| 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 |
| 888 |
> |
applying the damped {\sc sf} or {\sc sp} techniques to systems |
| 889 |
|
containing point multipoles. The SSD family of water models is the |
| 890 |
|
perfect test case because of the dipole-dipole (and |
| 891 |
|
charge-dipole/quadrupole) interactions that are present. The {\sc sf} |
| 924 |
|
surprisingly well. The average density shows a modest increase when |
| 925 |
|
using damped electrostatics in place of the reaction field. This comes |
| 926 |
|
about because we neglect the pressure effect due to the surroundings |
| 927 |
< |
outside of the cuttoff, instead relying on screening effects to |
| 927 |
> |
outside of the cutoff, instead relying on screening effects to |
| 928 |
|
neutralize electrostatic interactions at long distances. The $C_p$ |
| 929 |
|
also shows a slight increase, indicating greater fluctuation in the |
| 930 |
|
enthalpy at constant pressure. The only other differences between the |
| 1104 |
|
\ref{fig:dielectricMap}D) highlights the similarities in how these |
| 1105 |
|
models respond to dielectric damping and how the dipolar and monopolar |
| 1106 |
|
electrostatic damping act in an equivalent fashion. Both these |
| 1107 |
< |
dielectric maps span a larger range than the 3, 4, and 5 point-charge |
| 1108 |
< |
water models; however, the SSD/RF range is greater than TRED, |
| 1109 |
< |
indicating that multipoles are a little more sensitive to damping than |
| 1110 |
< |
monopoles. |
| 1107 |
> |
dielectric maps span a larger range than the three-, four-, and |
| 1108 |
> |
five-point water models; however, the SSD/RF range is greater than |
| 1109 |
> |
TRED, indicating that multipoles are more sensitive to damping than |
| 1110 |
> |
monopoles. |
| 1111 |
|
|
| 1112 |
|
The final dielectric comparison comes through the Debye relaxation |
| 1113 |
|
time ($\tau_D$) or the collective dipolar relaxation time when |
| 1114 |
|
assuming a Debye treatment for the dielectric |
| 1115 |
|
relaxation.\cite{Chandra99,Kindt96} This value is calculated through |
| 1116 |
|
equation (\ref{eq:reorientCorr}) applied to the total system dipole |
| 1117 |
< |
moment. The values for both of the models are a slower than the |
| 1117 |
> |
moment. The values for both of the models are slower than the |
| 1118 |
|
experimental relaxation; however, they compare compare very well to |
| 1119 |
|
experiment considering the Debye relaxation times calculated for the |
| 1120 |
|
original SSD (11.95~ps) and the SPC/E (6.95~ps) and TIP3P (6.1~ps) |
| 1126 |
|
|
| 1127 |
|
In the above sections, the density maximum and temperature dependence |
| 1128 |
|
of the self-diffusion constant were studied for the SSD water model, |
| 1129 |
< |
both with and without the use of reaction field, via a series of {\it |
| 1130 |
< |
NPT} and {\it NVE} simulations. The constant pressure simulations |
| 1131 |
< |
showed a density maximum near 260~K. In most cases, the calculated |
| 1132 |
< |
densities were significantly lower than the densities obtained from |
| 1133 |
< |
other water models (and experiment). Analysis of self-diffusion showed |
| 1134 |
< |
SSD to capture the transport properties of water well in both the |
| 1135 |
< |
liquid and supercooled liquid regimes. |
| 1129 |
> |
both with and without the use of reaction field, via a series of $NPT$ |
| 1130 |
> |
and $NVE$ simulations. The constant pressure simulations showed a |
| 1131 |
> |
density maximum near 260~K. In most cases, the calculated densities |
| 1132 |
> |
were significantly lower than the densities obtained from other water |
| 1133 |
> |
models (and experiment). Analysis of self-diffusion showed SSD to |
| 1134 |
> |
capture the transport properties of water well in both the liquid and |
| 1135 |
> |
supercooled liquid regimes. |
| 1136 |
|
|
| 1137 |
|
In order to correct the density behavior, we reparametrized the |
| 1138 |
|
original SSD model for use both with and without a reaction field |
| 1145 |
|
|
| 1146 |
|
We also showed that SSD/RF performs well under the alternative damped |
| 1147 |
|
electrostatic conditions, validating the multipolar damping work in |
| 1148 |
< |
the previous chapter. To improve the modeling of water when {\sc sf}, |
| 1149 |
< |
the TRED water model was developed. This model maintains improves upon |
| 1150 |
< |
the thermodynamic properties of SSD/RF using damped electrostatics |
| 1151 |
< |
while maintaining the exceptional depiction of water dynamics. |
| 1148 |
> |
the previous chapter. To improve the modeling of water with the {\sc |
| 1149 |
> |
sf} technique, the TRED water model was developed. This model improves |
| 1150 |
> |
upon the thermodynamic properties of SSD/RF using damped |
| 1151 |
> |
electrostatics while maintaining the exceptional depiction of water |
| 1152 |
> |
dynamics. |
| 1153 |
|
|
| 1154 |
|
The simple water models investigated here are excellent choices for |
| 1155 |
|
representing water in large scale simulations of biochemical |