| 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} |
| 9 |
|
\usepackage{graphicx} |
| 10 |
|
\usepackage{booktabs} |
| 11 |
|
\usepackage{cite} |
| 12 |
+ |
\usepackage{enumitem} |
| 13 |
|
|
| 14 |
|
\begin{document} |
| 15 |
|
|
| 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, |
| 935 |
|
|
| 936 |
|
\footnotesize |
| 937 |
|
\begin{center} |
| 938 |
< |
\begin{tabular}{@{} ccrrrrrrrr @{}} \\ |
| 938 |
> |
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 939 |
|
\toprule |
| 940 |
|
\toprule |
| 939 |
– |
|
| 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 & |
| 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 @{}} |
| 1030 |
– |
\\ |
| 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 @{}} |
| 1087 |
– |
\\ |
| 1087 |
|
\toprule |
| 1088 |
|
\toprule |
| 1089 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1180 |
|
|
| 1181 |
|
\footnotesize |
| 1182 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1184 |
– |
\\ |
| 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 @{}} |
| 1240 |
– |
\\ |
| 1238 |
|
\toprule |
| 1239 |
|
\toprule |
| 1240 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque |
| 1304 |
|
|
| 1305 |
|
\footnotesize |
| 1306 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1310 |
– |
\\ |
| 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 @{}} |
| 1350 |
– |
\\ |
| 1346 |
|
\toprule |
| 1347 |
|
\toprule |
| 1348 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1394 |
|
|
| 1395 |
|
\footnotesize |
| 1396 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1402 |
– |
\\ |
| 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 @{}} |
| 1443 |
– |
\\ |
| 1437 |
|
\toprule |
| 1438 |
|
\toprule |
| 1439 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1492 |
|
|
| 1493 |
|
\footnotesize |
| 1494 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1502 |
– |
\\ |
| 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 @{}} |
| 1559 |
– |
\\ |
| 1551 |
|
\toprule |
| 1552 |
|
\toprule |
| 1553 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1611 |
|
|
| 1612 |
|
\footnotesize |
| 1613 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1623 |
– |
\\ |
| 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 @{}} |
| 1680 |
– |
\\ |
| 1670 |
|
\toprule |
| 1671 |
|
\toprule |
| 1672 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1731 |
|
|
| 1732 |
|
\footnotesize |
| 1733 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1745 |
– |
\\ |
| 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 @{}} |
| 1802 |
– |
\\ |
| 1790 |
|
\toprule |
| 1791 |
|
\toprule |
| 1792 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 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 |
| 2000 |
> |
molecular pressure ($p(t)$) is given by |
| 2001 |
|
\begin{equation} |
| 2002 |
< |
P(t) = \frac{1}{\textrm{d}V}\sum_\mu |
| 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} |
| 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 transform of |
| 2103 |
< |
$g_\textrm{OO}(r)$. |
| 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, |
| 2130 |
|
that the $g_\textrm{OO}(r)$ is insensitive to the choice of |
| 2131 |
|
electrostatic correction. |
| 2132 |
|
|
| 2133 |
< |
\subsection{Thermodynamic Properties} |
| 2133 |
> |
\subsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 2134 |
|
|
| 2135 |
< |
\subsection{Dynamic Properties} |
| 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 |
| 2507 |
|
required to obtain the level of accuracy most researchers have come to |
| 2508 |
|
expect. |
| 2509 |
|
|
| 2204 |
– |
|
| 2205 |
– |
|
| 2206 |
– |
|
| 2510 |
|
\chapter{\label{chap:water}SIMPLE MODELS FOR WATER} |
| 2511 |
|
|
| 2512 |
|
\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS} |
