| 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 |
|
|
| 43 |
|
|
| 44 |
|
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} |
| 45 |
|
|
| 46 |
< |
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION |
| 44 |
< |
TECHNIQUES} |
| 45 |
< |
|
| 46 |
< |
In molecular simulations, proper accumulation of the electrostatic |
| 47 |
< |
interactions is essential and is one of the most |
| 48 |
< |
computationally-demanding tasks. The common molecular mechanics force |
| 49 |
< |
fields represent atomic sites with full or partial charges protected |
| 50 |
< |
by Lennard-Jones (short range) interactions. This means that nearly |
| 51 |
< |
every pair interaction involves a calculation of charge-charge forces. |
| 52 |
< |
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
| 53 |
< |
interactions quickly become the most expensive part of molecular |
| 54 |
< |
simulations. Historically, the electrostatic pair interaction would |
| 55 |
< |
not have decayed appreciably within the typical box lengths that could |
| 56 |
< |
be feasibly simulated. In the larger systems that are more typical of |
| 57 |
< |
modern simulations, large cutoffs should be used to incorporate |
| 58 |
< |
electrostatics correctly. |
| 59 |
< |
|
| 60 |
< |
There have been many efforts to address the proper and practical |
| 61 |
< |
handling of electrostatic interactions, and these have resulted in a |
| 62 |
< |
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are |
| 63 |
< |
typically classified as implicit methods (i.e., continuum dielectrics, |
| 64 |
< |
static dipolar fields),\cite{Born20,Grossfield00} explicit methods |
| 65 |
< |
(i.e., Ewald summations, interaction shifting or |
| 66 |
< |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 67 |
< |
reaction field type methods, fast multipole |
| 68 |
< |
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 69 |
< |
often preferred because they physically incorporate solvent molecules |
| 70 |
< |
in the system of interest, but these methods are sometimes difficult |
| 71 |
< |
to utilize because of their high computational cost.\cite{Roux99} In |
| 72 |
< |
addition to the computational cost, there have been some questions |
| 73 |
< |
regarding possible artifacts caused by the inherent periodicity of the |
| 74 |
< |
explicit Ewald summation.\cite{Tobias01} |
| 75 |
< |
|
| 76 |
< |
In this chapter, we focus on a new set of pairwise methods devised by |
| 77 |
< |
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 78 |
< |
methods along with a few other mixed methods (i.e. reaction field) are |
| 79 |
< |
compared with the smooth particle mesh Ewald |
| 80 |
< |
sum,\cite{Onsager36,Essmann99} which is our reference method for |
| 81 |
< |
handling long-range electrostatic interactions. The new methods for |
| 82 |
< |
handling electrostatics have the potential to scale linearly with |
| 83 |
< |
increasing system size since they involve only a simple modification |
| 84 |
< |
to the direct pairwise sum. They also lack the added periodicity of |
| 85 |
< |
the Ewald sum, so they can be used for systems which are non-periodic |
| 86 |
< |
or which have one- or two-dimensional periodicity. Below, these |
| 87 |
< |
methods are evaluated using a variety of model systems to |
| 88 |
< |
establish their usability in molecular simulations. |
| 89 |
< |
|
| 90 |
< |
\section{The Ewald Sum} |
| 91 |
< |
|
| 92 |
< |
The complete accumulation of the electrostatic interactions in a system with |
| 93 |
< |
periodic boundary conditions (PBC) requires the consideration of the |
| 94 |
< |
effect of all charges within a (cubic) simulation box as well as those |
| 95 |
< |
in the periodic replicas, |
| 96 |
< |
\begin{equation} |
| 97 |
< |
V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime |
| 98 |
< |
\left[ \sum_{i=1}^N\sum_{j=1}^N \phi |
| 99 |
< |
\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) |
| 100 |
< |
\right], |
| 101 |
< |
\label{eq:PBCSum} |
| 102 |
< |
\end{equation} |
| 103 |
< |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
| 104 |
< |
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
| 105 |
< |
prime indicates $i = j$ are neglected for $\mathbf{n} = |
| 106 |
< |
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
| 107 |
< |
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
| 108 |
< |
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
| 109 |
< |
$j$, and $\phi$ is the solution to Poisson's equation |
| 110 |
< |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for |
| 111 |
< |
charge-charge interactions). In the case of monopole electrostatics, |
| 112 |
< |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for |
| 113 |
< |
non-neutral systems. |
| 114 |
< |
|
| 115 |
< |
The electrostatic summation problem was originally studied by Ewald |
| 116 |
< |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 117 |
< |
took was to convert this conditionally convergent sum into two |
| 118 |
< |
absolutely convergent summations: a short-ranged real-space summation |
| 119 |
< |
and a long-ranged reciprocal-space summation, |
| 120 |
< |
\begin{equation} |
| 121 |
< |
\begin{split} |
| 122 |
< |
V_\textrm{elec} = \frac{1}{2}& |
| 123 |
< |
\sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime |
| 124 |
< |
\frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)} |
| 125 |
< |
{|\mathbf{r}_{ij}+\mathbf{n}|} \\ |
| 126 |
< |
&+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} |
| 127 |
< |
\exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) |
| 128 |
< |
\cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ |
| 129 |
< |
&- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 |
| 130 |
< |
+ \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3} |
| 131 |
< |
\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
| 132 |
< |
\end{split} |
| 133 |
< |
\label{eq:EwaldSum} |
| 134 |
< |
\end{equation} |
| 135 |
< |
where $\alpha$ is the damping or convergence parameter with units of |
| 136 |
< |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 137 |
< |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 138 |
< |
constant of the surrounding medium. The final two terms of |
| 139 |
< |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 140 |
< |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 141 |
< |
dipolar term was neglected in early applications in molecular |
| 142 |
< |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 143 |
< |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 144 |
< |
dipole moment which is magnified through replication of the periodic |
| 145 |
< |
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
| 146 |
< |
system is said to be using conducting (or ``tin-foil'') boundary |
| 147 |
< |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 148 |
< |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 149 |
< |
time. Initially, due to the small system sizes that could be |
| 150 |
< |
simulated feasibly, the entire simulation box was replicated to |
| 151 |
< |
convergence. In more modern simulations, the systems have grown large |
| 152 |
< |
enough that a real-space cutoff could potentially give convergent |
| 153 |
< |
behavior. Indeed, it has been observed that with the choice of a |
| 154 |
< |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 155 |
< |
rapidly convergent and small relative to the real-space |
| 156 |
< |
portion.\cite{Karasawa89,Kolafa92} |
| 157 |
< |
|
| 158 |
< |
\begin{figure} |
| 159 |
< |
\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf} |
| 160 |
< |
\caption{The change in the need for the Ewald sum with |
| 161 |
< |
increasing computational power. A:~Initially, only small systems |
| 162 |
< |
could be studied, and the Ewald sum replicated the simulation box to |
| 163 |
< |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 164 |
< |
convergence for the larger systems of charges that are common today.} |
| 165 |
< |
\label{fig:ewaldTime} |
| 166 |
< |
\end{figure} |
| 167 |
< |
|
| 168 |
< |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 169 |
< |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 170 |
< |
the computational cost between the direct and reciprocal-space |
| 171 |
< |
portions of the summation. The choice of this value allows one to |
| 172 |
< |
select whether the real-space or reciprocal space portion of the |
| 173 |
< |
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
| 174 |
< |
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 175 |
< |
$\alpha$ and thoughtful algorithm development, this cost can be |
| 176 |
< |
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 177 |
< |
taken to reduce the cost of the Ewald summation even further is to set |
| 178 |
< |
$\alpha$ such that the real-space interactions decay rapidly, allowing |
| 179 |
< |
for a short spherical cutoff. Then the reciprocal space summation is |
| 180 |
< |
optimized. These optimizations usually involve utilization of the |
| 181 |
< |
fast Fourier transform (FFT),\cite{Hockney81} leading to the |
| 182 |
< |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 183 |
< |
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 184 |
< |
methods, the cost of the reciprocal-space portion of the Ewald |
| 185 |
< |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
| 186 |
< |
\log N)$. |
| 187 |
< |
|
| 188 |
< |
These developments and optimizations have made the use of the Ewald |
| 189 |
< |
summation routine in simulations with periodic boundary |
| 190 |
< |
conditions. However, in certain systems, such as vapor-liquid |
| 191 |
< |
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 192 |
< |
can prove problematic. The Ewald sum has been reformulated to handle |
| 193 |
< |
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these |
| 194 |
< |
methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 195 |
< |
recently, there have been several successful efforts toward reducing |
| 196 |
< |
the computational cost of 2-D lattice |
| 197 |
< |
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 198 |
< |
bringing them more in line with the cost of the full 3-D summation. |
| 199 |
< |
|
| 200 |
< |
Several studies have recognized that the inherent periodicity in the |
| 201 |
< |
Ewald sum can also have an effect on three-dimensional |
| 202 |
< |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 203 |
< |
Solvated proteins are essentially kept at high concentration due to |
| 204 |
< |
the periodicity of the electrostatic summation method. In these |
| 205 |
< |
systems, the more compact folded states of a protein can be |
| 206 |
< |
artificially stabilized by the periodic replicas introduced by the |
| 207 |
< |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 208 |
< |
considering the use of the Ewald summation where the assumed |
| 209 |
< |
periodicity would introduce spurious effects in the system dynamics. |
| 210 |
< |
|
| 211 |
< |
|
| 212 |
< |
\section{The Wolf and Zahn Methods} |
| 213 |
< |
|
| 214 |
< |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 215 |
< |
for the accurate accumulation of electrostatic interactions in an |
| 216 |
< |
efficient pairwise fashion. This procedure lacks the inherent |
| 217 |
< |
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
| 218 |
< |
observed that the electrostatic interaction is effectively |
| 219 |
< |
short-ranged in condensed phase systems and that neutralization of the |
| 220 |
< |
charge contained within the cutoff radius is crucial for potential |
| 221 |
< |
stability. They devised a pairwise summation method that ensures |
| 222 |
< |
charge neutrality and gives results similar to those obtained with the |
| 223 |
< |
Ewald summation. The resulting shifted Coulomb potential includes |
| 224 |
< |
image-charges subtracted out through placement on the cutoff sphere |
| 225 |
< |
and a distance-dependent damping function (identical to that seen in |
| 226 |
< |
the real-space portion of the Ewald sum) to aid convergence |
| 227 |
< |
\begin{equation} |
| 228 |
< |
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}} |
| 229 |
< |
- \lim_{r_{ij}\rightarrow R_\textrm{c}} |
| 230 |
< |
\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 231 |
< |
\label{eq:WolfPot} |
| 232 |
< |
\end{equation} |
| 233 |
< |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 234 |
< |
potential. However, neutralizing the charge contained within each |
| 235 |
< |
cutoff sphere requires the placement of a self-image charge on the |
| 236 |
< |
surface of the cutoff sphere. This additional self-term in the total |
| 237 |
< |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
| 238 |
< |
Madelung energies for many crystals. |
| 239 |
< |
|
| 240 |
< |
In order to use their charge-neutralized potential in molecular |
| 241 |
< |
dynamics simulations, Wolf \textit{et al.} suggested taking the |
| 242 |
< |
derivative of this potential prior to evaluation of the limit. This |
| 243 |
< |
procedure gives an expression for the forces, |
| 244 |
< |
\begin{equation} |
| 245 |
< |
\begin{split} |
| 246 |
< |
F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{& |
| 247 |
< |
\Biggr[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}} |
| 248 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}} |
| 249 |
< |
\Biggr]\\ |
| 250 |
< |
&-\Biggr[ |
| 251 |
< |
\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2} |
| 252 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 253 |
< |
\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} |
| 254 |
< |
\Biggr]\Biggr\}, |
| 255 |
< |
\end{split} |
| 256 |
< |
\label{eq:WolfForces} |
| 257 |
< |
\end{equation} |
| 258 |
< |
that incorporates both image charges and damping of the electrostatic |
| 259 |
< |
interaction. |
| 260 |
< |
|
| 261 |
< |
More recently, Zahn \textit{et al.} investigated these potential and |
| 262 |
< |
force expressions for use in simulations involving water.\cite{Zahn02} |
| 263 |
< |
In their work, they pointed out that the forces and derivative of |
| 264 |
< |
the potential are not commensurate. Attempts to use both |
| 265 |
< |
eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 266 |
< |
to poor energy conservation. They correctly observed that taking the |
| 267 |
< |
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 268 |
< |
derivatives gives forces for a different potential energy function |
| 269 |
< |
than the one shown in eq. (\ref{eq:WolfPot}). |
| 270 |
< |
|
| 271 |
< |
Zahn \textit{et al.} introduced a modified form of this summation |
| 272 |
< |
method as a way to use the technique in Molecular Dynamics |
| 273 |
< |
simulations. They proposed a new damped Coulomb potential, |
| 274 |
< |
\begin{equation} |
| 275 |
< |
\begin{split} |
| 276 |
< |
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{& |
| 277 |
< |
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} \\ |
| 278 |
< |
&- \left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
| 279 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 280 |
< |
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
| 281 |
< |
\right]\left(r_{ij}-R_\mathrm{c}\right)\Biggr\}, |
| 282 |
< |
\end{split} |
| 283 |
< |
\label{eq:ZahnPot} |
| 284 |
< |
\end{equation} |
| 285 |
< |
and showed that this potential does fairly well at capturing the |
| 286 |
< |
structural and dynamic properties of water compared the same |
| 287 |
< |
properties obtained using the Ewald sum. |
| 288 |
< |
|
| 289 |
< |
\section{Simple Forms for Pairwise Electrostatics} |
| 290 |
< |
|
| 291 |
< |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 292 |
< |
al.} are constructed using two different (and separable) computational |
| 293 |
< |
tricks: |
| 294 |
< |
|
| 295 |
< |
\begin{enumerate} |
| 296 |
< |
\item shifting through the use of image charges, and |
| 297 |
< |
\item damping the electrostatic interaction. |
| 298 |
< |
\end{enumerate} |
| 299 |
< |
Wolf \textit{et al.} treated the |
| 300 |
< |
development of their summation method as a progressive application of |
| 301 |
< |
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 302 |
< |
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 303 |
< |
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
| 304 |
< |
both techniques. It is possible, however, to separate these |
| 305 |
< |
tricks and study their effects independently. |
| 306 |
< |
|
| 307 |
< |
Starting with the original observation that the effective range of the |
| 308 |
< |
electrostatic interaction in condensed phases is considerably less |
| 309 |
< |
than $r^{-1}$, either the cutoff sphere neutralization or the |
| 310 |
< |
distance-dependent damping technique could be used as a foundation for |
| 311 |
< |
a new pairwise summation method. Wolf \textit{et al.} made the |
| 312 |
< |
observation that charge neutralization within the cutoff sphere plays |
| 313 |
< |
a significant role in energy convergence; therefore we will begin our |
| 314 |
< |
analysis with the various shifted forms that maintain this charge |
| 315 |
< |
neutralization. We can evaluate the methods of Wolf |
| 316 |
< |
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 317 |
< |
shifted potential, |
| 318 |
< |
\begin{equation} |
| 319 |
< |
V_\textrm{SP}(r) = \begin{cases} |
| 320 |
< |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 321 |
< |
R_\textrm{c} |
| 322 |
< |
\end{cases}, |
| 323 |
< |
\label{eq:shiftingPotForm} |
| 324 |
< |
\end{equation} |
| 325 |
< |
and shifted force, |
| 326 |
< |
\begin{equation} |
| 327 |
< |
V_\textrm{SF}(r) = \begin{cases} |
| 328 |
< |
v(r) - v_\textrm{c} |
| 329 |
< |
- \left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 330 |
< |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 331 |
< |
\end{cases}, |
| 332 |
< |
\label{eq:shiftingForm} |
| 333 |
< |
\end{equation} |
| 334 |
< |
functions where $v(r)$ is the unshifted form of the potential, and |
| 335 |
< |
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
| 336 |
< |
that both the potential and the forces goes to zero at the cutoff |
| 337 |
< |
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
| 338 |
< |
potential is smooth at the cutoff radius |
| 339 |
< |
($R_\textrm{c}$).\cite{Allen87} |
| 340 |
< |
|
| 341 |
< |
The forces associated with the shifted potential are simply the forces |
| 342 |
< |
of the unshifted potential itself (when inside the cutoff sphere), |
| 343 |
< |
\begin{equation} |
| 344 |
< |
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 345 |
< |
\end{equation} |
| 346 |
< |
and are zero outside. Inside the cutoff sphere, the forces associated |
| 347 |
< |
with the shifted force form can be written, |
| 348 |
< |
\begin{equation} |
| 349 |
< |
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 350 |
< |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 351 |
< |
\end{equation} |
| 352 |
< |
|
| 353 |
< |
If the potential, $v(r)$, is taken to be the normal Coulomb potential, |
| 354 |
< |
\begin{equation} |
| 355 |
< |
v(r) = \frac{q_i q_j}{r}, |
| 356 |
< |
\label{eq:Coulomb} |
| 357 |
< |
\end{equation} |
| 358 |
< |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 359 |
< |
al.}'s undamped prescription: |
| 360 |
< |
\begin{equation} |
| 361 |
< |
V_\textrm{SP}(r) = |
| 362 |
< |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 363 |
< |
r\leqslant R_\textrm{c}, |
| 364 |
< |
\label{eq:SPPot} |
| 365 |
< |
\end{equation} |
| 366 |
< |
with associated forces, |
| 367 |
< |
\begin{equation} |
| 368 |
< |
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) |
| 369 |
< |
\quad r\leqslant R_\textrm{c}. |
| 370 |
< |
\label{eq:SPForces} |
| 371 |
< |
\end{equation} |
| 372 |
< |
These forces are identical to the forces of the standard Coulomb |
| 373 |
< |
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 374 |
< |
\textit{et al.} as undesirable. They pointed out that the effect of |
| 375 |
< |
the image charges is neglected in the forces when this form is |
| 376 |
< |
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
| 377 |
< |
molecular dynamics. Additionally, there is a discontinuity in the |
| 378 |
< |
forces at the cutoff radius which results in energy drift during MD |
| 379 |
< |
simulations. |
| 380 |
< |
|
| 381 |
< |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 382 |
< |
will give, |
| 383 |
< |
\begin{equation} |
| 384 |
< |
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}} |
| 385 |
< |
+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] |
| 386 |
< |
\quad r\leqslant R_\textrm{c}. |
| 387 |
< |
\label{eq:SFPot} |
| 388 |
< |
\end{equation} |
| 389 |
< |
with associated forces, |
| 390 |
< |
\begin{equation} |
| 391 |
< |
F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) |
| 392 |
< |
\quad r\leqslant R_\textrm{c}. |
| 393 |
< |
\label{eq:SFForces} |
| 394 |
< |
\end{equation} |
| 395 |
< |
This formulation has the benefits that there are no discontinuities at |
| 396 |
< |
the cutoff radius, while the neutralizing image charges are present in |
| 397 |
< |
both the energy and force expressions. It would be simple to add the |
| 398 |
< |
self-neutralizing term back when computing the total energy of the |
| 399 |
< |
system, thereby maintaining the agreement with the Madelung energies. |
| 400 |
< |
A side effect of this treatment is the alteration in the shape of the |
| 401 |
< |
potential that comes from the derivative term. Thus, a degree of |
| 402 |
< |
clarity about agreement with the empirical potential is lost in order |
| 403 |
< |
to gain functionality in dynamics simulations. |
| 404 |
< |
|
| 405 |
< |
Wolf \textit{et al.} originally discussed the energetics of the |
| 406 |
< |
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
| 407 |
< |
insufficient for accurate determination of the energy with reasonable |
| 408 |
< |
cutoff distances. The calculated Madelung energies fluctuated around |
| 409 |
< |
the expected value as the cutoff radius was increased, but the |
| 410 |
< |
oscillations converged toward the correct value.\cite{Wolf99} A |
| 411 |
< |
damping function was incorporated to accelerate the convergence; and |
| 412 |
< |
though alternative forms for the damping function could be |
| 413 |
< |
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 414 |
< |
chosen to mirror the effective screening used in the Ewald summation. |
| 415 |
< |
Incorporating this error function damping into the simple Coulomb |
| 416 |
< |
potential, |
| 417 |
< |
\begin{equation} |
| 418 |
< |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 419 |
< |
\label{eq:dampCoulomb} |
| 420 |
< |
\end{equation} |
| 421 |
< |
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 422 |
< |
\begin{equation} |
| 423 |
< |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} |
| 424 |
< |
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) |
| 425 |
< |
\quad r\leqslant R_\textrm{c}, |
| 426 |
< |
\label{eq:DSPPot} |
| 427 |
< |
\end{equation} |
| 428 |
< |
with associated forces, |
| 429 |
< |
\begin{equation} |
| 430 |
< |
F_{\textrm{DSP}}(r) = q_iq_j |
| 431 |
< |
\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} |
| 432 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) |
| 433 |
< |
\quad r\leqslant R_\textrm{c}. |
| 434 |
< |
\label{eq:DSPForces} |
| 435 |
< |
\end{equation} |
| 436 |
< |
Again, this damped shifted potential suffers from a |
| 437 |
< |
force-discontinuity at the cutoff radius, and the image charges play |
| 438 |
< |
no role in the forces. To remedy these concerns, one may derive a |
| 439 |
< |
{\sc sf} variant by including the derivative term in |
| 440 |
< |
eq. (\ref{eq:shiftingForm}), |
| 441 |
< |
\begin{equation} |
| 442 |
< |
\begin{split} |
| 443 |
< |
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 444 |
< |
\frac{\mathrm{erfc}\left(\alpha r\right)}{r} |
| 445 |
< |
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ |
| 446 |
< |
&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
| 447 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 448 |
< |
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
| 449 |
< |
\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} |
| 450 |
< |
\quad r\leqslant R_\textrm{c}. |
| 451 |
< |
\label{eq:DSFPot} |
| 452 |
< |
\end{split} |
| 453 |
< |
\end{equation} |
| 454 |
< |
The derivative of the above potential will lead to the following forces, |
| 455 |
< |
\begin{equation} |
| 456 |
< |
\begin{split} |
| 457 |
< |
F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 458 |
< |
\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} |
| 459 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ |
| 460 |
< |
&- \left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)} |
| 461 |
< |
{R_{\textrm{c}}^2} |
| 462 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 463 |
< |
\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} |
| 464 |
< |
\right)\Biggr{]} |
| 465 |
< |
\quad r\leqslant R_\textrm{c}. |
| 466 |
< |
\label{eq:DSFForces} |
| 467 |
< |
\end{split} |
| 468 |
< |
\end{equation} |
| 469 |
< |
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
| 470 |
< |
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
| 471 |
< |
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
| 472 |
< |
|
| 473 |
< |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 474 |
< |
derived by Zahn \textit{et al.}; however, there are two important |
| 475 |
< |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 476 |
< |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 477 |
< |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 478 |
< |
in the Zahn potential, resulting in a potential discontinuity as |
| 479 |
< |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 480 |
< |
portion is different. The missing $v_\textrm{c}$ term would not |
| 481 |
< |
affect molecular dynamics simulations (although the computed energy |
| 482 |
< |
would be expected to have sudden jumps as particle distances crossed |
| 483 |
< |
$R_c$). The sign problem is a potential source of errors, however. |
| 484 |
< |
In fact, it introduces a discontinuity in the forces at the cutoff, |
| 485 |
< |
because the force function is shifted in the wrong direction and |
| 486 |
< |
doesn't cross zero at $R_\textrm{c}$. |
| 487 |
< |
|
| 488 |
< |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 489 |
< |
electrostatic summation method in which the potential and forces are |
| 490 |
< |
continuous at the cutoff radius and which incorporates the damping |
| 491 |
< |
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
| 492 |
< |
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
| 493 |
< |
{\sc sf}, damping) are at reproducing the correct electrostatic |
| 494 |
< |
summation performed by the Ewald sum. |
| 495 |
< |
|
| 496 |
< |
|
| 497 |
< |
\section{Evaluating Pairwise Summation Techniques} |
| 498 |
< |
|
| 499 |
< |
In classical molecular mechanics simulations, there are two primary |
| 500 |
< |
techniques utilized to obtain information about the system of |
| 501 |
< |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 502 |
< |
techniques utilize pairwise summations of interactions between |
| 503 |
< |
particle sites, but they use these summations in different ways. |
| 504 |
< |
|
| 505 |
< |
In MC, the potential energy difference between configurations dictates |
| 506 |
< |
the progression of MC sampling. Going back to the origins of this |
| 507 |
< |
method, the acceptance criterion for the canonical ensemble laid out |
| 508 |
< |
by Metropolis \textit{et al.} states that a subsequent configuration |
| 509 |
< |
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
| 510 |
< |
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
| 511 |
< |
Maintaining the correct $\Delta E$ when using an alternate method for |
| 512 |
< |
handling the long-range electrostatics will ensure proper sampling |
| 513 |
< |
from the ensemble. |
| 514 |
< |
|
| 515 |
< |
In MD, the derivative of the potential governs how the system will |
| 516 |
< |
progress in time. Consequently, the force and torque vectors on each |
| 517 |
< |
body in the system dictate how the system evolves. If the magnitude |
| 518 |
< |
and direction of these vectors are similar when using alternate |
| 519 |
< |
electrostatic summation techniques, the dynamics in the short term |
| 520 |
< |
will be indistinguishable. Because error in MD calculations is |
| 521 |
< |
cumulative, one should expect greater deviation at longer times, |
| 522 |
< |
although methods which have large differences in the force and torque |
| 523 |
< |
vectors will diverge from each other more rapidly. |
| 524 |
< |
|
| 525 |
< |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 526 |
< |
|
| 527 |
< |
The pairwise summation techniques (outlined in section |
| 528 |
< |
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 529 |
< |
studying the energy differences between conformations. We took the |
| 530 |
< |
{\sc spme}-computed energy difference between two conformations to be the |
| 531 |
< |
correct behavior. An ideal performance by an alternative method would |
| 532 |
< |
reproduce these energy differences exactly (even if the absolute |
| 533 |
< |
energies calculated by the methods are different). Since none of the |
| 534 |
< |
methods provide exact energy differences, we used linear least squares |
| 535 |
< |
regressions of energy gap data to evaluate how closely the methods |
| 536 |
< |
mimicked the Ewald energy gaps. Unitary results for both the |
| 537 |
< |
correlation (slope) and correlation coefficient for these regressions |
| 538 |
< |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 539 |
< |
Sample correlation plots for two alternate methods are shown in |
| 540 |
< |
Fig. \ref{fig:linearFit}. |
| 541 |
< |
|
| 542 |
< |
\begin{figure} |
| 543 |
< |
\centering |
| 544 |
< |
\includegraphics[width = \linewidth]{./figures/dualLinear.pdf} |
| 545 |
< |
\caption{Example least squares regressions of the configuration energy |
| 546 |
< |
differences for SPC/E water systems. The upper plot shows a data set |
| 547 |
< |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 548 |
< |
shows a data set with a good correlation coefficient.} |
| 549 |
< |
\label{fig:linearFit} |
| 550 |
< |
\end{figure} |
| 551 |
< |
|
| 552 |
< |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 553 |
< |
were represented using 500 independent configurations. Thus, each of |
| 554 |
< |
the alternative (non-Ewald) electrostatic summation methods was |
| 555 |
< |
evaluated using an accumulated 873,250 configurational energy |
| 556 |
< |
differences. |
| 557 |
< |
|
| 558 |
< |
Results and discussion for the individual analysis of each of the |
| 559 |
< |
system types appear in sections \ref{sec:SystemResults}, while the |
| 560 |
< |
cumulative results over all the investigated systems appear below in |
| 561 |
< |
sections \ref{sec:EnergyResults}. |
| 562 |
< |
|
| 563 |
< |
\subsection{Molecular Dynamics and the Force and Torque |
| 564 |
< |
Vectors}\label{sec:MDMethods} We evaluated the pairwise methods |
| 565 |
< |
(outlined in section \ref{sec:ESMethods}) for use in MD simulations by |
| 566 |
< |
comparing the force and torque vectors with those obtained using the |
| 567 |
< |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 568 |
< |
direction of these vectors on each of the bodies in the system were |
| 569 |
< |
analyzed. For the magnitude of these vectors, linear least squares |
| 570 |
< |
regression analyses were performed as described previously for |
| 571 |
< |
comparing $\Delta E$ values. Instead of a single energy difference |
| 572 |
< |
between two system configurations, we compared the magnitudes of the |
| 573 |
< |
forces (and torques) on each molecule in each configuration. For a |
| 574 |
< |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 575 |
< |
vectors and 1000 torque vectors. With 500 configurations, this |
| 576 |
< |
results in 520,000 force and 500,000 torque vector comparisons. |
| 577 |
< |
Additionally, data from seven different system types was aggregated |
| 578 |
< |
before the comparison was made. |
| 579 |
< |
|
| 580 |
< |
The {\it directionality} of the force and torque vectors was |
| 581 |
< |
investigated through measurement of the angle ($\theta$) formed |
| 582 |
< |
between those computed from the particular method and those from {\sc spme}, |
| 583 |
< |
\begin{equation} |
| 584 |
< |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} |
| 585 |
< |
\cdot \hat{F}_\textrm{M}\right), |
| 586 |
< |
\end{equation} |
| 587 |
< |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 588 |
< |
vector computed using method M. Each of these $\theta$ values was |
| 589 |
< |
accumulated in a distribution function and weighted by the area on the |
| 590 |
< |
unit sphere. Since this distribution is a measure of angular error |
| 591 |
< |
between two different electrostatic summation methods, there is no |
| 592 |
< |
{\it a priori} reason for the profile to adhere to any specific |
| 593 |
< |
shape. Thus, gaussian fits were used to measure the width of the |
| 594 |
< |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 595 |
< |
each of these fits and was used to compare distribution widths. |
| 596 |
< |
Values of $\sigma^2$ near zero indicate vector directions |
| 597 |
< |
indistinguishable from those calculated when using the reference |
| 598 |
< |
method ({\sc spme}). |
| 599 |
< |
|
| 600 |
< |
\subsection{Short-time Dynamics} |
| 601 |
< |
|
| 602 |
< |
The effects of the alternative electrostatic summation methods on the |
| 603 |
< |
short-time dynamics of charged systems were evaluated by considering a |
| 604 |
< |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 605 |
< |
performing pairwise methods was used in this comparison. The NaCl |
| 606 |
< |
crystal was chosen to avoid possible complications from the treatment |
| 607 |
< |
of orientational motion in molecular systems. All systems were |
| 608 |
< |
started with the same initial positions and velocities. Simulations |
| 609 |
< |
were performed under the microcanonical ensemble, and velocity |
| 610 |
< |
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 611 |
< |
of the trajectories, |
| 612 |
< |
\begin{equation} |
| 613 |
< |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 614 |
< |
\label{eq:vCorr} |
| 615 |
< |
\end{equation} |
| 616 |
< |
Velocity autocorrelation functions require detailed short time data, |
| 617 |
< |
thus velocity information was saved every 2 fs over 10 ps |
| 618 |
< |
trajectories. Because the NaCl crystal is composed of two different |
| 619 |
< |
atom types, the average of the two resulting velocity autocorrelation |
| 620 |
< |
functions was used for comparisons. |
| 621 |
< |
|
| 622 |
< |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 623 |
< |
|
| 624 |
< |
The effects of the same subset of alternative electrostatic methods on |
| 625 |
< |
the {\it long-time} dynamics of charged systems were evaluated using |
| 626 |
< |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 627 |
< |
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 628 |
< |
autocorrelation function, \begin{equation} I(\omega) = |
| 629 |
< |
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 630 |
< |
\label{eq:powerSpec} |
| 631 |
< |
\end{equation} |
| 632 |
< |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 633 |
< |
NaCl crystal is composed of two different atom types, the average of |
| 634 |
< |
the two resulting power spectra was used for comparisons. Simulations |
| 635 |
< |
were performed under the microcanonical ensemble, and velocity |
| 636 |
< |
information was saved every 5~fs over 100~ps trajectories. |
| 637 |
< |
|
| 638 |
< |
\subsection{Representative Simulations}\label{sec:RepSims} |
| 639 |
< |
A variety of representative molecular simulations were analyzed to |
| 640 |
< |
determine the relative effectiveness of the pairwise summation |
| 641 |
< |
techniques in reproducing the energetics and dynamics exhibited by |
| 642 |
< |
{\sc spme}. We wanted to span the space of typical molecular |
| 643 |
< |
simulations (i.e. from liquids of neutral molecules to ionic |
| 644 |
< |
crystals), so the systems studied were: |
| 645 |
< |
|
| 646 |
< |
\begin{enumerate} |
| 647 |
< |
\item liquid water (SPC/E),\cite{Berendsen87} |
| 648 |
< |
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 649 |
< |
\item NaCl crystals, |
| 650 |
< |
\item NaCl melts, |
| 651 |
< |
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 652 |
< |
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 653 |
< |
\item a 6\AA\ radius sphere of Argon in water. |
| 654 |
< |
\end{enumerate} |
| 655 |
< |
|
| 656 |
< |
By utilizing the pairwise techniques (outlined in section |
| 657 |
< |
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 658 |
< |
charged particles, and mixtures of the two, we hope to discern under |
| 659 |
< |
which conditions it will be possible to use one of the alternative |
| 660 |
< |
summation methodologies instead of the Ewald sum. |
| 661 |
< |
|
| 662 |
< |
For the solid and liquid water configurations, configurations were |
| 663 |
< |
taken at regular intervals from high temperature trajectories of 1000 |
| 664 |
< |
SPC/E water molecules. Each configuration was equilibrated |
| 665 |
< |
independently at a lower temperature (300K for the liquid, 200K for |
| 666 |
< |
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 667 |
< |
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 668 |
< |
these systems were selected and equilibrated in the same manner as the |
| 669 |
< |
water systems. In order to introduce measurable fluctuations in the |
| 670 |
< |
configuration energy differences, the crystalline simulations were |
| 671 |
< |
equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 672 |
< |
NaCl configurations needed to represent a fully disordered array of |
| 673 |
< |
point charges, so the high temperature of 7000K was selected for |
| 674 |
< |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 675 |
< |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 676 |
< |
water positions were then randomly swapped, and the resulting |
| 677 |
< |
configurations were again equilibrated individually. Finally, for the |
| 678 |
< |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 679 |
< |
waters within 6\AA\ of the center of the equilibrated water |
| 680 |
< |
configurations were converted to argon. |
| 681 |
< |
|
| 682 |
< |
These procedures guaranteed us a set of representative configurations |
| 683 |
< |
from chemically-relevant systems sampled from appropriate |
| 684 |
< |
ensembles. Force field parameters for the ions and Argon were taken |
| 685 |
< |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 686 |
< |
|
| 687 |
< |
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 688 |
< |
We compared the following alternative summation methods with results |
| 689 |
< |
from the reference method ({\sc spme}): |
| 690 |
< |
|
| 691 |
< |
\begin{enumerate} |
| 692 |
< |
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 693 |
< |
and 0.3\AA$^{-1}$, |
| 694 |
< |
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 695 |
< |
and 0.3\AA$^{-1}$, |
| 696 |
< |
\item reaction field with an infinite dielectric constant, and |
| 697 |
< |
\item an unmodified cutoff. |
| 698 |
< |
\end{enumerate} |
| 699 |
< |
|
| 700 |
< |
Group-based cutoffs with a fifth-order polynomial switching function |
| 701 |
< |
were utilized for the reaction field simulations. Additionally, we |
| 702 |
< |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 703 |
< |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 704 |
< |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 705 |
< |
were performed using the {\sc oopse} molecular mechanics |
| 706 |
< |
package.\cite{Meineke05} All other portions of the energy calculation |
| 707 |
< |
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 708 |
< |
manner across all systems and configurations. |
| 709 |
< |
|
| 710 |
< |
The alternative methods were also evaluated with three different |
| 711 |
< |
cutoff radii (9, 12, and 15\AA). As noted previously, the |
| 712 |
< |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 713 |
< |
real-space and reciprocal-space portions of the Ewald calculation. |
| 714 |
< |
Typical molecular mechanics packages set this to a value dependent on |
| 715 |
< |
the cutoff radius and a tolerance (typically less than $1 \times |
| 716 |
< |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 717 |
< |
increasing accuracy at the expense of computational time spent on the |
| 718 |
< |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 719 |
< |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 720 |
< |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 721 |
< |
0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\ |
| 722 |
< |
respectively. |
| 723 |
< |
|
| 724 |
< |
|
| 725 |
< |
\section{Discussion on the Pairwise Technique Evaluation} |
| 726 |
< |
|
| 727 |
< |
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 728 |
< |
In order to evaluate the performance of the pairwise electrostatic |
| 729 |
< |
summation methods for Monte Carlo simulations, the energy differences |
| 730 |
< |
between configurations were compared to the values obtained when using |
| 731 |
< |
{\sc spme}. The results for the combined regression analysis of all |
| 732 |
< |
of the systems are shown in figure \ref{fig:delE}. |
| 733 |
< |
|
| 734 |
< |
\begin{figure} |
| 735 |
< |
\centering |
| 736 |
< |
\includegraphics[width=4.75in]{./figures/delEplot.pdf} |
| 737 |
< |
\caption{Statistical analysis of the quality of configurational energy |
| 738 |
< |
differences for a given electrostatic method compared with the |
| 739 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 740 |
< |
indicate $\Delta E$ values indistinguishable from those obtained using |
| 741 |
< |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 742 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 743 |
< |
inverted triangles).} |
| 744 |
< |
\label{fig:delE} |
| 745 |
< |
\end{figure} |
| 746 |
< |
|
| 747 |
< |
The most striking feature of this plot is how well the Shifted Force |
| 748 |
< |
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 749 |
< |
differences. For the undamped {\sc sf} method, and the |
| 750 |
< |
moderately-damped {\sc sp} methods, the results are nearly |
| 751 |
< |
indistinguishable from the Ewald results. The other common methods do |
| 752 |
< |
significantly less well. |
| 753 |
< |
|
| 754 |
< |
The unmodified cutoff method is essentially unusable. This is not |
| 755 |
< |
surprising since hard cutoffs give large energy fluctuations as atoms |
| 756 |
< |
or molecules move in and out of the cutoff |
| 757 |
< |
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
| 758 |
< |
some degree by using group based cutoffs with a switching |
| 759 |
< |
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 760 |
< |
significant improvement using the group-switched cutoff because the |
| 761 |
< |
salt and salt solution systems contain non-neutral groups. Section |
| 762 |
< |
\ref{sec:SystemResults} includes results for systems comprised entirely |
| 763 |
< |
of neutral groups. |
| 764 |
< |
|
| 765 |
< |
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 766 |
< |
the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$ |
| 767 |
< |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 768 |
< |
results, particularly with a cutoff radius greater than 12 |
| 769 |
< |
\AA . Use of a larger damping parameter is more helpful for the |
| 770 |
< |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 771 |
< |
with larger cutoffs. |
| 772 |
< |
|
| 773 |
< |
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 774 |
< |
worse} correlation with Ewald. Overall, the undamped case is the best |
| 775 |
< |
performing set, as the correlation and quality of fits are |
| 776 |
< |
consistently superior regardless of the cutoff distance. The undamped |
| 777 |
< |
case is also less computationally demanding (because no evaluation of |
| 778 |
< |
the complementary error function is required). |
| 779 |
< |
|
| 780 |
< |
The reaction field results illustrates some of that method's |
| 781 |
< |
limitations, primarily that it was developed for use in homogenous |
| 782 |
< |
systems; although it does provide results that are an improvement over |
| 783 |
< |
those from an unmodified cutoff. |
| 784 |
< |
|
| 785 |
< |
\sub |
| 786 |
< |
|
| 787 |
< |
\subsection{Magnitudes of the Force and Torque Vectors} |
| 788 |
< |
|
| 789 |
< |
Evaluation of pairwise methods for use in Molecular Dynamics |
| 790 |
< |
simulations requires consideration of effects on the forces and |
| 791 |
< |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
| 792 |
< |
regression results for the force and torque vector magnitudes, |
| 793 |
< |
respectively. The data in these figures was generated from an |
| 794 |
< |
accumulation of the statistics from all of the system types. |
| 795 |
< |
|
| 796 |
< |
\begin{figure} |
| 797 |
< |
\centering |
| 798 |
< |
\includegraphics[width=4.75in]{./figures/frcMagplot.pdf} |
| 799 |
< |
\caption{Statistical analysis of the quality of the force vector |
| 800 |
< |
magnitudes for a given electrostatic method compared with the |
| 801 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 802 |
< |
indicate force magnitude values indistinguishable from those obtained |
| 803 |
< |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 804 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 805 |
< |
inverted triangles).} |
| 806 |
< |
\label{fig:frcMag} |
| 807 |
< |
\end{figure} |
| 808 |
< |
|
| 809 |
< |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 810 |
< |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 811 |
< |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 812 |
< |
This is perhaps expected because this method explicitly incorporates a |
| 813 |
< |
smooth transition in the forces at the cutoff radius as well as the |
| 814 |
< |
neutralizing image charges. |
| 815 |
< |
|
| 816 |
< |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 817 |
< |
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 818 |
< |
poor, but using group based cutoffs and a switching function provides |
| 819 |
< |
an improvement much more significant than what was seen with $\Delta |
| 820 |
< |
E$. |
| 821 |
< |
|
| 822 |
< |
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 823 |
< |
method is generating usable forces. Further increases in damping, |
| 824 |
< |
while beneficial for simulations with a cutoff radius of 9\AA\ , is |
| 825 |
< |
detrimental to simulations with larger cutoff radii. |
| 826 |
< |
|
| 827 |
< |
The reaction field results are surprisingly good, considering the poor |
| 828 |
< |
quality of the fits for the $\Delta E$ results. There is still a |
| 829 |
< |
considerable degree of scatter in the data, but the forces correlate |
| 830 |
< |
well with the Ewald forces in general. We note that the reaction |
| 831 |
< |
field calculations do not include the pure NaCl systems, so these |
| 832 |
< |
results are partly biased towards conditions in which the method |
| 833 |
< |
performs more favorably. |
| 834 |
< |
|
| 835 |
< |
\begin{figure} |
| 836 |
< |
\centering |
| 837 |
< |
\includegraphics[width=4.75in]{./figures/trqMagplot.pdf} |
| 838 |
< |
\caption{Statistical analysis of the quality of the torque vector |
| 839 |
< |
magnitudes for a given electrostatic method compared with the |
| 840 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 841 |
< |
indicate torque magnitude values indistinguishable from those obtained |
| 842 |
< |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 843 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 844 |
< |
inverted triangles).} |
| 845 |
< |
\label{fig:trqMag} |
| 846 |
< |
\end{figure} |
| 46 |
> |
\input{electrostaticsChapter} |
| 47 |
|
|
| 848 |
– |
Molecular torques were only available from the systems which contained |
| 849 |
– |
rigid molecules (i.e. the systems containing water). The data in |
| 850 |
– |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 851 |
– |
|
| 852 |
– |
Torques appear to be much more sensitive to charges at a longer |
| 853 |
– |
distance. The striking feature in comparing the new electrostatic |
| 854 |
– |
methods with {\sc spme} is how much the agreement improves with increasing |
| 855 |
– |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 856 |
– |
appears to be reproducing the {\sc spme} torques most accurately. |
| 857 |
– |
|
| 858 |
– |
Water molecules are dipolar, and the reaction field method reproduces |
| 859 |
– |
the effect of the surrounding polarized medium on each of the |
| 860 |
– |
molecular bodies. Therefore it is not surprising that reaction field |
| 861 |
– |
performs best of all of the methods on molecular torques. |
| 862 |
– |
|
| 863 |
– |
\subsection{Directionality of the Force and Torque Vectors} |
| 864 |
– |
|
| 865 |
– |
It is clearly important that a new electrostatic method can reproduce |
| 866 |
– |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 867 |
– |
sum. However, the {\it directionality} of these vectors will also be |
| 868 |
– |
vital in calculating dynamical quantities accurately. Force and |
| 869 |
– |
torque directionalities were investigated by measuring the angles |
| 870 |
– |
formed between these vectors and the same vectors calculated using |
| 871 |
– |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 872 |
– |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 873 |
– |
distributions of the combined set over all system types. |
| 874 |
– |
|
| 875 |
– |
\begin{figure} |
| 876 |
– |
\centering |
| 877 |
– |
\includegraphics[width=4.75in]{./figures/frcTrqAngplot.pdf} |
| 878 |
– |
\caption{Statistical analysis of the width of the angular distribution |
| 879 |
– |
that the force and torque vectors from a given electrostatic method |
| 880 |
– |
make with their counterparts obtained using the reference Ewald sum. |
| 881 |
– |
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 882 |
– |
indicate force and torque directions indistinguishable from those |
| 883 |
– |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 884 |
– |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 885 |
– |
and 15\AA\ = inverted triangles).} |
| 886 |
– |
\label{fig:frcTrqAng} |
| 887 |
– |
\end{figure} |
| 888 |
– |
|
| 889 |
– |
Both the force and torque $\sigma^2$ results from the analysis of the |
| 890 |
– |
total accumulated system data are tabulated in figure |
| 891 |
– |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 892 |
– |
sp}) method would be essentially unusable for molecular dynamics |
| 893 |
– |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 894 |
– |
method, however, is generating force and torque vectors which are |
| 895 |
– |
within a few degrees of the Ewald results even with weak (or no) |
| 896 |
– |
damping. |
| 897 |
– |
|
| 898 |
– |
All of the sets (aside from the over-damped case) show the improvement |
| 899 |
– |
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 900 |
– |
from 9 to 12\AA\ typically results in a halving of the width of the |
| 901 |
– |
distribution, with a similar improvement when going from 12 to 15 |
| 902 |
– |
\AA . |
| 903 |
– |
|
| 904 |
– |
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 905 |
– |
all do equivalently well at capturing the direction of both the force |
| 906 |
– |
and torque vectors. Using the electrostatic damping improves the |
| 907 |
– |
angular behavior significantly for the {\sc sp} and moderately for the |
| 908 |
– |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 909 |
– |
it is important to recognize that the force vectors cover all |
| 910 |
– |
particles in all seven systems, while torque vectors are only |
| 911 |
– |
available for neutral molecular groups. Damping is more beneficial to |
| 912 |
– |
charged bodies, and this observation is investigated further in |
| 913 |
– |
section \ref{SystemResults}. |
| 914 |
– |
|
| 915 |
– |
Although not discussed previously, group based cutoffs can be applied |
| 916 |
– |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 917 |
– |
will reintroduce small discontinuities at the cutoff radius, but the |
| 918 |
– |
effects of these can be minimized by utilizing a switching function. |
| 919 |
– |
Though there are no significant benefits or drawbacks observed in |
| 920 |
– |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 921 |
– |
is a measurable improvement in the directionality of the forces and |
| 922 |
– |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 923 |
– |
obtained both without (N) and with (Y) group based cutoffs and a |
| 924 |
– |
switching function. Note that the $\alpha$ values have units of |
| 925 |
– |
\AA$^{-1}$ and the variance values have units of degrees$^2$. The |
| 926 |
– |
{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much |
| 927 |
– |
narrower angular distributions when using group-based cutoffs. The |
| 928 |
– |
{\sc sf} method likewise shows improvement in the undamped and lightly |
| 929 |
– |
damped cases. |
| 930 |
– |
|
| 931 |
– |
\begin{table} |
| 932 |
– |
\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$) |
| 933 |
– |
THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A |
| 934 |
– |
GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING |
| 935 |
– |
THE REFERENCE EWALD SUMMATION} |
| 936 |
– |
|
| 937 |
– |
\footnotesize |
| 938 |
– |
\begin{center} |
| 939 |
– |
\begin{tabular}{@{} ccrrrrrrrr @{}} \\ |
| 940 |
– |
\toprule |
| 941 |
– |
\toprule |
| 942 |
– |
|
| 943 |
– |
& &\multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted |
| 944 |
– |
Force} \\ |
| 945 |
– |
\cmidrule(lr){3-6} \cmidrule(l){7-10} $R_\textrm{c}$ & Groups & |
| 946 |
– |
$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & |
| 947 |
– |
$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\ |
| 948 |
– |
|
| 949 |
– |
\midrule |
| 950 |
– |
|
| 951 |
– |
9\AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
| 952 |
– |
& \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\ |
| 953 |
– |
12\AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
| 954 |
– |
& \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\ |
| 955 |
– |
15\AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
| 956 |
– |
& \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\ |
| 957 |
– |
|
| 958 |
– |
\midrule |
| 959 |
– |
|
| 960 |
– |
9\AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
| 961 |
– |
& \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\ |
| 962 |
– |
12\AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
| 963 |
– |
& \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\ |
| 964 |
– |
15\AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
| 965 |
– |
& \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\ |
| 966 |
– |
|
| 967 |
– |
\bottomrule |
| 968 |
– |
\end{tabular} |
| 969 |
– |
\end{center} |
| 970 |
– |
\label{tab:groupAngle} |
| 971 |
– |
\end{table} |
| 972 |
– |
|
| 973 |
– |
One additional trend in table \ref{tab:groupAngle} is that the |
| 974 |
– |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 975 |
– |
increases, something that is more obvious with group-based cutoffs. |
| 976 |
– |
The complimentary error function inserted into the potential weakens |
| 977 |
– |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 978 |
– |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 979 |
– |
electrostatic interaction and to remove it completely. Kast |
| 980 |
– |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 981 |
– |
values for these types of electrostatic summation methods by fitting |
| 982 |
– |
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 983 |
– |
0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\ |
| 984 |
– |
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 985 |
– |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 986 |
– |
these findings, choices this high would introduce error in the |
| 987 |
– |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 988 |
– |
observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial, |
| 989 |
– |
but damping may be unnecessary when using the {\sc sf} method. |
| 990 |
– |
|
| 991 |
– |
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 992 |
– |
|
| 993 |
– |
Zahn {\it et al.} investigated the structure and dynamics of water |
| 994 |
– |
using eqs. (\ref{eq:ZahnPot}) and |
| 995 |
– |
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 996 |
– |
that a method similar (but not identical with) the damped {\sc sf} |
| 997 |
– |
method resulted in properties very similar to those obtained when |
| 998 |
– |
using the Ewald summation. The properties they studied (pair |
| 999 |
– |
distribution functions, diffusion constants, and velocity and |
| 1000 |
– |
orientational correlation functions) may not be particularly sensitive |
| 1001 |
– |
to the long-range and collective behavior that governs the |
| 1002 |
– |
low-frequency behavior in crystalline systems. Additionally, the |
| 1003 |
– |
ionic crystals are the worst case scenario for the pairwise methods |
| 1004 |
– |
because they lack the reciprocal space contribution contained in the |
| 1005 |
– |
Ewald summation. |
| 1006 |
– |
|
| 1007 |
– |
We are using two separate measures to probe the effects of these |
| 1008 |
– |
alternative electrostatic methods on the dynamics in crystalline |
| 1009 |
– |
materials. For short- and intermediate-time dynamics, we are |
| 1010 |
– |
computing the velocity autocorrelation function, and for long-time |
| 1011 |
– |
and large length-scale collective motions, we are looking at the |
| 1012 |
– |
low-frequency portion of the power spectrum. |
| 1013 |
– |
|
| 1014 |
– |
\begin{figure} |
| 1015 |
– |
\centering |
| 1016 |
– |
\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf} |
| 1017 |
– |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1018 |
– |
1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1019 |
– |
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 1020 |
– |
the first minimum. The times to first collision are nearly identical, |
| 1021 |
– |
but differences can be seen in the peaks and troughs, where the |
| 1022 |
– |
undamped and weakly damped methods are stiffer than the moderately |
| 1023 |
– |
damped and {\sc spme} methods.} |
| 1024 |
– |
\label{fig:vCorrPlot} |
| 1025 |
– |
\end{figure} |
| 1026 |
– |
|
| 1027 |
– |
The short-time decay of the velocity autocorrelation function through |
| 1028 |
– |
the first collision are nearly identical in figure |
| 1029 |
– |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1030 |
– |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1031 |
– |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1032 |
– |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1033 |
– |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1034 |
– |
approaches the {\sc spme} results. With $\alpha$ values of 0.2\AA$^{-1}$, |
| 1035 |
– |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1036 |
– |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1037 |
– |
and {\sc sp} potentials become nearly identical with increased |
| 1038 |
– |
damping. However, this appears to indicate that once damping is |
| 1039 |
– |
utilized, the details of the form of the potential (and forces) |
| 1040 |
– |
constructed out of the damped electrostatic interaction are less |
| 1041 |
– |
important. |
| 1042 |
– |
|
| 1043 |
– |
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 1044 |
– |
|
| 1045 |
– |
To evaluate how the differences between the methods affect the |
| 1046 |
– |
collective long-time motion, we computed power spectra from long-time |
| 1047 |
– |
traces of the velocity autocorrelation function. The power spectra for |
| 1048 |
– |
the best-performing alternative methods are shown in |
| 1049 |
– |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1050 |
– |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1051 |
– |
ringing resulting from data truncation. This procedure had no |
| 1052 |
– |
noticeable effect on peak location or magnitude. |
| 1053 |
– |
|
| 1054 |
– |
\begin{figure} |
| 1055 |
– |
\centering |
| 1056 |
– |
\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf} |
| 1057 |
– |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1058 |
– |
functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf} |
| 1059 |
– |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
| 1060 |
– |
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
| 1061 |
– |
spectra differ.} |
| 1062 |
– |
\label{fig:methodPS} |
| 1063 |
– |
\end{figure} |
| 1064 |
– |
|
| 1065 |
– |
While the high frequency regions of the power spectra for the |
| 1066 |
– |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1067 |
– |
the low frequency region shows how the summation methods differ. |
| 1068 |
– |
Considering the low-frequency inset (expanded in the upper frame of |
| 1069 |
– |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1070 |
– |
correlated motions are blue-shifted when using undamped or weakly |
| 1071 |
– |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1072 |
– |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1073 |
– |
correlated motion to the Ewald method (which has a convergence |
| 1074 |
– |
parameter of 0.3119\AA$^{-1}$). This weakening of the electrostatic |
| 1075 |
– |
interaction with increased damping explains why the long-ranged |
| 1076 |
– |
correlated motions are at lower frequencies for the moderately damped |
| 1077 |
– |
methods than for undamped or weakly damped methods. |
| 1078 |
– |
|
| 1079 |
– |
To isolate the role of the damping constant, we have computed the |
| 1080 |
– |
spectra for a single method ({\sc sf}) with a range of damping |
| 1081 |
– |
constants and compared this with the {\sc spme} spectrum. |
| 1082 |
– |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1083 |
– |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1084 |
– |
However, even without any electrostatic damping, the {\sc sf} method |
| 1085 |
– |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1086 |
– |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1087 |
– |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1088 |
– |
Most} of the collective behavior in the crystal is accurately captured |
| 1089 |
– |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1090 |
– |
obtained using moderate damping in addition to the shifting at the |
| 1091 |
– |
cutoff distance. |
| 1092 |
– |
|
| 1093 |
– |
\begin{figure} |
| 1094 |
– |
\centering |
| 1095 |
– |
\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf} |
| 1096 |
– |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1097 |
– |
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1098 |
– |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1099 |
– |
electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement |
| 1100 |
– |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1101 |
– |
result in underestimates of frequencies of the long-wavelength |
| 1102 |
– |
motions.} |
| 1103 |
– |
\label{fig:dampInc} |
| 1104 |
– |
\end{figure} |
| 1105 |
– |
|
| 1106 |
– |
|
| 48 |
|
\chapter{\label{chap:water}SIMPLE MODELS FOR WATER} |
| 49 |
|
|
| 50 |
|
\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS} |
| 57 |
|
\backmatter |
| 58 |
|
|
| 59 |
|
\bibliographystyle{ndthesis} |
| 60 |
< |
\bibliography{dissertation} |
| 60 |
> |
\bibliography{dissertation} |
| 61 |
|
|
| 62 |
|
\end{document} |
| 63 |
|
|
