| 1 |
< |
\documentclass[12pt]{ndthesis} |
| 1 |
> |
\documentclass[11pt]{ndthesis} |
| 2 |
|
|
| 3 |
|
% some packages for things like equations and graphics |
| 4 |
|
\usepackage{amsmath,bm} |
| 7 |
|
\usepackage{tabularx} |
| 8 |
|
\usepackage{graphicx} |
| 9 |
|
\usepackage{booktabs} |
| 10 |
+ |
\usepackage{cite} |
| 11 |
|
|
| 12 |
|
\begin{document} |
| 13 |
|
|
| 287 |
|
structural and dynamic properties of water compared the same |
| 288 |
|
properties obtained using the Ewald sum. |
| 289 |
|
|
| 290 |
< |
\section{Simple Forms for Pairwise Electrostatics} |
| 290 |
> |
\section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation} |
| 291 |
|
|
| 292 |
|
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 293 |
|
al.} are constructed using two different (and separable) computational |
| 777 |
|
the complementary error function is required). |
| 778 |
|
|
| 779 |
|
The reaction field results illustrates some of that method's |
| 780 |
< |
limitations, primarily that it was developed for use in homogenous |
| 780 |
> |
limitations, primarily that it was developed for use in homogeneous |
| 781 |
|
systems; although it does provide results that are an improvement over |
| 782 |
|
those from an unmodified cutoff. |
| 783 |
|
|
| 902 |
|
all do equivalently well at capturing the direction of both the force |
| 903 |
|
and torque vectors. Using the electrostatic damping improves the |
| 904 |
|
angular behavior significantly for the {\sc sp} and moderately for the |
| 905 |
< |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 905 |
> |
{\sc sf} methods. Over-damping is detrimental to both methods. Again |
| 906 |
|
it is important to recognize that the force vectors cover all |
| 907 |
|
particles in all seven systems, while torque vectors are only |
| 908 |
|
available for neutral molecular groups. Damping is more beneficial to |
| 909 |
|
charged bodies, and this observation is investigated further in |
| 910 |
< |
section \ref{IndividualResults}. |
| 910 |
> |
section \ref{sec:IndividualResults}. |
| 911 |
|
|
| 912 |
|
Although not discussed previously, group based cutoffs can be applied |
| 913 |
|
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 972 |
|
increases, something that is more obvious with group-based cutoffs. |
| 973 |
|
The complimentary error function inserted into the potential weakens |
| 974 |
|
the electrostatic interaction as the value of $\alpha$ is increased. |
| 975 |
< |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 975 |
> |
However, at larger values of $\alpha$, it is possible to over-damp the |
| 976 |
|
electrostatic interaction and to remove it completely. Kast |
| 977 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 978 |
|
values for these types of electrostatic summation methods by fitting |
| 1122 |
|
agreement with {\sc spme} in both energetic and dynamic behavior when |
| 1123 |
|
using the {\sc sf} method with and without damping. The {\sc sp} |
| 1124 |
|
method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly |
| 1125 |
< |
with cutoff radii greater than 12\AA. Overdamping the electrostatics |
| 1125 |
> |
with cutoff radii greater than 12\AA. Over-damping the electrostatics |
| 1126 |
|
reduces the agreement between both these methods and {\sc spme}. |
| 1127 |
|
|
| 1128 |
|
The pure cutoff ({\sc pc}) method performs poorly, again mirroring the |
| 1155 |
|
no damping and only modest improvement for the recommended conditions |
| 1156 |
|
($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA). The |
| 1157 |
|
{\sc sf} method shows modest narrowing across all damping and cutoff |
| 1158 |
< |
ranges of interest. When overdamping these methods, group cutoffs and |
| 1158 |
> |
ranges of interest. When over-damping these methods, group cutoffs and |
| 1159 |
|
the switching function do not improve the force and torque |
| 1160 |
|
directionalities. |
| 1161 |
|
|
| 1834 |
|
|
| 1835 |
|
This system does not appear to show any significant deviations from |
| 1836 |
|
the previously observed results. The {\sc sp} and {\sc sf} methods |
| 1837 |
< |
have aggrements similar to those observed in section |
| 1837 |
> |
have agreements similar to those observed in section |
| 1838 |
|
\ref{sec:WaterResults}. The only significant difference is the |
| 1839 |
|
improvement in the configuration energy differences for the {\sc rf} |
| 1840 |
|
method. This is surprising in that we are introducing an inhomogeneity |
| 1954 |
|
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1955 |
|
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1956 |
|
electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement |
| 1957 |
< |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1957 |
> |
with the power spectrum obtained using the Ewald sum. Over-damping can |
| 1958 |
|
result in underestimates of frequencies of the long-wavelength |
| 1959 |
|
motions.} |
| 1960 |
|
\label{fig:dampInc} |
| 1961 |
|
\end{figure} |
| 1962 |
|
|
| 1963 |
< |
\section{Synopsis of the Pairwise Method Evaluation}\label{sec:PairwiseSynopsis} |
| 1963 |
> |
\section{An Application: TIP5P-E Water}\label{sec:t5peApplied} |
| 1964 |
> |
|
| 1965 |
> |
The above sections focused on the energetics and dynamics of a variety |
| 1966 |
> |
of systems when utilizing the {\sc sp} and {\sc sf} pairwise |
| 1967 |
> |
techniques. A unitary correlation with results obtained using the |
| 1968 |
> |
Ewald summation should result in a successful reproduction of both the |
| 1969 |
> |
static and dynamic properties of any selected system. To test this, |
| 1970 |
> |
we decided to calculate a series of properties for the TIP5P-E water |
| 1971 |
> |
model when using the {\sc sf} technique. |
| 1972 |
> |
|
| 1973 |
> |
The TIP5P-E water model is a variant of Mahoney and Jorgensen's |
| 1974 |
> |
five-point transferable intermolecular potential (TIP5P) model for |
| 1975 |
> |
water.\cite{Mahoney00} TIP5P was developed to reproduce the density |
| 1976 |
> |
maximum anomaly present in liquid water near 4$^\circ$C. As with many |
| 1977 |
> |
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC, |
| 1978 |
> |
and SPC/E), TIP5P was parametrized using a simple cutoff with no |
| 1979 |
> |
long-range electrostatic |
| 1980 |
> |
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
| 1981 |
> |
Without this correction, the pressure term on the central particle |
| 1982 |
> |
from the surroundings is missing. Because they expand to compensate |
| 1983 |
> |
for this added pressure term when this correction is included, systems |
| 1984 |
> |
composed of these particles tend to underpredict the density of water |
| 1985 |
> |
under standard conditions. When using any form of long-range |
| 1986 |
> |
electrostatic correction, it has become common practice to develop or |
| 1987 |
> |
utilize a reparametrized water model that corrects for this |
| 1988 |
> |
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 1989 |
> |
this practice and was optimized specifically for use with the Ewald |
| 1990 |
> |
summation.\cite{Rick04} In his publication, Rick preserved the |
| 1991 |
> |
geometry and point charge magnitudes in TIP5P and focused on altering |
| 1992 |
> |
the Lennard-Jones parameters to correct the density at |
| 1993 |
> |
298K.\cite{Rick04} With the density corrected, he compared common |
| 1994 |
> |
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1995 |
> |
9\AA\ cutoff. |
| 1996 |
> |
|
| 1997 |
> |
In the following sections, we compared these same water properties |
| 1998 |
> |
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1999 |
> |
{\sc sf} technique. In the above evaluation of the pairwise |
| 2000 |
> |
techniques, we observed some flexibility in the choice of parameters. |
| 2001 |
> |
Because of this, the following comparisons include the {\sc sf} |
| 2002 |
> |
technique with a 12\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 2003 |
> |
0.2\AA$^{-1}$, as well as a 9\AA\ cutoff with an $\alpha$ = |
| 2004 |
> |
0.2\AA$^{-1}$. |
| 2005 |
> |
|
| 2006 |
> |
\subsection{Density}\label{sec:t5peDensity} |
| 2007 |
> |
|
| 2008 |
> |
As stated previously, the property that prompted the development of |
| 2009 |
> |
TIP5P-E was the density at 1 atm. The density depends upon the |
| 2010 |
> |
internal pressure of the system in the $NPT$ ensemble, and the |
| 2011 |
> |
calculation of the pressure includes a components from both the |
| 2012 |
> |
kinetic energy and the virial. More specifically, the instantaneous |
| 2013 |
> |
molecular pressure ($P(t)$) is given by |
| 2014 |
> |
\begin{equation} |
| 2015 |
> |
P(t) = \frac{1}{\textrm{d}V}\sum_\mu |
| 2016 |
> |
\left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} |
| 2017 |
> |
+ \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right], |
| 2018 |
> |
\label{eq:MolecularPressure} |
| 2019 |
> |
\end{equation} |
| 2020 |
> |
where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of |
| 2021 |
> |
molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass |
| 2022 |
> |
($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on |
| 2023 |
> |
atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the |
| 2024 |
> |
right term in the brackets of eq. \ref{eq:MolecularPressure}) is |
| 2025 |
> |
directly dependent on the interatomic forces. Since the {\sc sp} |
| 2026 |
> |
method does not modify the forces (see |
| 2027 |
> |
sec. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will |
| 2028 |
> |
be identical to that obtained without an electrostatic correction. |
| 2029 |
> |
The {\sc sf} method does alter the virial component and, by way of the |
| 2030 |
> |
modified pressures, should provide densities more in line with those |
| 2031 |
> |
obtained using the Ewald summation. |
| 2032 |
> |
|
| 2033 |
> |
To compare densities, $NPT$ simulations were performed with the same |
| 2034 |
> |
temperatures as those selected by Rick in his Ewald summation |
| 2035 |
> |
simulations.\cite{Rick04} In order to improve statistics around the |
| 2036 |
> |
density maximum, 3ns trajectories were accumulated at 0, 12.5, and |
| 2037 |
> |
25$^\circ$C, while 2ns trajectories were obtained at all other |
| 2038 |
> |
temperatures. The average densities were calculated from the later |
| 2039 |
> |
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 2040 |
> |
method for accumulating statistics, these sequences were spliced into |
| 2041 |
> |
200 segments to calculate the average density and standard deviation |
| 2042 |
> |
at each temperature.\cite{Mahoney00} |
| 2043 |
> |
|
| 2044 |
> |
\begin{figure} |
| 2045 |
> |
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 2046 |
> |
\caption{Density versus temperature for the TIP5P-E water model when |
| 2047 |
> |
using the Ewald summation (Ref. \cite{Rick04}) and the {\sc sf} method |
| 2048 |
> |
with various parameters. The pressure term from the image-charge shell |
| 2049 |
> |
is larger than that provided by the reciprocal-space portion of the |
| 2050 |
> |
Ewald summation, leading to slightly lower densities. This effect is |
| 2051 |
> |
more visible with the 9\AA\ cutoff, where the image charges exert a |
| 2052 |
> |
greater force on the central particle. The error bars for the {\sc sf} |
| 2053 |
> |
methods show plus or minus the standard deviation of the density |
| 2054 |
> |
measurement at each temperature.} |
| 2055 |
> |
\label{fig:t5peDensities} |
| 2056 |
> |
\end{figure} |
| 2057 |
|
|
| 2058 |
+ |
Figure \ref{fig:t5peDensities} shows the densities calculated for |
| 2059 |
+ |
TIP5P-E using differing electrostatic corrections overlaid on the |
| 2060 |
+ |
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 2061 |
+ |
technique are close to, though typically lower than, those calculated |
| 2062 |
+ |
while using the Ewald summation. These slightly reduced densities |
| 2063 |
+ |
indicate that the pressure component from the image charges at |
| 2064 |
+ |
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 2065 |
+ |
portion of the Ewald summation. Bringing the image charges closer to |
| 2066 |
+ |
the central particle by choosing a 9\AA\ R$_\textrm{c}$ (rather than |
| 2067 |
+ |
the preferred 12\AA\ R$_\textrm{c}$) increases the strength of their |
| 2068 |
+ |
interactions, resulting in a further reduction of the densities. |
| 2069 |
+ |
|
| 2070 |
+ |
Because the strength of the image charge interactions has a noticable |
| 2071 |
+ |
effect on the density, we would expect the use of electrostatic |
| 2072 |
+ |
damping to also play a role in these calculations. Larger values of |
| 2073 |
+ |
$\alpha$ weaken the pair-interactions; and since electrostatic damping |
| 2074 |
+ |
is distance-dependent, force components from the image charges will be |
| 2075 |
+ |
reduced more than those from particles close the the central |
| 2076 |
+ |
charge. This effect is visible in figure \ref{fig:t5peDensities} with |
| 2077 |
+ |
the damped {\sc sf} sums showing slightly higher densities; however, |
| 2078 |
+ |
it is apparent that the choice of cutoff radius plays a much more |
| 2079 |
+ |
important role in the resulting densities. |
| 2080 |
+ |
|
| 2081 |
+ |
As a final note, all of the above density calculations were performed |
| 2082 |
+ |
with systems of 512 water molecules. Rick observed a system sized |
| 2083 |
+ |
dependence of the computed densities when using the Ewald summation, |
| 2084 |
+ |
most likely due to his tying of the convergence parameter to the box |
| 2085 |
+ |
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
| 2086 |
+ |
calculated densities were on average 0.002 to 0.003 g/cm$^3$ lower. A |
| 2087 |
+ |
system size of 256 molecules would force the use of a shorter |
| 2088 |
+ |
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also |
| 2089 |
+ |
lower the densities. Moving to larger systems, as long as the |
| 2090 |
+ |
R$_\textrm{c}$ remains at a fixed value, we would expect the densities |
| 2091 |
+ |
to remain constant. |
| 2092 |
+ |
|
| 2093 |
+ |
\subsection{Liquid Structure}\label{sec:t5peLiqStructure} |
| 2094 |
+ |
|
| 2095 |
+ |
A common function considered when developing and comparing water |
| 2096 |
+ |
models is the oxygen-oxygen radial distribution function |
| 2097 |
+ |
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of |
| 2098 |
+ |
finding a pair of oxygen atoms some distance ($r$) apart relative to a |
| 2099 |
+ |
random distribution at the same density.\cite{Allen87} It is |
| 2100 |
+ |
calculated via |
| 2101 |
+ |
\begin{equation} |
| 2102 |
+ |
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i} |
| 2103 |
+ |
\delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle, |
| 2104 |
+ |
\label{eq:GOOofR} |
| 2105 |
+ |
\end{equation} |
| 2106 |
+ |
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen |
| 2107 |
+ |
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or |
| 2108 |
+ |
neutron scattering experiments through the oxygen-oxygen structure |
| 2109 |
+ |
factor ($S_\textrm{OO}(k)$) by the following relationship: |
| 2110 |
+ |
\begin{equation} |
| 2111 |
+ |
S_\textrm{OO}(k) = 1 + 4\pi\rho |
| 2112 |
+ |
\int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r. |
| 2113 |
+ |
\label{eq:SOOofK} |
| 2114 |
+ |
\end{equation} |
| 2115 |
+ |
Thus, $S_\textrm{OO}(k)$ is related to the Fourier transform of |
| 2116 |
+ |
$g_\textrm{OO}(r)$. |
| 2117 |
+ |
|
| 2118 |
+ |
The expermentally determined $g_\textrm{OO}(r)$ for liquid water has |
| 2119 |
+ |
been compared in great detail with the various common water models, |
| 2120 |
+ |
and TIP5P was found to be in better agreement than other rigid, |
| 2121 |
+ |
non-polarizable models.\cite{Sorenson00} This excellent agreement with |
| 2122 |
+ |
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To |
| 2123 |
+ |
check whether the choice of using the Ewald summation or the {\sc sf} |
| 2124 |
+ |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K |
| 2125 |
+ |
and 1atm were determined for the systems compared in the previous |
| 2126 |
+ |
section. |
| 2127 |
+ |
|
| 2128 |
+ |
\begin{figure} |
| 2129 |
+ |
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 2130 |
+ |
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298K and |
| 2131 |
+ |
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc |
| 2132 |
+ |
sf} technique with varying parameters. Even with the reduced densities |
| 2133 |
+ |
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially |
| 2134 |
+ |
identical.} |
| 2135 |
+ |
\label{fig:t5peGofRs} |
| 2136 |
+ |
\end{figure} |
| 2137 |
+ |
|
| 2138 |
+ |
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc |
| 2139 |
+ |
sf} technique with a various parameters are overlaid on the |
| 2140 |
+ |
$g_\textrm{OO}(r)$ while using the Ewald summation. The differences in |
| 2141 |
+ |
density do not appear to have any effect on the liquid structure as |
| 2142 |
+ |
the $g_\textrm{OO}(r)$s are indistinquishable. These results indicate |
| 2143 |
+ |
that the $g_\textrm{OO}(r)$ is insensitive to the choice of |
| 2144 |
+ |
electrostatic correction. |
| 2145 |
+ |
|
| 2146 |
+ |
\subsection{Thermodynamic Properties} |
| 2147 |
+ |
|
| 2148 |
+ |
\subsection{Dynamic Properties} |
| 2149 |
+ |
|
| 2150 |
+ |
\section{Conclusions}\label{sec:PairwiseConclusions} |
| 2151 |
+ |
|
| 2152 |
|
The above investigation of pairwise electrostatic summation techniques |
| 2153 |
|
shows that there are viable and computationally efficient alternatives |
| 2154 |
|
to the Ewald summation. These methods are derived from the damped and |
| 2201 |
|
required to obtain the level of accuracy most researchers have come to |
| 2202 |
|
expect. |
| 2203 |
|
|
| 2016 |
– |
\section{An Application: TIP5P-E Water} |
| 2204 |
|
|
| 2205 |
|
|
| 2206 |
+ |
|
| 2207 |
|
\chapter{\label{chap:water}SIMPLE MODELS FOR WATER} |
| 2208 |
|
|
| 2209 |
|
\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS} |
