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 |
592 |
|
important properties. In this case, it would be ideal to correct the |
593 |
|
densities while maintaining the accurate transport behavior. |
594 |
|
|
595 |
< |
The parameters available for tuning include the $\sigma$ and $\epsilon$ |
596 |
< |
Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky |
597 |
< |
attractive and dipole repulsive terms with their respective |
598 |
< |
cutoffs. To alter the attractive and repulsive terms of the sticky |
599 |
< |
potential independently, it is necessary to separate the terms as |
600 |
< |
follows: |
601 |
< |
\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 |
595 |
> |
The parameters available for tuning include the $\sigma$ and |
596 |
> |
$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the |
597 |
> |
strength of the sticky potential ($\nu_0$), and the sticky attractive |
598 |
> |
and dipole repulsive cubic switching function cutoffs ($r_l$, $r_u$ |
599 |
> |
and $r_l^\prime$, $r_u^\prime$ respectively). The results of the |
600 |
> |
reparameterizations are shown in table \ref{params}. We are calling |
601 |
> |
these reparameterizations the Soft Sticky Dipole / Reaction Field |
602 |
|
(SSD/RF - for use with a reaction field) and Soft Sticky Dipole |
603 |
< |
Enhanced (SSD/E - an attempt to improve the liquid structure in |
603 |
> |
Extended (SSD/E - an attempt to improve the liquid structure in |
604 |
|
simulations without a long-range correction). |
605 |
|
|
606 |
|
\begin{table} |
615 |
|
\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ |
616 |
|
\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\ |
617 |
|
\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ |
618 |
+ |
\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\ |
619 |
|
\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ |
620 |
|
\ \ \ $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\\ |
621 |
|
\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ |
622 |
|
\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ |
623 |
|
\end{tabular} |
793 |
|
\begin{center} |
794 |
|
\epsfxsize=6in |
795 |
|
\epsfbox{ssdeDiffuse.epsi} |
796 |
< |
\caption{Plots of the diffusion constants calculated from SSD/E and SSD1, |
797 |
< |
both without a reaction field, along with experimental results |
798 |
< |
[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations were |
799 |
< |
performed at the average densities observed in the 1 atm NPT |
800 |
< |
simulations for the respective models. SSD/E is slightly more fluid |
801 |
< |
than experiment at all of the temperatures, but it is closer than SSD1 |
802 |
< |
without a long-range correction.} |
796 |
> |
\caption{The diffusion constants calculated from SSD/E and SSD1, |
797 |
> |
both without a reaction field, along with experimental results |
798 |
> |
[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations |
799 |
> |
were performed at the average densities observed in the 1 atm NPT |
800 |
> |
simulations for the respective models. SSD/E is slightly more mobile |
801 |
> |
than experiment at all of the temperatures, but it is closer to |
802 |
> |
experiment at biologically relavent temperatures than SSD1 without a |
803 |
> |
long-range correction.} |
804 |
|
\label{ssdediffuse} |
805 |
|
\end{center} |
806 |
|
\end{figure} |
811 |
|
the densities, it is important that the excellent diffusive behavior |
812 |
|
of SSD be maintained or improved. Figure \ref{ssdediffuse} compares |
813 |
|
the temperature dependence of the diffusion constant of SSD/E to SSD1 |
814 |
< |
without an active reaction field, both at the densities calculated at |
815 |
< |
1 atm and at the experimentally calculated densities for super-cooled |
816 |
< |
and liquid water. The diffusion constant for SSD/E is consistently |
817 |
< |
higher than experiment, while SSD1 remains lower than experiment until |
818 |
< |
relatively high temperatures (greater than 330 K). Both models follow |
819 |
< |
the shape of the experimental curve well below 300 K but tend to |
820 |
< |
diffuse too rapidly at higher temperatures, something that is |
821 |
< |
especially apparent with SSD1. 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 |
826 |
< |
changes made to improve the liquid structure may have had an adverse |
827 |
< |
affect on the density maximum, but they improve the transport behavior |
828 |
< |
of SSD/E relative to SSD1. |
814 |
> |
without an active reaction field at the densities calculated from the |
815 |
> |
NPT simulations at 1 atm. The diffusion constant for SSD/E is |
816 |
> |
consistently higher than experiment, while SSD1 remains lower than |
817 |
> |
experiment until relatively high temperatures (around 360 K). Both |
818 |
> |
models follow the shape of the experimental curve well below 300 K but |
819 |
> |
tend to diffuse too rapidly at higher temperatures, as seen in SSD1's |
820 |
> |
crossing above 360 K. This increasing diffusion relative to the |
821 |
> |
experimental values is caused by the rapidly decreasing system density |
822 |
> |
with increasing temperature. Both SSD1 and SSD/E show this deviation |
823 |
> |
in diffusive behavior, but this trend has different implications on |
824 |
> |
the diffusive behavior of the models. While SSD1 shows more |
825 |
> |
experimentally accurate diffusive behavior in the high temperature |
826 |
> |
regimes, SSD/E shows more accurate behavior in the supercooled and |
827 |
> |
biologically relavent temperature ranges. Thus, the changes made to |
828 |
> |
improve the liquid structure may have had an adverse affect on the |
829 |
> |
density maximum, but they improve the transport behavior of SSD/E |
830 |
> |
relative to SSD1 under the most commonly simulated conditions. |
831 |
|
|
832 |
|
\begin{figure} |
833 |
|
\begin{center} |
834 |
|
\epsfxsize=6in |
835 |
|
\epsfbox{ssdrfDiffuse.epsi} |
836 |
< |
\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1, |
836 |
> |
\caption{The diffusion constants calculated from SSD/RF and SSD1, |
837 |
|
both with an active reaction field, along with experimental results |
838 |
< |
[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations |
838 |
> |
[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations |
839 |
|
were performed at the average densities observed in the 1 atm NPT |
840 |
|
simulations for both of the models. Note how accurately SSD/RF |
841 |
|
simulates the diffusion of water throughout this temperature |
842 |
|
range. The more rapidly increasing diffusion constants at high |
843 |
< |
temperatures for both models is attributed to the significantly lower |
844 |
< |
densities than observed in experiment.} |
843 |
> |
temperatures for both models is attributed to lower calculated |
844 |
> |
densities than those observed in experiment.} |
845 |
|
\label{ssdrfdiffuse} |
846 |
|
\end{center} |
847 |
|
\end{figure} |
849 |
|
In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are |
850 |
|
compared to SSD1 with an active reaction field. Note that SSD/RF |
851 |
|
tracks the experimental results quantitatively, identical within error |
852 |
< |
throughout the temperature range shown and with only a slight |
853 |
< |
increasing trend at higher temperatures. SSD1 tends to diffuse more |
854 |
< |
slowly at low temperatures and deviates to diffuse too rapidly at |
852 |
> |
throughout most of the temperature range shown and exhibiting only a |
853 |
> |
slight increasing trend at higher temperatures. SSD1 tends to diffuse |
854 |
> |
more slowly at low temperatures and deviates to diffuse too rapidly at |
855 |
|
temperatures greater than 330 K. As stated above, this deviation away |
856 |
|
from the ideal trend is due to a rapid decrease in density at higher |
857 |
|
temperatures. SSD/RF does not suffer from this problem as much as SSD1 |
859 |
|
values. These results again emphasize the importance of careful |
860 |
|
reparameterization when using an altered long-range correction. |
861 |
|
|
862 |
+ |
\begin{table} |
863 |
+ |
\begin{center} |
864 |
+ |
\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}.} |
865 |
+ |
\begin{tabular}{ l c c c c c } |
866 |
+ |
\hline \\[-3mm] |
867 |
+ |
\ \ \ \ \ \ & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ |
868 |
+ |
\ & \ SSD/RF \ \ \ & \ Expt. \\ |
869 |
+ |
\hline \\[-3mm] |
870 |
+ |
\ \ \ $\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 \\ |
871 |
+ |
\ \ \ $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 \\ |
872 |
+ |
\ \ \ $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}$ \\ |
873 |
+ |
\ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\ |
874 |
+ |
\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.4$^\text{c}$ \\ |
875 |
+ |
\ \ \ $\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}$ \\ |
876 |
+ |
\ \ \ $\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}$ \\ |
877 |
+ |
\end{tabular} |
878 |
+ |
\label{liquidproperties} |
879 |
+ |
\end{center} |
880 |
+ |
\end{table} |
881 |
+ |
|
882 |
+ |
Table \ref{liquidproperties} gives a synopsis of the liquid state |
883 |
+ |
properties of the water models compared in this study along with the |
884 |
+ |
experimental values for liquid water at ambient conditions. The |
885 |
+ |
coordination number and hydrogen bonds per particle were calculated by |
886 |
+ |
integrating the following relation: |
887 |
+ |
\begin{equation} |
888 |
+ |
4\pi\rho\int_{0}^{a}r^2\text{g}(r)dr, |
889 |
+ |
\end{equation} |
890 |
+ |
where $\rho$ is the number density of pair interactions, $a$ is the |
891 |
+ |
radial location of the minima following the first solvation shell |
892 |
+ |
peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for |
893 |
+ |
calculation of the coordination number or hydrogen bonds per particle |
894 |
+ |
respectively. |
895 |
+ |
|
896 |
+ |
The time constants for the self orientational autocorrelation function |
897 |
+ |
are also displayed in Table \ref{liquidproperties}. The dipolar |
898 |
+ |
orientational time correlation function ($\Gamma_{l}$) is described |
899 |
+ |
by: |
900 |
+ |
\begin{equation} |
901 |
+ |
\Gamma_{l}(t) = \langle P_l[\mathbf{u}_j(0)\cdot\mathbf{u}_j(t)]\rangle, |
902 |
+ |
\end{equation} |
903 |
+ |
where $P_l$ is a Legendre polynomial of order $l$ and $\mathbf{u}_j$ |
904 |
+ |
is the unit vector of the particle dipole.\cite{Rahman71} From these |
905 |
+ |
correlation functions, the orientational relaxation time of the dipole |
906 |
+ |
vector can be calculated from an exponential fit in the long-time |
907 |
+ |
regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculation of these |
908 |
+ |
time constants were averaged from five detailed NVE simulations |
909 |
+ |
performed at the STP density for each of the respective models. |
910 |
+ |
|
911 |
|
\subsection{Additional Observations} |
912 |
|
|
913 |
|
\begin{figure} |