| 1 |
< |
\documentclass[12pt]{ndthesis} |
| 1 |
> |
\documentclass[11pt]{ndthesis} |
| 2 |
|
|
| 3 |
|
% some packages for things like equations and graphics |
| 4 |
+ |
\usepackage[tbtags]{amsmath} |
| 5 |
|
\usepackage{amsmath,bm} |
| 6 |
|
\usepackage{amssymb} |
| 7 |
|
\usepackage{mathrsfs} |
| 8 |
|
\usepackage{tabularx} |
| 9 |
|
\usepackage{graphicx} |
| 10 |
|
\usepackage{booktabs} |
| 11 |
+ |
\usepackage{cite} |
| 12 |
+ |
\usepackage{enumitem} |
| 13 |
|
|
| 14 |
|
\begin{document} |
| 15 |
|
|
| 289 |
|
structural and dynamic properties of water compared the same |
| 290 |
|
properties obtained using the Ewald sum. |
| 291 |
|
|
| 292 |
< |
\section{Simple Forms for Pairwise Electrostatics} |
| 292 |
> |
\section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation} |
| 293 |
|
|
| 294 |
|
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 295 |
|
al.} are constructed using two different (and separable) computational |
| 296 |
|
tricks: |
| 297 |
|
|
| 298 |
< |
\begin{enumerate} |
| 298 |
> |
\begin{enumerate}[itemsep=0pt] |
| 299 |
|
\item shifting through the use of image charges, and |
| 300 |
|
\item damping the electrostatic interaction. |
| 301 |
|
\end{enumerate} |
| 302 |
< |
Wolf \textit{et al.} treated the |
| 303 |
< |
development of their summation method as a progressive application of |
| 304 |
< |
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 305 |
< |
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 306 |
< |
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
| 307 |
< |
both techniques. It is possible, however, to separate these |
| 308 |
< |
tricks and study their effects independently. |
| 302 |
> |
Wolf \textit{et al.} treated the development of their summation method |
| 303 |
> |
as a progressive application of these techniques,\cite{Wolf99} while |
| 304 |
> |
Zahn \textit{et al.} founded their damped Coulomb modification |
| 305 |
> |
(Eq. (\ref{eq:ZahnPot})) on the post-limit forces |
| 306 |
> |
(Eq. (\ref{eq:WolfForces})) which were derived using both techniques. |
| 307 |
> |
It is possible, however, to separate these tricks and study their |
| 308 |
> |
effects independently. |
| 309 |
|
|
| 310 |
|
Starting with the original observation that the effective range of the |
| 311 |
|
electrostatic interaction in condensed phases is considerably less |
| 646 |
|
simulations (i.e. from liquids of neutral molecules to ionic |
| 647 |
|
crystals), so the systems studied were: |
| 648 |
|
|
| 649 |
< |
\begin{enumerate} |
| 649 |
> |
\begin{enumerate}[itemsep=0pt] |
| 650 |
|
\item liquid water (SPC/E),\cite{Berendsen87} |
| 651 |
|
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 652 |
|
\item NaCl crystals, |
| 691 |
|
We compared the following alternative summation methods with results |
| 692 |
|
from the reference method ({\sc spme}): |
| 693 |
|
|
| 694 |
< |
\begin{enumerate} |
| 694 |
> |
\begin{enumerate}[itemsep=0pt] |
| 695 |
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 696 |
|
and 0.3\AA$^{-1}$, |
| 697 |
|
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 779 |
|
the complementary error function is required). |
| 780 |
|
|
| 781 |
|
The reaction field results illustrates some of that method's |
| 782 |
< |
limitations, primarily that it was developed for use in homogenous |
| 782 |
> |
limitations, primarily that it was developed for use in homogeneous |
| 783 |
|
systems; although it does provide results that are an improvement over |
| 784 |
|
those from an unmodified cutoff. |
| 785 |
|
|
| 904 |
|
all do equivalently well at capturing the direction of both the force |
| 905 |
|
and torque vectors. Using the electrostatic damping improves the |
| 906 |
|
angular behavior significantly for the {\sc sp} and moderately for the |
| 907 |
< |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 907 |
> |
{\sc sf} methods. Over-damping is detrimental to both methods. Again |
| 908 |
|
it is important to recognize that the force vectors cover all |
| 909 |
|
particles in all seven systems, while torque vectors are only |
| 910 |
|
available for neutral molecular groups. Damping is more beneficial to |
| 911 |
|
charged bodies, and this observation is investigated further in |
| 912 |
< |
section \ref{IndividualResults}. |
| 912 |
> |
section \ref{sec:IndividualResults}. |
| 913 |
|
|
| 914 |
|
Although not discussed previously, group based cutoffs can be applied |
| 915 |
|
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 935 |
|
|
| 936 |
|
\footnotesize |
| 937 |
|
\begin{center} |
| 938 |
< |
\begin{tabular}{@{} ccrrrrrrrr @{}} \\ |
| 938 |
> |
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 939 |
|
\toprule |
| 940 |
|
\toprule |
| 938 |
– |
|
| 941 |
|
& &\multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted |
| 942 |
|
Force} \\ |
| 943 |
|
\cmidrule(lr){3-6} \cmidrule(l){7-10} $R_\textrm{c}$ & Groups & |
| 973 |
|
increases, something that is more obvious with group-based cutoffs. |
| 974 |
|
The complimentary error function inserted into the potential weakens |
| 975 |
|
the electrostatic interaction as the value of $\alpha$ is increased. |
| 976 |
< |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 976 |
> |
However, at larger values of $\alpha$, it is possible to over-damp the |
| 977 |
|
electrostatic interaction and to remove it completely. Kast |
| 978 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 979 |
|
values for these types of electrostatic summation methods by fitting |
| 1010 |
|
investigated. In all of the individual results table, the method |
| 1011 |
|
abbreviations are as follows: |
| 1012 |
|
|
| 1013 |
< |
\begin{itemize} |
| 1013 |
> |
\begin{itemize}[itemsep=0pt] |
| 1014 |
|
\item PC = Pure Cutoff, |
| 1015 |
|
\item SP = Shifted Potential, |
| 1016 |
|
\item SF = Shifted Force, |
| 1028 |
|
|
| 1029 |
|
\footnotesize |
| 1030 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1029 |
– |
\\ |
| 1031 |
|
\toprule |
| 1032 |
|
\toprule |
| 1033 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1084 |
|
|
| 1085 |
|
\footnotesize |
| 1086 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1086 |
– |
\\ |
| 1087 |
|
\toprule |
| 1088 |
|
\toprule |
| 1089 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1121 |
|
agreement with {\sc spme} in both energetic and dynamic behavior when |
| 1122 |
|
using the {\sc sf} method with and without damping. The {\sc sp} |
| 1123 |
|
method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly |
| 1124 |
< |
with cutoff radii greater than 12\AA. Overdamping the electrostatics |
| 1124 |
> |
with cutoff radii greater than 12\AA. Over-damping the electrostatics |
| 1125 |
|
reduces the agreement between both these methods and {\sc spme}. |
| 1126 |
|
|
| 1127 |
|
The pure cutoff ({\sc pc}) method performs poorly, again mirroring the |
| 1154 |
|
no damping and only modest improvement for the recommended conditions |
| 1155 |
|
($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA). The |
| 1156 |
|
{\sc sf} method shows modest narrowing across all damping and cutoff |
| 1157 |
< |
ranges of interest. When overdamping these methods, group cutoffs and |
| 1157 |
> |
ranges of interest. When over-damping these methods, group cutoffs and |
| 1158 |
|
the switching function do not improve the force and torque |
| 1159 |
|
directionalities. |
| 1160 |
|
|
| 1180 |
|
|
| 1181 |
|
\footnotesize |
| 1182 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1183 |
– |
\\ |
| 1183 |
|
\toprule |
| 1184 |
|
\toprule |
| 1185 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1235 |
|
|
| 1236 |
|
\footnotesize |
| 1237 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1239 |
– |
\\ |
| 1238 |
|
\toprule |
| 1239 |
|
\toprule |
| 1240 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque |
| 1304 |
|
|
| 1305 |
|
\footnotesize |
| 1306 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1309 |
– |
\\ |
| 1307 |
|
\toprule |
| 1308 |
|
\toprule |
| 1309 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1343 |
|
|
| 1344 |
|
\footnotesize |
| 1345 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1349 |
– |
\\ |
| 1346 |
|
\toprule |
| 1347 |
|
\toprule |
| 1348 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1394 |
|
|
| 1395 |
|
\footnotesize |
| 1396 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1401 |
– |
\\ |
| 1397 |
|
\toprule |
| 1398 |
|
\toprule |
| 1399 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1434 |
|
|
| 1435 |
|
\footnotesize |
| 1436 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1442 |
– |
\\ |
| 1437 |
|
\toprule |
| 1438 |
|
\toprule |
| 1439 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1492 |
|
|
| 1493 |
|
\footnotesize |
| 1494 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1501 |
– |
\\ |
| 1495 |
|
\toprule |
| 1496 |
|
\toprule |
| 1497 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1548 |
|
|
| 1549 |
|
\footnotesize |
| 1550 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1558 |
– |
\\ |
| 1551 |
|
\toprule |
| 1552 |
|
\toprule |
| 1553 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1611 |
|
|
| 1612 |
|
\footnotesize |
| 1613 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1622 |
– |
\\ |
| 1614 |
|
\toprule |
| 1615 |
|
\toprule |
| 1616 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1667 |
|
|
| 1668 |
|
\footnotesize |
| 1669 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1679 |
– |
\\ |
| 1670 |
|
\toprule |
| 1671 |
|
\toprule |
| 1672 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1731 |
|
|
| 1732 |
|
\footnotesize |
| 1733 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1744 |
– |
\\ |
| 1734 |
|
\toprule |
| 1735 |
|
\toprule |
| 1736 |
|
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1787 |
|
|
| 1788 |
|
\footnotesize |
| 1789 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1801 |
– |
\\ |
| 1790 |
|
\toprule |
| 1791 |
|
\toprule |
| 1792 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1821 |
|
|
| 1822 |
|
This system does not appear to show any significant deviations from |
| 1823 |
|
the previously observed results. The {\sc sp} and {\sc sf} methods |
| 1824 |
< |
have aggrements similar to those observed in section |
| 1824 |
> |
have agreements similar to those observed in section |
| 1825 |
|
\ref{sec:WaterResults}. The only significant difference is the |
| 1826 |
|
improvement in the configuration energy differences for the {\sc rf} |
| 1827 |
|
method. This is surprising in that we are introducing an inhomogeneity |
| 1941 |
|
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1942 |
|
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1943 |
|
electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement |
| 1944 |
< |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1944 |
> |
with the power spectrum obtained using the Ewald sum. Over-damping can |
| 1945 |
|
result in underestimates of frequencies of the long-wavelength |
| 1946 |
|
motions.} |
| 1947 |
|
\label{fig:dampInc} |
| 1948 |
|
\end{figure} |
| 1949 |
|
|
| 1950 |
< |
\section{Synopsis of the Pairwise Method Evaluation}\label{sec:PairwiseSynopsis} |
| 1950 |
> |
\section{An Application: TIP5P-E Water}\label{sec:t5peApplied} |
| 1951 |
> |
|
| 1952 |
> |
The above sections focused on the energetics and dynamics of a variety |
| 1953 |
> |
of systems when utilizing the {\sc sp} and {\sc sf} pairwise |
| 1954 |
> |
techniques. A unitary correlation with results obtained using the |
| 1955 |
> |
Ewald summation should result in a successful reproduction of both the |
| 1956 |
> |
static and dynamic properties of any selected system. To test this, |
| 1957 |
> |
we decided to calculate a series of properties for the TIP5P-E water |
| 1958 |
> |
model when using the {\sc sf} technique. |
| 1959 |
> |
|
| 1960 |
> |
The TIP5P-E water model is a variant of Mahoney and Jorgensen's |
| 1961 |
> |
five-point transferable intermolecular potential (TIP5P) model for |
| 1962 |
> |
water.\cite{Mahoney00} TIP5P was developed to reproduce the density |
| 1963 |
> |
maximum anomaly present in liquid water near 4$^\circ$C. As with many |
| 1964 |
> |
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC, |
| 1965 |
> |
and SPC/E), TIP5P was parametrized using a simple cutoff with no |
| 1966 |
> |
long-range electrostatic |
| 1967 |
> |
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
| 1968 |
> |
Without this correction, the pressure term on the central particle |
| 1969 |
> |
from the surroundings is missing. Because they expand to compensate |
| 1970 |
> |
for this added pressure term when this correction is included, systems |
| 1971 |
> |
composed of these particles tend to underpredict the density of water |
| 1972 |
> |
under standard conditions. When using any form of long-range |
| 1973 |
> |
electrostatic correction, it has become common practice to develop or |
| 1974 |
> |
utilize a reparametrized water model that corrects for this |
| 1975 |
> |
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 1976 |
> |
this practice and was optimized specifically for use with the Ewald |
| 1977 |
> |
summation.\cite{Rick04} In his publication, Rick preserved the |
| 1978 |
> |
geometry and point charge magnitudes in TIP5P and focused on altering |
| 1979 |
> |
the Lennard-Jones parameters to correct the density at |
| 1980 |
> |
298K.\cite{Rick04} With the density corrected, he compared common |
| 1981 |
> |
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1982 |
> |
9\AA\ cutoff. |
| 1983 |
> |
|
| 1984 |
> |
In the following sections, we compared these same water properties |
| 1985 |
> |
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1986 |
> |
{\sc sf} technique. In the above evaluation of the pairwise |
| 1987 |
> |
techniques, we observed some flexibility in the choice of parameters. |
| 1988 |
> |
Because of this, the following comparisons include the {\sc sf} |
| 1989 |
> |
technique with a 12\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 1990 |
> |
0.2\AA$^{-1}$, as well as a 9\AA\ cutoff with an $\alpha$ = |
| 1991 |
> |
0.2\AA$^{-1}$. |
| 1992 |
> |
|
| 1993 |
> |
\subsection{Density}\label{sec:t5peDensity} |
| 1994 |
> |
|
| 1995 |
> |
As stated previously, the property that prompted the development of |
| 1996 |
> |
TIP5P-E was the density at 1 atm. The density depends upon the |
| 1997 |
> |
internal pressure of the system in the $NPT$ ensemble, and the |
| 1998 |
> |
calculation of the pressure includes a components from both the |
| 1999 |
> |
kinetic energy and the virial. More specifically, the instantaneous |
| 2000 |
> |
molecular pressure ($p(t)$) is given by |
| 2001 |
> |
\begin{equation} |
| 2002 |
> |
p(t) = \frac{1}{\textrm{d}V}\sum_\mu |
| 2003 |
> |
\left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} |
| 2004 |
> |
+ \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right], |
| 2005 |
> |
\label{eq:MolecularPressure} |
| 2006 |
> |
\end{equation} |
| 2007 |
> |
where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of |
| 2008 |
> |
molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass |
| 2009 |
> |
($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on |
| 2010 |
> |
atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the |
| 2011 |
> |
right term in the brackets of eq. \ref{eq:MolecularPressure}) is |
| 2012 |
> |
directly dependent on the interatomic forces. Since the {\sc sp} |
| 2013 |
> |
method does not modify the forces (see |
| 2014 |
> |
sec. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will |
| 2015 |
> |
be identical to that obtained without an electrostatic correction. |
| 2016 |
> |
The {\sc sf} method does alter the virial component and, by way of the |
| 2017 |
> |
modified pressures, should provide densities more in line with those |
| 2018 |
> |
obtained using the Ewald summation. |
| 2019 |
> |
|
| 2020 |
> |
To compare densities, $NPT$ simulations were performed with the same |
| 2021 |
> |
temperatures as those selected by Rick in his Ewald summation |
| 2022 |
> |
simulations.\cite{Rick04} In order to improve statistics around the |
| 2023 |
> |
density maximum, 3ns trajectories were accumulated at 0, 12.5, and |
| 2024 |
> |
25$^\circ$C, while 2ns trajectories were obtained at all other |
| 2025 |
> |
temperatures. The average densities were calculated from the later |
| 2026 |
> |
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 2027 |
> |
method for accumulating statistics, these sequences were spliced into |
| 2028 |
> |
200 segments to calculate the average density and standard deviation |
| 2029 |
> |
at each temperature.\cite{Mahoney00} |
| 2030 |
> |
|
| 2031 |
> |
\begin{figure} |
| 2032 |
> |
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 2033 |
> |
\caption{Density versus temperature for the TIP5P-E water model when |
| 2034 |
> |
using the Ewald summation (Ref. \cite{Rick04}) and the {\sc sf} method |
| 2035 |
> |
with various parameters. The pressure term from the image-charge shell |
| 2036 |
> |
is larger than that provided by the reciprocal-space portion of the |
| 2037 |
> |
Ewald summation, leading to slightly lower densities. This effect is |
| 2038 |
> |
more visible with the 9\AA\ cutoff, where the image charges exert a |
| 2039 |
> |
greater force on the central particle. The error bars for the {\sc sf} |
| 2040 |
> |
methods show plus or minus the standard deviation of the density |
| 2041 |
> |
measurement at each temperature.} |
| 2042 |
> |
\label{fig:t5peDensities} |
| 2043 |
> |
\end{figure} |
| 2044 |
|
|
| 2045 |
+ |
Figure \ref{fig:t5peDensities} shows the densities calculated for |
| 2046 |
+ |
TIP5P-E using differing electrostatic corrections overlaid on the |
| 2047 |
+ |
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 2048 |
+ |
technique are close to, though typically lower than, those calculated |
| 2049 |
+ |
while using the Ewald summation. These slightly reduced densities |
| 2050 |
+ |
indicate that the pressure component from the image charges at |
| 2051 |
+ |
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 2052 |
+ |
portion of the Ewald summation. Bringing the image charges closer to |
| 2053 |
+ |
the central particle by choosing a 9\AA\ R$_\textrm{c}$ (rather than |
| 2054 |
+ |
the preferred 12\AA\ R$_\textrm{c}$) increases the strength of their |
| 2055 |
+ |
interactions, resulting in a further reduction of the densities. |
| 2056 |
+ |
|
| 2057 |
+ |
Because the strength of the image charge interactions has a noticable |
| 2058 |
+ |
effect on the density, we would expect the use of electrostatic |
| 2059 |
+ |
damping to also play a role in these calculations. Larger values of |
| 2060 |
+ |
$\alpha$ weaken the pair-interactions; and since electrostatic damping |
| 2061 |
+ |
is distance-dependent, force components from the image charges will be |
| 2062 |
+ |
reduced more than those from particles close the the central |
| 2063 |
+ |
charge. This effect is visible in figure \ref{fig:t5peDensities} with |
| 2064 |
+ |
the damped {\sc sf} sums showing slightly higher densities; however, |
| 2065 |
+ |
it is apparent that the choice of cutoff radius plays a much more |
| 2066 |
+ |
important role in the resulting densities. |
| 2067 |
+ |
|
| 2068 |
+ |
As a final note, all of the above density calculations were performed |
| 2069 |
+ |
with systems of 512 water molecules. Rick observed a system sized |
| 2070 |
+ |
dependence of the computed densities when using the Ewald summation, |
| 2071 |
+ |
most likely due to his tying of the convergence parameter to the box |
| 2072 |
+ |
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
| 2073 |
+ |
calculated densities were on average 0.002 to 0.003 g/cm$^3$ lower. A |
| 2074 |
+ |
system size of 256 molecules would force the use of a shorter |
| 2075 |
+ |
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also |
| 2076 |
+ |
lower the densities. Moving to larger systems, as long as the |
| 2077 |
+ |
R$_\textrm{c}$ remains at a fixed value, we would expect the densities |
| 2078 |
+ |
to remain constant. |
| 2079 |
+ |
|
| 2080 |
+ |
\subsection{Liquid Structure}\label{sec:t5peLiqStructure} |
| 2081 |
+ |
|
| 2082 |
+ |
A common function considered when developing and comparing water |
| 2083 |
+ |
models is the oxygen-oxygen radial distribution function |
| 2084 |
+ |
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of |
| 2085 |
+ |
finding a pair of oxygen atoms some distance ($r$) apart relative to a |
| 2086 |
+ |
random distribution at the same density.\cite{Allen87} It is |
| 2087 |
+ |
calculated via |
| 2088 |
+ |
\begin{equation} |
| 2089 |
+ |
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i} |
| 2090 |
+ |
\delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle, |
| 2091 |
+ |
\label{eq:GOOofR} |
| 2092 |
+ |
\end{equation} |
| 2093 |
+ |
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen |
| 2094 |
+ |
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or |
| 2095 |
+ |
neutron scattering experiments through the oxygen-oxygen structure |
| 2096 |
+ |
factor ($S_\textrm{OO}(k)$) by the following relationship: |
| 2097 |
+ |
\begin{equation} |
| 2098 |
+ |
S_\textrm{OO}(k) = 1 + 4\pi\rho |
| 2099 |
+ |
\int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r. |
| 2100 |
+ |
\label{eq:SOOofK} |
| 2101 |
+ |
\end{equation} |
| 2102 |
+ |
Thus, $S_\textrm{OO}(k)$ is related to the Fourier-Laplace transform |
| 2103 |
+ |
of $g_\textrm{OO}(r)$. |
| 2104 |
+ |
|
| 2105 |
+ |
The expermentally determined $g_\textrm{OO}(r)$ for liquid water has |
| 2106 |
+ |
been compared in great detail with the various common water models, |
| 2107 |
+ |
and TIP5P was found to be in better agreement than other rigid, |
| 2108 |
+ |
non-polarizable models.\cite{Sorenson00} This excellent agreement with |
| 2109 |
+ |
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To |
| 2110 |
+ |
check whether the choice of using the Ewald summation or the {\sc sf} |
| 2111 |
+ |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K |
| 2112 |
+ |
and 1atm were determined for the systems compared in the previous |
| 2113 |
+ |
section. |
| 2114 |
+ |
|
| 2115 |
+ |
\begin{figure} |
| 2116 |
+ |
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 2117 |
+ |
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298K and |
| 2118 |
+ |
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc |
| 2119 |
+ |
sf} technique with varying parameters. Even with the reduced densities |
| 2120 |
+ |
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially |
| 2121 |
+ |
identical.} |
| 2122 |
+ |
\label{fig:t5peGofRs} |
| 2123 |
+ |
\end{figure} |
| 2124 |
+ |
|
| 2125 |
+ |
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc |
| 2126 |
+ |
sf} technique with a various parameters are overlaid on the |
| 2127 |
+ |
$g_\textrm{OO}(r)$ while using the Ewald summation. The differences in |
| 2128 |
+ |
density do not appear to have any effect on the liquid structure as |
| 2129 |
+ |
the $g_\textrm{OO}(r)$s are indistinquishable. These results indicate |
| 2130 |
+ |
that the $g_\textrm{OO}(r)$ is insensitive to the choice of |
| 2131 |
+ |
electrostatic correction. |
| 2132 |
+ |
|
| 2133 |
+ |
\subsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 2134 |
+ |
|
| 2135 |
+ |
In addition to the density, there are a variety of thermodynamic |
| 2136 |
+ |
quantities that can be calculated for water and compared directly to |
| 2137 |
+ |
experimental values. Some of these additional quatities include the |
| 2138 |
+ |
latent heat of vaporization ($\Delta H_\textrm{vap}$), the constant |
| 2139 |
+ |
pressure heat capacity ($C_p$), the isothermal compressibility |
| 2140 |
+ |
($\kappa_T$), the thermal expansivity ($\alpha_p$), and the static |
| 2141 |
+ |
dielectric constant ($\epsilon$). All of these properties were |
| 2142 |
+ |
calculated for TIP5P-E with the Ewald summation, so they provide a |
| 2143 |
+ |
good set for comparisons involving the {\sc sf} technique. |
| 2144 |
+ |
|
| 2145 |
+ |
The $\Delta H_\textrm{vap}$ is the enthalpy change required to |
| 2146 |
+ |
transform one mol of substance from the liquid phase to the gas |
| 2147 |
+ |
phase.\cite{Berry00} In molecular simulations, this quantity can be |
| 2148 |
+ |
determined via |
| 2149 |
+ |
\begin{equation} |
| 2150 |
+ |
\begin{split} |
| 2151 |
+ |
\Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq.} \\ |
| 2152 |
+ |
&= E_\textrm{gas} - E_\textrm{liq.} |
| 2153 |
+ |
+ p(V_\textrm{gas} - V_\textrm{liq.}) \\ |
| 2154 |
+ |
&\approx -\frac{\langle U_\textrm{liq.}\rangle}{N}+ RT, |
| 2155 |
+ |
\end{split} |
| 2156 |
+ |
\label{eq:DeltaHVap} |
| 2157 |
+ |
\end{equation} |
| 2158 |
+ |
where $E$ is the total energy, $U$ is the potential energy, $p$ is the |
| 2159 |
+ |
pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is |
| 2160 |
+ |
the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As |
| 2161 |
+ |
seen in the last line of equation \ref{eq:DeltaHVap}, we can |
| 2162 |
+ |
approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas |
| 2163 |
+ |
state. This allows us to cancel the kinetic energy terms, leaving only |
| 2164 |
+ |
the $U_\textrm{liq.}$ term. Additionally, the $pV$ term for the gas is |
| 2165 |
+ |
several orders of magnitude larger than that of the liquid, so we can |
| 2166 |
+ |
neglect the liquid $pV$ term. |
| 2167 |
+ |
|
| 2168 |
+ |
The remaining thermodynamic properties can all be calculated from |
| 2169 |
+ |
fluctuations of the enthalpy, volume, and system dipole |
| 2170 |
+ |
moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the |
| 2171 |
+ |
enthalpy in constant pressure simulations via |
| 2172 |
+ |
\begin{equation} |
| 2173 |
+ |
\begin{split} |
| 2174 |
+ |
C_p = \left(\frac{\partial H}{\partial T}\right)_{N,p} |
| 2175 |
+ |
= \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2}, |
| 2176 |
+ |
\end{split} |
| 2177 |
+ |
\label{eq:Cp} |
| 2178 |
+ |
\end{equation} |
| 2179 |
+ |
where $k_B$ is Boltzmann's constant. $\kappa_T$ can be calculated via |
| 2180 |
+ |
\begin{equation} |
| 2181 |
+ |
\begin{split} |
| 2182 |
+ |
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T} |
| 2183 |
+ |
= \frac{(\langle V^2\rangle_{N,P,T} - \langle V\rangle^2_{N,P,T})} |
| 2184 |
+ |
{k_BT\langle V\rangle_{N,P,T}}, |
| 2185 |
+ |
\end{split} |
| 2186 |
+ |
\label{eq:kappa} |
| 2187 |
+ |
\end{equation} |
| 2188 |
+ |
and $\alpha_p$ can be calculated via |
| 2189 |
+ |
\begin{equation} |
| 2190 |
+ |
\begin{split} |
| 2191 |
+ |
\alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P} |
| 2192 |
+ |
= \frac{(\langle VH\rangle_{N,P,T} |
| 2193 |
+ |
- \langle V\rangle_{N,P,T}\langle H\rangle_{N,P,T})} |
| 2194 |
+ |
{k_BT^2\langle V\rangle_{N,P,T}}. |
| 2195 |
+ |
\end{split} |
| 2196 |
+ |
\label{eq:alpha} |
| 2197 |
+ |
\end{equation} |
| 2198 |
+ |
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can |
| 2199 |
+ |
be calculated for systems of non-polarizable substances via |
| 2200 |
+ |
\begin{equation} |
| 2201 |
+ |
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
| 2202 |
+ |
\label{eq:staticDielectric} |
| 2203 |
+ |
\end{equation} |
| 2204 |
+ |
where $\epsilon_0$ is the permittivity of free space and $\langle |
| 2205 |
+ |
M^2\rangle$ is the fluctuation of the system dipole |
| 2206 |
+ |
moment.\cite{Allen87} The numerator in the fractional term in equation |
| 2207 |
+ |
\ref{eq:staticDielectric} is the fluctuation of the simulation-box |
| 2208 |
+ |
dipole moment, identical to the quantity calculated in the |
| 2209 |
+ |
finite-system Kirkwood $g$ factor ($G_k$): |
| 2210 |
+ |
\begin{equation} |
| 2211 |
+ |
G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
| 2212 |
+ |
\label{eq:KirkwoodFactor} |
| 2213 |
+ |
\end{equation} |
| 2214 |
+ |
where $\mu$ is the dipole moment of a single molecule of the |
| 2215 |
+ |
homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The |
| 2216 |
+ |
fluctuation term in both equation \ref{eq:staticDielectric} and |
| 2217 |
+ |
\ref{eq:KirkwoodFactor} is calculated as follows, |
| 2218 |
+ |
\begin{equation} |
| 2219 |
+ |
\begin{split} |
| 2220 |
+ |
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
| 2221 |
+ |
- \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
| 2222 |
+ |
&= \langle M_x^2+M_y^2+M_z^2\rangle |
| 2223 |
+ |
- (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
| 2224 |
+ |
+ \langle M_x\rangle^2). |
| 2225 |
+ |
\end{split} |
| 2226 |
+ |
\label{eq:fluctBoxDipole} |
| 2227 |
+ |
\end{equation} |
| 2228 |
+ |
This fluctuation term can be accumulated during the simulation; |
| 2229 |
+ |
however, it converges rather slowly, thus requiring multi-nanosecond |
| 2230 |
+ |
simulation times.\cite{Horn04} In the case of tin-foil boundary |
| 2231 |
+ |
conditions, the dielectric/surface term of equation \ref{eq:EwaldSum} |
| 2232 |
+ |
is equal to zero. Since the {\sc sf} method also lacks this |
| 2233 |
+ |
dielectric/surface term, equation \ref{eq:staticDielectric} is still |
| 2234 |
+ |
valid for determining static dielectric constants. |
| 2235 |
+ |
|
| 2236 |
+ |
All of the above properties were calculated from the same trajectories |
| 2237 |
+ |
used to determine the densities in section \ref{sec:t5peDensity} |
| 2238 |
+ |
except for the static dielectric constants. The $\epsilon$ values were |
| 2239 |
+ |
accumulated from 2ns $NVE$ ensemble trajectories with system densities |
| 2240 |
+ |
fixed at the average values from the $NPT$ simulations at each of the |
| 2241 |
+ |
temperatures. The resulting values are displayed in figure |
| 2242 |
+ |
\ref{fig:t5peThermo}. |
| 2243 |
+ |
\begin{figure} |
| 2244 |
+ |
\centering |
| 2245 |
+ |
\includegraphics[width=5.5in]{./figures/t5peThermo.pdf} |
| 2246 |
+ |
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation |
| 2247 |
+ |
and the {\sc sf} techniques along with the experimental values. Units |
| 2248 |
+ |
for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$, |
| 2249 |
+ |
cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$, |
| 2250 |
+ |
and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from |
| 2251 |
+ |
reference \cite{Rick04}. Experimental values for $\Delta |
| 2252 |
+ |
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference |
| 2253 |
+ |
\cite{Kell75}. Experimental values for $C_p$ are from reference |
| 2254 |
+ |
\cite{Wagner02}. Experimental values for $\epsilon$ are from reference |
| 2255 |
+ |
\cite{Malmberg56}.} |
| 2256 |
+ |
\label{fig:t5peThermo} |
| 2257 |
+ |
\end{figure} |
| 2258 |
+ |
|
| 2259 |
+ |
As observed for the density in section \ref{sec:t5peDensity}, the |
| 2260 |
+ |
property trends with temperature seen when using the Ewald summation |
| 2261 |
+ |
are reproduced with the {\sc sf} technique. Differences include the |
| 2262 |
+ |
calculated values of $\Delta H_\textrm{vap}$ underpredicting the Ewald |
| 2263 |
+ |
values. This is to be expected due to the direct weakening of the |
| 2264 |
+ |
electrostatic interaction through forced neutralization in {\sc |
| 2265 |
+ |
sf}. This results in an increase of the intermolecular potential |
| 2266 |
+ |
producing lower values from equation \ref{eq:DeltaHVap}. The slopes of |
| 2267 |
+ |
these values with temperature are similar to that seen using the Ewald |
| 2268 |
+ |
summation; however, they are both steeper than the experimental trend, |
| 2269 |
+ |
indirectly resulting in the inflated $C_p$ values at all temperatures. |
| 2270 |
+ |
|
| 2271 |
+ |
Above the supercooled regim\'{e}, $C_p$, $\kappa_T$, and $\alpha_p$ |
| 2272 |
+ |
values all overlap within error. As indicated for the $\Delta |
| 2273 |
+ |
H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph, |
| 2274 |
+ |
the deviations between experiment and simulation in this region are |
| 2275 |
+ |
not the fault of the electrostatic summation methods but are due to |
| 2276 |
+ |
the TIP5P class model itself. Like most rigid, non-polarizable, |
| 2277 |
+ |
point-charge water models, the density decreases with temperature at a |
| 2278 |
+ |
much faster rate than experiment (see figure |
| 2279 |
+ |
\ref{fig:t5peDensities}). The reduced density leads to the inflated |
| 2280 |
+ |
compressibility and expansivity values at higher temperatures seen |
| 2281 |
+ |
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability |
| 2282 |
+ |
and many-body effects are required in order for simulation to overcome |
| 2283 |
+ |
these differences with experiment.\cite{Laasonen93,Donchev06} |
| 2284 |
+ |
|
| 2285 |
+ |
At temperatures below the freezing point for experimental water, the |
| 2286 |
+ |
differences between {\sc sf} and the Ewald summation results are more |
| 2287 |
+ |
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region |
| 2288 |
+ |
indicate a more pronounced transition in the supercooled regim\'{e}, |
| 2289 |
+ |
particularly in the case of {\sc sf} without damping. This points to |
| 2290 |
+ |
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 2291 |
+ |
temperatures below 250K in these simulations. Because the systems are |
| 2292 |
+ |
locked in different regions of phase-space, comparisons between |
| 2293 |
+ |
properties at these temperatures are not exactly fair. This |
| 2294 |
+ |
observation is explored in more detail in section |
| 2295 |
+ |
\ref{sec:t5peDynamics}. |
| 2296 |
+ |
|
| 2297 |
+ |
The final thermodynamic property displayed in figure |
| 2298 |
+ |
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
| 2299 |
+ |
between the Ewald summation and the {\sc sf} technique (and experiment |
| 2300 |
+ |
for that matter). It is known that the dielectric constant is |
| 2301 |
+ |
dependent upon and quite sensitive to the imposed boundary |
| 2302 |
+ |
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the |
| 2303 |
+ |
converged $\epsilon$ values accumulated for the {\sc sf} |
| 2304 |
+ |
simulations. Lack of a damping function results in dielectric |
| 2305 |
+ |
constants significantly smaller than that obtained using the Ewald |
| 2306 |
+ |
sum. Increasing the damping coefficient to 0.2\AA$^{-1}$ improves the |
| 2307 |
+ |
agreement considerably. It should be noted that the choice of the |
| 2308 |
+ |
``Ewald coefficient'' value also has a significant effect on the |
| 2309 |
+ |
calculated value when using the Ewald summation. In the simulations of |
| 2310 |
+ |
TIP5P-E with the Ewald sum, this screening parameter was tethered to |
| 2311 |
+ |
the simulation box size (as was the $R_\textrm{c}$).\cite{Rick04} |
| 2312 |
+ |
Systems with larger screening parameters reported larger dielectric |
| 2313 |
+ |
constant values, the same behavior we see here with {\sc sf}. In |
| 2314 |
+ |
section \ref{sec:dampingDielectric}, this connection is further |
| 2315 |
+ |
explored as optimal damping coefficients are determined for {\sc |
| 2316 |
+ |
sf} for capturing the dielectric behavior. |
| 2317 |
+ |
|
| 2318 |
+ |
\subsection{Dynamic Properties}\label{sec:t5peDynamics} |
| 2319 |
+ |
|
| 2320 |
+ |
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 2321 |
+ |
method, 200ps $NVE$ simulations were performed for each temperature at |
| 2322 |
+ |
the average density reported by the $NPT$ simulations. The |
| 2323 |
+ |
self-diffusion constants ($D$) were calculated with the Einstein |
| 2324 |
+ |
relation using the mean square displacement (MSD), |
| 2325 |
+ |
\begin{equation} |
| 2326 |
+ |
D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t}, |
| 2327 |
+ |
\label{eq:MSD} |
| 2328 |
+ |
\end{equation} |
| 2329 |
+ |
where $t$ is time, and $\mathbf{r}_i$ is the position of particle |
| 2330 |
+ |
$i$.\cite{Allen87} Figure \ref{fig:ExampleMSD} shows an example MSD |
| 2331 |
+ |
plot. As labeled in the figure, MSD plots consist of three distinct |
| 2332 |
+ |
regions: |
| 2333 |
+ |
|
| 2334 |
+ |
\begin{enumerate}[itemsep=0pt] |
| 2335 |
+ |
\item parabolic short-time ballistic motion, |
| 2336 |
+ |
\item linear diffusive regime, and |
| 2337 |
+ |
\item poor statistic region at long-time. |
| 2338 |
+ |
\end{enumerate} |
| 2339 |
+ |
The slope from the linear region (region 2) is used to calculate $D$. |
| 2340 |
+ |
\begin{figure} |
| 2341 |
+ |
\centering |
| 2342 |
+ |
\includegraphics[width=3.5in]{./figures/ExampleMSD.pdf} |
| 2343 |
+ |
\caption{Example plot of mean square displacement verses time. The |
| 2344 |
+ |
left red region is the ballistic motion regime, the middle green |
| 2345 |
+ |
region is the linear diffusive regime, and the right blue region is |
| 2346 |
+ |
the region with poor statistics.} |
| 2347 |
+ |
\label{fig:ExampleMSD} |
| 2348 |
+ |
\end{figure} |
| 2349 |
+ |
|
| 2350 |
+ |
\begin{figure} |
| 2351 |
+ |
\centering |
| 2352 |
+ |
\includegraphics[width=3.5in]{./figures/waterFrame.pdf} |
| 2353 |
+ |
\caption{Body-fixed coordinate frame for a water molecule. The |
| 2354 |
+ |
respective molecular principle axes point in the direction of the |
| 2355 |
+ |
labeled frame axes.} |
| 2356 |
+ |
\label{fig:waterFrame} |
| 2357 |
+ |
\end{figure} |
| 2358 |
+ |
In addition to translational diffusion, reorientational time constants |
| 2359 |
+ |
were calculated for comparisons with the Ewald simulations and with |
| 2360 |
+ |
experiments. These values were determined from 25ps $NVE$ trajectories |
| 2361 |
+ |
through calculation of the orientational time correlation function, |
| 2362 |
+ |
\begin{equation} |
| 2363 |
+ |
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t) |
| 2364 |
+ |
\cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle, |
| 2365 |
+ |
\label{eq:OrientCorr} |
| 2366 |
+ |
\end{equation} |
| 2367 |
+ |
where $P_l$ is the Legendre polynomial of order $l$ and |
| 2368 |
+ |
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along |
| 2369 |
+ |
principle axis $\alpha$. The principle axis frame for these water |
| 2370 |
+ |
molecules is shown in figure \ref{fig:waterFrame}. As an example, |
| 2371 |
+ |
$C_l^y$ is calculated from the time evolution of the unit vector |
| 2372 |
+ |
connecting the two hydrogen atoms. |
| 2373 |
+ |
|
| 2374 |
+ |
\begin{figure} |
| 2375 |
+ |
\centering |
| 2376 |
+ |
\includegraphics[width=3.5in]{./figures/ExampleOrientCorr.pdf} |
| 2377 |
+ |
\caption{Example plots of the orientational autocorrelation functions |
| 2378 |
+ |
for the first and second Legendre polynomials. These curves show the |
| 2379 |
+ |
time decay of the unit vector along the $y$ principle axis.} |
| 2380 |
+ |
\label{fig:OrientCorr} |
| 2381 |
+ |
\end{figure} |
| 2382 |
+ |
From the orientation autocorrelation functions, we can obtain time |
| 2383 |
+ |
constants for rotational relaxation. Figure \ref{fig:OrientCorr} shows |
| 2384 |
+ |
some example plots of orientational autocorrelation functions for the |
| 2385 |
+ |
first and second Legendre polynomials. The relatively short time |
| 2386 |
+ |
portions (between 1 and 3ps for water) of these curves can be fit to |
| 2387 |
+ |
an exponential decay to obtain these constants, and they are directly |
| 2388 |
+ |
comparable to water orientational relaxation times from nuclear |
| 2389 |
+ |
magnetic resonance (NMR). The relaxation constant obtained from |
| 2390 |
+ |
$C_2^y(t)$ is of particular interest because it is about the principle |
| 2391 |
+ |
axis with the minimum moment of inertia and should thereby dominate |
| 2392 |
+ |
the orientational relaxation of the molecule.\cite{Impey82} This means |
| 2393 |
+ |
that $C_2^y(t)$ should provide the best comparison to the NMR |
| 2394 |
+ |
relaxation data. |
| 2395 |
+ |
|
| 2396 |
+ |
\begin{figure} |
| 2397 |
+ |
\centering |
| 2398 |
+ |
\includegraphics[width=5.5in]{./figures/t5peDynamics.pdf} |
| 2399 |
+ |
\caption{Diffusion constants ({\it upper}) and reorientational time |
| 2400 |
+ |
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf} |
| 2401 |
+ |
technique compared with experiment. Data at temperatures less that |
| 2402 |
+ |
0$^\circ$C were not included in the $\tau_2^y$ plot to allow for |
| 2403 |
+ |
easier comparisons in the more relavent temperature regime.} |
| 2404 |
+ |
\label{fig:t5peDynamics} |
| 2405 |
+ |
\end{figure} |
| 2406 |
+ |
Results for the diffusion constants and reorientational time constants |
| 2407 |
+ |
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is |
| 2408 |
+ |
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using |
| 2409 |
+ |
the Ewald sum are reproduced with the {\sc sf} techinque. The enhanced |
| 2410 |
+ |
diffusion at high temperatures are again the product of the lower |
| 2411 |
+ |
densities in comparison with experiment and do not provide any special |
| 2412 |
+ |
insight into differences between the electrostatic summation |
| 2413 |
+ |
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to |
| 2414 |
+ |
diffuse a little faster than with the Ewald sum; however, use of light |
| 2415 |
+ |
to moderate damping results in indistiguishable $D$ values. Though not |
| 2416 |
+ |
apparent in this figure, {\sc sf} values at the lowest temperature are |
| 2417 |
+ |
approximately an order of magnitude lower than with Ewald. These |
| 2418 |
+ |
values support the observation from section \ref{sec:t5peThermo} that |
| 2419 |
+ |
there appeared to be a change to a more glassy-like phase with the |
| 2420 |
+ |
{\sc sf} technique at these lower temperatures. |
| 2421 |
+ |
|
| 2422 |
+ |
The $\tau_2^y$ results in the lower frame of figure |
| 2423 |
+ |
\ref{fig:t5peDynamics} show a much greater difference between the {\sc |
| 2424 |
+ |
sf} results and the Ewald results. At all temperatures shown, TIP5P-E |
| 2425 |
+ |
relaxes faster than experiment with the Ewald sum while tracking |
| 2426 |
+ |
experiment fairly well when using the {\sc sf} technique, independent |
| 2427 |
+ |
of the choice of damping constant. Their are several possible reasons |
| 2428 |
+ |
for this deviation between techniques. The Ewald results were taken |
| 2429 |
+ |
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick |
| 2430 |
+ |
calculation from a 10ps trajectory with {\sc sf} with an $\alpha$ of |
| 2431 |
+ |
0.2\AA$^-1$ at 25$^\circ$C showed a 0.4ps drop in $\tau_2^y$, placing |
| 2432 |
+ |
the result more in line with that obtained using the Ewald sum. These |
| 2433 |
+ |
results support this explanation; however, recomputing the results to |
| 2434 |
+ |
meet a poorer statistical standard is counter-productive. Assuming the |
| 2435 |
+ |
Ewald results are not the product of poor statistics, differences in |
| 2436 |
+ |
techniques to integrate the orientational motion could also play a |
| 2437 |
+ |
role. {\sc shake} is the most commonly used technique for |
| 2438 |
+ |
approximating rigid-body orientational motion,\cite{Ryckaert77} where |
| 2439 |
+ |
as in {\sc oopse}, we maintain and integrate the entire rotation |
| 2440 |
+ |
matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake} |
| 2441 |
+ |
is an iterative constraint technique, if the convergence tolerances |
| 2442 |
+ |
are raised for increased performance, error will accumulate in the |
| 2443 |
+ |
orientational motion. Finally, the Ewald results were calculated using |
| 2444 |
+ |
the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf} |
| 2445 |
+ |
calculations. The additional mode of motion due to the thermostat will |
| 2446 |
+ |
alter the dynamics, resulting in differences between $NVT$ and $NVE$ |
| 2447 |
+ |
results. These differences are increasingly noticable as the |
| 2448 |
+ |
thermostat time constant decreases. |
| 2449 |
+ |
|
| 2450 |
+ |
\section{Damping of Point Multipoles}\label{sec:dampingMultipoles} |
| 2451 |
+ |
|
| 2452 |
+ |
|
| 2453 |
+ |
|
| 2454 |
+ |
\section{Damping and Dielectric Constants}\label{sec:dampingDielectric} |
| 2455 |
+ |
|
| 2456 |
+ |
\section{Conclusions}\label{sec:PairwiseConclusions} |
| 2457 |
+ |
|
| 2458 |
|
The above investigation of pairwise electrostatic summation techniques |
| 2459 |
|
shows that there are viable and computationally efficient alternatives |
| 2460 |
|
to the Ewald summation. These methods are derived from the damped and |
| 2507 |
|
required to obtain the level of accuracy most researchers have come to |
| 2508 |
|
expect. |
| 2509 |
|
|
| 2016 |
– |
\section{An Application: TIP5P-E Water} |
| 2017 |
– |
|
| 2018 |
– |
|
| 2510 |
|
\chapter{\label{chap:water}SIMPLE MODELS FOR WATER} |
| 2511 |
|
|
| 2512 |
|
\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS} |
