| 43 |
|
|
| 44 |
|
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} |
| 45 |
|
|
| 46 |
< |
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION |
| 47 |
< |
TECHNIQUES} |
| 48 |
< |
|
| 49 |
< |
In molecular simulations, proper accumulation of the electrostatic |
| 50 |
< |
interactions is essential and is one of the most |
| 51 |
< |
computationally-demanding tasks. The common molecular mechanics force |
| 52 |
< |
fields represent atomic sites with full or partial charges protected |
| 53 |
< |
by Lennard-Jones (short range) interactions. This means that nearly |
| 54 |
< |
every pair interaction involves a calculation of charge-charge forces. |
| 55 |
< |
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
| 56 |
< |
interactions quickly become the most expensive part of molecular |
| 57 |
< |
simulations. Historically, the electrostatic pair interaction would |
| 58 |
< |
not have decayed appreciably within the typical box lengths that could |
| 59 |
< |
be feasibly simulated. In the larger systems that are more typical of |
| 60 |
< |
modern simulations, large cutoffs should be used to incorporate |
| 61 |
< |
electrostatics correctly. |
| 62 |
< |
|
| 63 |
< |
There have been many efforts to address the proper and practical |
| 64 |
< |
handling of electrostatic interactions, and these have resulted in a |
| 65 |
< |
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are |
| 66 |
< |
typically classified as implicit methods (i.e., continuum dielectrics, |
| 67 |
< |
static dipolar fields),\cite{Born20,Grossfield00} explicit methods |
| 68 |
< |
(i.e., Ewald summations, interaction shifting or |
| 69 |
< |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 70 |
< |
reaction field type methods, fast multipole |
| 71 |
< |
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 72 |
< |
often preferred because they physically incorporate solvent molecules |
| 73 |
< |
in the system of interest, but these methods are sometimes difficult |
| 74 |
< |
to utilize because of their high computational cost.\cite{Roux99} In |
| 75 |
< |
addition to the computational cost, there have been some questions |
| 76 |
< |
regarding possible artifacts caused by the inherent periodicity of the |
| 77 |
< |
explicit Ewald summation.\cite{Tobias01} |
| 78 |
< |
|
| 79 |
< |
In this chapter, we focus on a new set of pairwise methods devised by |
| 80 |
< |
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 81 |
< |
methods along with a few other mixed methods (i.e. reaction field) are |
| 82 |
< |
compared with the smooth particle mesh Ewald |
| 83 |
< |
sum,\cite{Onsager36,Essmann99} which is our reference method for |
| 84 |
< |
handling long-range electrostatic interactions. The new methods for |
| 85 |
< |
handling electrostatics have the potential to scale linearly with |
| 86 |
< |
increasing system size since they involve only a simple modification |
| 87 |
< |
to the direct pairwise sum. They also lack the added periodicity of |
| 88 |
< |
the Ewald sum, so they can be used for systems which are non-periodic |
| 89 |
< |
or which have one- or two-dimensional periodicity. Below, these |
| 90 |
< |
methods are evaluated using a variety of model systems to |
| 91 |
< |
establish their usability in molecular simulations. |
| 92 |
< |
|
| 93 |
< |
\section{The Ewald Sum} |
| 94 |
< |
|
| 95 |
< |
The complete accumulation of the electrostatic interactions in a system with |
| 96 |
< |
periodic boundary conditions (PBC) requires the consideration of the |
| 97 |
< |
effect of all charges within a (cubic) simulation box as well as those |
| 98 |
< |
in the periodic replicas, |
| 99 |
< |
\begin{equation} |
| 100 |
< |
V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime |
| 101 |
< |
\left[ \sum_{i=1}^N\sum_{j=1}^N \phi |
| 102 |
< |
\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) |
| 103 |
< |
\right], |
| 104 |
< |
\label{eq:PBCSum} |
| 105 |
< |
\end{equation} |
| 106 |
< |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
| 107 |
< |
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
| 108 |
< |
prime indicates $i = j$ are neglected for $\mathbf{n} = |
| 109 |
< |
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
| 110 |
< |
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
| 111 |
< |
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
| 112 |
< |
$j$, and $\phi$ is the solution to Poisson's equation |
| 113 |
< |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for |
| 114 |
< |
charge-charge interactions). In the case of monopole electrostatics, |
| 115 |
< |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for |
| 116 |
< |
non-neutral systems. |
| 117 |
< |
|
| 118 |
< |
The electrostatic summation problem was originally studied by Ewald |
| 119 |
< |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 120 |
< |
took was to convert this conditionally convergent sum into two |
| 121 |
< |
absolutely convergent summations: a short-ranged real-space summation |
| 122 |
< |
and a long-ranged reciprocal-space summation, |
| 123 |
< |
\begin{equation} |
| 124 |
< |
\begin{split} |
| 125 |
< |
V_\textrm{elec} = \frac{1}{2}& |
| 126 |
< |
\sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime |
| 127 |
< |
\frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)} |
| 128 |
< |
{|\mathbf{r}_{ij}+\mathbf{n}|} \\ |
| 129 |
< |
&+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} |
| 130 |
< |
\exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) |
| 131 |
< |
\cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ |
| 132 |
< |
&- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 |
| 133 |
< |
+ \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3} |
| 134 |
< |
\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
| 135 |
< |
\end{split} |
| 136 |
< |
\label{eq:EwaldSum} |
| 137 |
< |
\end{equation} |
| 138 |
< |
where $\alpha$ is the damping or convergence parameter with units of |
| 139 |
< |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 140 |
< |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 141 |
< |
constant of the surrounding medium. The final two terms of |
| 142 |
< |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 143 |
< |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 144 |
< |
dipolar term was neglected in early applications in molecular |
| 145 |
< |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 146 |
< |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 147 |
< |
dipole moment which is magnified through replication of the periodic |
| 148 |
< |
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
| 149 |
< |
system is said to be using conducting (or ``tin-foil'') boundary |
| 150 |
< |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 151 |
< |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 152 |
< |
time. Initially, due to the small system sizes that could be |
| 153 |
< |
simulated feasibly, the entire simulation box was replicated to |
| 154 |
< |
convergence. In more modern simulations, the systems have grown large |
| 155 |
< |
enough that a real-space cutoff could potentially give convergent |
| 156 |
< |
behavior. Indeed, it has been observed that with the choice of a |
| 157 |
< |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 158 |
< |
rapidly convergent and small relative to the real-space |
| 159 |
< |
portion.\cite{Karasawa89,Kolafa92} |
| 160 |
< |
|
| 161 |
< |
\begin{figure} |
| 162 |
< |
\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf} |
| 163 |
< |
\caption{The change in the need for the Ewald sum with |
| 164 |
< |
increasing computational power. A:~Initially, only small systems |
| 165 |
< |
could be studied, and the Ewald sum replicated the simulation box to |
| 166 |
< |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 167 |
< |
convergence for the larger systems of charges that are common today.} |
| 168 |
< |
\label{fig:ewaldTime} |
| 169 |
< |
\end{figure} |
| 170 |
< |
|
| 171 |
< |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 172 |
< |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 173 |
< |
the computational cost between the direct and reciprocal-space |
| 174 |
< |
portions of the summation. The choice of this value allows one to |
| 175 |
< |
select whether the real-space or reciprocal space portion of the |
| 176 |
< |
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
| 177 |
< |
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 178 |
< |
$\alpha$ and thoughtful algorithm development, this cost can be |
| 179 |
< |
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 180 |
< |
taken to reduce the cost of the Ewald summation even further is to set |
| 181 |
< |
$\alpha$ such that the real-space interactions decay rapidly, allowing |
| 182 |
< |
for a short spherical cutoff. Then the reciprocal space summation is |
| 183 |
< |
optimized. These optimizations usually involve utilization of the |
| 184 |
< |
fast Fourier transform (FFT),\cite{Hockney81} leading to the |
| 185 |
< |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 186 |
< |
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 187 |
< |
methods, the cost of the reciprocal-space portion of the Ewald |
| 188 |
< |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
| 189 |
< |
\log N)$. |
| 190 |
< |
|
| 191 |
< |
These developments and optimizations have made the use of the Ewald |
| 192 |
< |
summation routine in simulations with periodic boundary |
| 193 |
< |
conditions. However, in certain systems, such as vapor-liquid |
| 194 |
< |
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 195 |
< |
can prove problematic. The Ewald sum has been reformulated to handle |
| 196 |
< |
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these |
| 197 |
< |
methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 198 |
< |
recently, there have been several successful efforts toward reducing |
| 199 |
< |
the computational cost of 2-D lattice |
| 200 |
< |
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 201 |
< |
bringing them more in line with the cost of the full 3-D summation. |
| 202 |
< |
|
| 203 |
< |
Several studies have recognized that the inherent periodicity in the |
| 204 |
< |
Ewald sum can also have an effect on three-dimensional |
| 205 |
< |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 206 |
< |
Solvated proteins are essentially kept at high concentration due to |
| 207 |
< |
the periodicity of the electrostatic summation method. In these |
| 208 |
< |
systems, the more compact folded states of a protein can be |
| 209 |
< |
artificially stabilized by the periodic replicas introduced by the |
| 210 |
< |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 211 |
< |
considering the use of the Ewald summation where the assumed |
| 212 |
< |
periodicity would introduce spurious effects in the system dynamics. |
| 213 |
< |
|
| 214 |
< |
|
| 215 |
< |
\section{The Wolf and Zahn Methods} |
| 216 |
< |
|
| 217 |
< |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 218 |
< |
for the accurate accumulation of electrostatic interactions in an |
| 219 |
< |
efficient pairwise fashion. This procedure lacks the inherent |
| 220 |
< |
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
| 221 |
< |
observed that the electrostatic interaction is effectively |
| 222 |
< |
short-ranged in condensed phase systems and that neutralization of the |
| 223 |
< |
charge contained within the cutoff radius is crucial for potential |
| 224 |
< |
stability. They devised a pairwise summation method that ensures |
| 225 |
< |
charge neutrality and gives results similar to those obtained with the |
| 226 |
< |
Ewald summation. The resulting shifted Coulomb potential includes |
| 227 |
< |
image-charges subtracted out through placement on the cutoff sphere |
| 228 |
< |
and a distance-dependent damping function (identical to that seen in |
| 229 |
< |
the real-space portion of the Ewald sum) to aid convergence |
| 230 |
< |
\begin{equation} |
| 231 |
< |
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}} |
| 232 |
< |
- \lim_{r_{ij}\rightarrow R_\textrm{c}} |
| 233 |
< |
\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 234 |
< |
\label{eq:WolfPot} |
| 235 |
< |
\end{equation} |
| 236 |
< |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 237 |
< |
potential. However, neutralizing the charge contained within each |
| 238 |
< |
cutoff sphere requires the placement of a self-image charge on the |
| 239 |
< |
surface of the cutoff sphere. This additional self-term in the total |
| 240 |
< |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
| 241 |
< |
Madelung energies for many crystals. |
| 242 |
< |
|
| 243 |
< |
In order to use their charge-neutralized potential in molecular |
| 244 |
< |
dynamics simulations, Wolf \textit{et al.} suggested taking the |
| 245 |
< |
derivative of this potential prior to evaluation of the limit. This |
| 246 |
< |
procedure gives an expression for the forces, |
| 247 |
< |
\begin{equation} |
| 248 |
< |
\begin{split} |
| 249 |
< |
F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{& |
| 250 |
< |
\Biggr[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}} |
| 251 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}} |
| 252 |
< |
\Biggr]\\ |
| 253 |
< |
&-\Biggr[ |
| 254 |
< |
\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2} |
| 255 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 256 |
< |
\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} |
| 257 |
< |
\Biggr]\Biggr\}, |
| 258 |
< |
\end{split} |
| 259 |
< |
\label{eq:WolfForces} |
| 260 |
< |
\end{equation} |
| 261 |
< |
that incorporates both image charges and damping of the electrostatic |
| 262 |
< |
interaction. |
| 263 |
< |
|
| 264 |
< |
More recently, Zahn \textit{et al.} investigated these potential and |
| 265 |
< |
force expressions for use in simulations involving water.\cite{Zahn02} |
| 266 |
< |
In their work, they pointed out that the forces and derivative of |
| 267 |
< |
the potential are not commensurate. Attempts to use both |
| 268 |
< |
eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 269 |
< |
to poor energy conservation. They correctly observed that taking the |
| 270 |
< |
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 271 |
< |
derivatives gives forces for a different potential energy function |
| 272 |
< |
than the one shown in eq. (\ref{eq:WolfPot}). |
| 273 |
< |
|
| 274 |
< |
Zahn \textit{et al.} introduced a modified form of this summation |
| 275 |
< |
method as a way to use the technique in Molecular Dynamics |
| 276 |
< |
simulations. They proposed a new damped Coulomb potential, |
| 277 |
< |
\begin{equation} |
| 278 |
< |
\begin{split} |
| 279 |
< |
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{& |
| 280 |
< |
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} \\ |
| 281 |
< |
&- \left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
| 282 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 283 |
< |
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
| 284 |
< |
\right]\left(r_{ij}-R_\mathrm{c}\right)\Biggr\}, |
| 285 |
< |
\end{split} |
| 286 |
< |
\label{eq:ZahnPot} |
| 287 |
< |
\end{equation} |
| 288 |
< |
and showed that this potential does fairly well at capturing the |
| 289 |
< |
structural and dynamic properties of water compared the same |
| 290 |
< |
properties obtained using the Ewald sum. |
| 291 |
< |
|
| 292 |
< |
\section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation} |
| 293 |
< |
|
| 294 |
< |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 295 |
< |
al.} are constructed using two different (and separable) computational |
| 296 |
< |
tricks: |
| 297 |
< |
|
| 298 |
< |
\begin{enumerate}[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 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 |
| 312 |
< |
than $r^{-1}$, either the cutoff sphere neutralization or the |
| 313 |
< |
distance-dependent damping technique could be used as a foundation for |
| 314 |
< |
a new pairwise summation method. Wolf \textit{et al.} made the |
| 315 |
< |
observation that charge neutralization within the cutoff sphere plays |
| 316 |
< |
a significant role in energy convergence; therefore we will begin our |
| 317 |
< |
analysis with the various shifted forms that maintain this charge |
| 318 |
< |
neutralization. We can evaluate the methods of Wolf |
| 319 |
< |
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 320 |
< |
shifted potential, |
| 321 |
< |
\begin{equation} |
| 322 |
< |
V_\textrm{SP}(r) = \begin{cases} |
| 323 |
< |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 324 |
< |
R_\textrm{c} |
| 325 |
< |
\end{cases}, |
| 326 |
< |
\label{eq:shiftingPotForm} |
| 327 |
< |
\end{equation} |
| 328 |
< |
and shifted force, |
| 329 |
< |
\begin{equation} |
| 330 |
< |
V_\textrm{SF}(r) = \begin{cases} |
| 331 |
< |
v(r) - v_\textrm{c} |
| 332 |
< |
- \left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 333 |
< |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 334 |
< |
\end{cases}, |
| 335 |
< |
\label{eq:shiftingForm} |
| 336 |
< |
\end{equation} |
| 337 |
< |
functions where $v(r)$ is the unshifted form of the potential, and |
| 338 |
< |
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
| 339 |
< |
that both the potential and the forces goes to zero at the cutoff |
| 340 |
< |
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
| 341 |
< |
potential is smooth at the cutoff radius |
| 342 |
< |
($R_\textrm{c}$).\cite{Allen87} |
| 343 |
< |
|
| 344 |
< |
The forces associated with the shifted potential are simply the forces |
| 345 |
< |
of the unshifted potential itself (when inside the cutoff sphere), |
| 346 |
< |
\begin{equation} |
| 347 |
< |
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 348 |
< |
\end{equation} |
| 349 |
< |
and are zero outside. Inside the cutoff sphere, the forces associated |
| 350 |
< |
with the shifted force form can be written, |
| 351 |
< |
\begin{equation} |
| 352 |
< |
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 353 |
< |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 354 |
< |
\end{equation} |
| 355 |
< |
|
| 356 |
< |
If the potential, $v(r)$, is taken to be the normal Coulomb potential, |
| 357 |
< |
\begin{equation} |
| 358 |
< |
v(r) = \frac{q_i q_j}{r}, |
| 359 |
< |
\label{eq:Coulomb} |
| 360 |
< |
\end{equation} |
| 361 |
< |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 362 |
< |
al.}'s undamped prescription: |
| 363 |
< |
\begin{equation} |
| 364 |
< |
V_\textrm{SP}(r) = |
| 365 |
< |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 366 |
< |
r\leqslant R_\textrm{c}, |
| 367 |
< |
\label{eq:SPPot} |
| 368 |
< |
\end{equation} |
| 369 |
< |
with associated forces, |
| 370 |
< |
\begin{equation} |
| 371 |
< |
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) |
| 372 |
< |
\quad r\leqslant R_\textrm{c}. |
| 373 |
< |
\label{eq:SPForces} |
| 374 |
< |
\end{equation} |
| 375 |
< |
These forces are identical to the forces of the standard Coulomb |
| 376 |
< |
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 377 |
< |
\textit{et al.} as undesirable. They pointed out that the effect of |
| 378 |
< |
the image charges is neglected in the forces when this form is |
| 379 |
< |
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
| 380 |
< |
molecular dynamics. Additionally, there is a discontinuity in the |
| 381 |
< |
forces at the cutoff radius which results in energy drift during MD |
| 382 |
< |
simulations. |
| 383 |
< |
|
| 384 |
< |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 385 |
< |
will give, |
| 386 |
< |
\begin{equation} |
| 387 |
< |
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}} |
| 388 |
< |
+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] |
| 389 |
< |
\quad r\leqslant R_\textrm{c}. |
| 390 |
< |
\label{eq:SFPot} |
| 391 |
< |
\end{equation} |
| 392 |
< |
with associated forces, |
| 393 |
< |
\begin{equation} |
| 394 |
< |
F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) |
| 395 |
< |
\quad r\leqslant R_\textrm{c}. |
| 396 |
< |
\label{eq:SFForces} |
| 397 |
< |
\end{equation} |
| 398 |
< |
This formulation has the benefits that there are no discontinuities at |
| 399 |
< |
the cutoff radius, while the neutralizing image charges are present in |
| 400 |
< |
both the energy and force expressions. It would be simple to add the |
| 401 |
< |
self-neutralizing term back when computing the total energy of the |
| 402 |
< |
system, thereby maintaining the agreement with the Madelung energies. |
| 403 |
< |
A side effect of this treatment is the alteration in the shape of the |
| 404 |
< |
potential that comes from the derivative term. Thus, a degree of |
| 405 |
< |
clarity about agreement with the empirical potential is lost in order |
| 406 |
< |
to gain functionality in dynamics simulations. |
| 407 |
< |
|
| 408 |
< |
Wolf \textit{et al.} originally discussed the energetics of the |
| 409 |
< |
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
| 410 |
< |
insufficient for accurate determination of the energy with reasonable |
| 411 |
< |
cutoff distances. The calculated Madelung energies fluctuated around |
| 412 |
< |
the expected value as the cutoff radius was increased, but the |
| 413 |
< |
oscillations converged toward the correct value.\cite{Wolf99} A |
| 414 |
< |
damping function was incorporated to accelerate the convergence; and |
| 415 |
< |
though alternative forms for the damping function could be |
| 416 |
< |
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 417 |
< |
chosen to mirror the effective screening used in the Ewald summation. |
| 418 |
< |
Incorporating this error function damping into the simple Coulomb |
| 419 |
< |
potential, |
| 420 |
< |
\begin{equation} |
| 421 |
< |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 422 |
< |
\label{eq:dampCoulomb} |
| 423 |
< |
\end{equation} |
| 424 |
< |
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 425 |
< |
\begin{equation} |
| 426 |
< |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} |
| 427 |
< |
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) |
| 428 |
< |
\quad r\leqslant R_\textrm{c}, |
| 429 |
< |
\label{eq:DSPPot} |
| 430 |
< |
\end{equation} |
| 431 |
< |
with associated forces, |
| 432 |
< |
\begin{equation} |
| 433 |
< |
F_{\textrm{DSP}}(r) = q_iq_j |
| 434 |
< |
\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} |
| 435 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) |
| 436 |
< |
\quad r\leqslant R_\textrm{c}. |
| 437 |
< |
\label{eq:DSPForces} |
| 438 |
< |
\end{equation} |
| 439 |
< |
Again, this damped shifted potential suffers from a |
| 440 |
< |
force-discontinuity at the cutoff radius, and the image charges play |
| 441 |
< |
no role in the forces. To remedy these concerns, one may derive a |
| 442 |
< |
{\sc sf} variant by including the derivative term in |
| 443 |
< |
eq. (\ref{eq:shiftingForm}), |
| 444 |
< |
\begin{equation} |
| 445 |
< |
\begin{split} |
| 446 |
< |
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 447 |
< |
\frac{\mathrm{erfc}\left(\alpha r\right)}{r} |
| 448 |
< |
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ |
| 449 |
< |
&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
| 450 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 451 |
< |
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
| 452 |
< |
\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} |
| 453 |
< |
\quad r\leqslant R_\textrm{c}. |
| 454 |
< |
\label{eq:DSFPot} |
| 455 |
< |
\end{split} |
| 456 |
< |
\end{equation} |
| 457 |
< |
The derivative of the above potential will lead to the following forces, |
| 458 |
< |
\begin{equation} |
| 459 |
< |
\begin{split} |
| 460 |
< |
F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 461 |
< |
\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} |
| 462 |
< |
+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ |
| 463 |
< |
&- \left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)} |
| 464 |
< |
{R_{\textrm{c}}^2} |
| 465 |
< |
+ \frac{2\alpha}{\pi^{1/2}} |
| 466 |
< |
\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} |
| 467 |
< |
\right)\Biggr{]} |
| 468 |
< |
\quad r\leqslant R_\textrm{c}. |
| 469 |
< |
\label{eq:DSFForces} |
| 470 |
< |
\end{split} |
| 471 |
< |
\end{equation} |
| 472 |
< |
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
| 473 |
< |
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
| 474 |
< |
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
| 475 |
< |
|
| 476 |
< |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 477 |
< |
derived by Zahn \textit{et al.}; however, there are two important |
| 478 |
< |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 479 |
< |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 480 |
< |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 481 |
< |
in the Zahn potential, resulting in a potential discontinuity as |
| 482 |
< |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 483 |
< |
portion is different. The missing $v_\textrm{c}$ term would not |
| 484 |
< |
affect molecular dynamics simulations (although the computed energy |
| 485 |
< |
would be expected to have sudden jumps as particle distances crossed |
| 486 |
< |
$R_c$). The sign problem is a potential source of errors, however. |
| 487 |
< |
In fact, it introduces a discontinuity in the forces at the cutoff, |
| 488 |
< |
because the force function is shifted in the wrong direction and |
| 489 |
< |
doesn't cross zero at $R_\textrm{c}$. |
| 490 |
< |
|
| 491 |
< |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 492 |
< |
electrostatic summation method in which the potential and forces are |
| 493 |
< |
continuous at the cutoff radius and which incorporates the damping |
| 494 |
< |
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
| 495 |
< |
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
| 496 |
< |
{\sc sf}, damping) are at reproducing the correct electrostatic |
| 497 |
< |
summation performed by the Ewald sum. |
| 498 |
< |
|
| 499 |
< |
|
| 500 |
< |
\section{Evaluating Pairwise Summation Techniques} |
| 501 |
< |
|
| 502 |
< |
In classical molecular mechanics simulations, there are two primary |
| 503 |
< |
techniques utilized to obtain information about the system of |
| 504 |
< |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 505 |
< |
techniques utilize pairwise summations of interactions between |
| 506 |
< |
particle sites, but they use these summations in different ways. |
| 507 |
< |
|
| 508 |
< |
In MC, the potential energy difference between configurations dictates |
| 509 |
< |
the progression of MC sampling. Going back to the origins of this |
| 510 |
< |
method, the acceptance criterion for the canonical ensemble laid out |
| 511 |
< |
by Metropolis \textit{et al.} states that a subsequent configuration |
| 512 |
< |
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
| 513 |
< |
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
| 514 |
< |
Maintaining the correct $\Delta E$ when using an alternate method for |
| 515 |
< |
handling the long-range electrostatics will ensure proper sampling |
| 516 |
< |
from the ensemble. |
| 517 |
< |
|
| 518 |
< |
In MD, the derivative of the potential governs how the system will |
| 519 |
< |
progress in time. Consequently, the force and torque vectors on each |
| 520 |
< |
body in the system dictate how the system evolves. If the magnitude |
| 521 |
< |
and direction of these vectors are similar when using alternate |
| 522 |
< |
electrostatic summation techniques, the dynamics in the short term |
| 523 |
< |
will be indistinguishable. Because error in MD calculations is |
| 524 |
< |
cumulative, one should expect greater deviation at longer times, |
| 525 |
< |
although methods which have large differences in the force and torque |
| 526 |
< |
vectors will diverge from each other more rapidly. |
| 527 |
< |
|
| 528 |
< |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 529 |
< |
|
| 530 |
< |
The pairwise summation techniques (outlined in section |
| 531 |
< |
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 532 |
< |
studying the energy differences between conformations. We took the |
| 533 |
< |
{\sc spme}-computed energy difference between two conformations to be the |
| 534 |
< |
correct behavior. An ideal performance by an alternative method would |
| 535 |
< |
reproduce these energy differences exactly (even if the absolute |
| 536 |
< |
energies calculated by the methods are different). Since none of the |
| 537 |
< |
methods provide exact energy differences, we used linear least squares |
| 538 |
< |
regressions of energy gap data to evaluate how closely the methods |
| 539 |
< |
mimicked the Ewald energy gaps. Unitary results for both the |
| 540 |
< |
correlation (slope) and correlation coefficient for these regressions |
| 541 |
< |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 542 |
< |
Sample correlation plots for two alternate methods are shown in |
| 543 |
< |
Fig. \ref{fig:linearFit}. |
| 544 |
< |
|
| 545 |
< |
\begin{figure} |
| 546 |
< |
\centering |
| 547 |
< |
\includegraphics[width = \linewidth]{./figures/dualLinear.pdf} |
| 548 |
< |
\caption{Example least squares regressions of the configuration energy |
| 549 |
< |
differences for SPC/E water systems. The upper plot shows a data set |
| 550 |
< |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 551 |
< |
shows a data set with a good correlation coefficient.} |
| 552 |
< |
\label{fig:linearFit} |
| 553 |
< |
\end{figure} |
| 554 |
< |
|
| 555 |
< |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 556 |
< |
were represented using 500 independent configurations. Thus, each of |
| 557 |
< |
the alternative (non-Ewald) electrostatic summation methods was |
| 558 |
< |
evaluated using an accumulated 873,250 configurational energy |
| 559 |
< |
differences. |
| 560 |
< |
|
| 561 |
< |
Results and discussion for the individual analysis of each of the |
| 562 |
< |
system types appear in sections \ref{sec:IndividualResults}, while the |
| 563 |
< |
cumulative results over all the investigated systems appear below in |
| 564 |
< |
sections \ref{sec:EnergyResults}. |
| 565 |
< |
|
| 566 |
< |
\subsection{Molecular Dynamics and the Force and Torque |
| 567 |
< |
Vectors}\label{sec:MDMethods} We evaluated the pairwise methods |
| 568 |
< |
(outlined in section \ref{sec:ESMethods}) for use in MD simulations by |
| 569 |
< |
comparing the force and torque vectors with those obtained using the |
| 570 |
< |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 571 |
< |
direction of these vectors on each of the bodies in the system were |
| 572 |
< |
analyzed. For the magnitude of these vectors, linear least squares |
| 573 |
< |
regression analyses were performed as described previously for |
| 574 |
< |
comparing $\Delta E$ values. Instead of a single energy difference |
| 575 |
< |
between two system configurations, we compared the magnitudes of the |
| 576 |
< |
forces (and torques) on each molecule in each configuration. For a |
| 577 |
< |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 578 |
< |
vectors and 1000 torque vectors. With 500 configurations, this |
| 579 |
< |
results in 520,000 force and 500,000 torque vector comparisons. |
| 580 |
< |
Additionally, data from seven different system types was aggregated |
| 581 |
< |
before the comparison was made. |
| 582 |
< |
|
| 583 |
< |
The {\it directionality} of the force and torque vectors was |
| 584 |
< |
investigated through measurement of the angle ($\theta$) formed |
| 585 |
< |
between those computed from the particular method and those from {\sc spme}, |
| 586 |
< |
\begin{equation} |
| 587 |
< |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} |
| 588 |
< |
\cdot \hat{F}_\textrm{M}\right), |
| 589 |
< |
\end{equation} |
| 590 |
< |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 591 |
< |
vector computed using method M. Each of these $\theta$ values was |
| 592 |
< |
accumulated in a distribution function and weighted by the area on the |
| 593 |
< |
unit sphere. Since this distribution is a measure of angular error |
| 594 |
< |
between two different electrostatic summation methods, there is no |
| 595 |
< |
{\it a priori} reason for the profile to adhere to any specific |
| 596 |
< |
shape. Thus, gaussian fits were used to measure the width of the |
| 597 |
< |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 598 |
< |
each of these fits and was used to compare distribution widths. |
| 599 |
< |
Values of $\sigma^2$ near zero indicate vector directions |
| 600 |
< |
indistinguishable from those calculated when using the reference |
| 601 |
< |
method ({\sc spme}). |
| 602 |
< |
|
| 603 |
< |
\subsection{Short-time Dynamics} |
| 604 |
< |
|
| 605 |
< |
The effects of the alternative electrostatic summation methods on the |
| 606 |
< |
short-time dynamics of charged systems were evaluated by considering a |
| 607 |
< |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 608 |
< |
performing pairwise methods was used in this comparison. The NaCl |
| 609 |
< |
crystal was chosen to avoid possible complications from the treatment |
| 610 |
< |
of orientational motion in molecular systems. All systems were |
| 611 |
< |
started with the same initial positions and velocities. Simulations |
| 612 |
< |
were performed under the microcanonical ensemble, and velocity |
| 613 |
< |
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 614 |
< |
of the trajectories, |
| 615 |
< |
\begin{equation} |
| 616 |
< |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 617 |
< |
\label{eq:vCorr} |
| 618 |
< |
\end{equation} |
| 619 |
< |
Velocity autocorrelation functions require detailed short time data, |
| 620 |
< |
thus velocity information was saved every 2 fs over 10 ps |
| 621 |
< |
trajectories. Because the NaCl crystal is composed of two different |
| 622 |
< |
atom types, the average of the two resulting velocity autocorrelation |
| 623 |
< |
functions was used for comparisons. |
| 624 |
< |
|
| 625 |
< |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 626 |
< |
|
| 627 |
< |
The effects of the same subset of alternative electrostatic methods on |
| 628 |
< |
the {\it long-time} dynamics of charged systems were evaluated using |
| 629 |
< |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 630 |
< |
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 631 |
< |
autocorrelation function, \begin{equation} I(\omega) = |
| 632 |
< |
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 633 |
< |
\label{eq:powerSpec} |
| 634 |
< |
\end{equation} |
| 635 |
< |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 636 |
< |
NaCl crystal is composed of two different atom types, the average of |
| 637 |
< |
the two resulting power spectra was used for comparisons. Simulations |
| 638 |
< |
were performed under the microcanonical ensemble, and velocity |
| 639 |
< |
information was saved every 5~fs over 100~ps trajectories. |
| 640 |
< |
|
| 641 |
< |
\subsection{Representative Simulations}\label{sec:RepSims} |
| 642 |
< |
A variety of representative molecular simulations were analyzed to |
| 643 |
< |
determine the relative effectiveness of the pairwise summation |
| 644 |
< |
techniques in reproducing the energetics and dynamics exhibited by |
| 645 |
< |
{\sc spme}. We wanted to span the space of typical molecular |
| 646 |
< |
simulations (i.e. from liquids of neutral molecules to ionic |
| 647 |
< |
crystals), so the systems studied were: |
| 648 |
< |
|
| 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, |
| 653 |
< |
\item NaCl melts, |
| 654 |
< |
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 655 |
< |
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 656 |
< |
\item a 6\AA\ radius sphere of Argon in water. |
| 657 |
< |
\end{enumerate} |
| 658 |
< |
|
| 659 |
< |
By utilizing the pairwise techniques (outlined in section |
| 660 |
< |
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 661 |
< |
charged particles, and mixtures of the two, we hope to discern under |
| 662 |
< |
which conditions it will be possible to use one of the alternative |
| 663 |
< |
summation methodologies instead of the Ewald sum. |
| 664 |
< |
|
| 665 |
< |
For the solid and liquid water configurations, configurations were |
| 666 |
< |
taken at regular intervals from high temperature trajectories of 1000 |
| 667 |
< |
SPC/E water molecules. Each configuration was equilibrated |
| 668 |
< |
independently at a lower temperature (300K for the liquid, 200K for |
| 669 |
< |
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 670 |
< |
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 671 |
< |
these systems were selected and equilibrated in the same manner as the |
| 672 |
< |
water systems. In order to introduce measurable fluctuations in the |
| 673 |
< |
configuration energy differences, the crystalline simulations were |
| 674 |
< |
equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 675 |
< |
NaCl configurations needed to represent a fully disordered array of |
| 676 |
< |
point charges, so the high temperature of 7000K was selected for |
| 677 |
< |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 678 |
< |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 679 |
< |
water positions were then randomly swapped, and the resulting |
| 680 |
< |
configurations were again equilibrated individually. Finally, for the |
| 681 |
< |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 682 |
< |
waters within 6\AA\ of the center of the equilibrated water |
| 683 |
< |
configurations were converted to argon. |
| 684 |
< |
|
| 685 |
< |
These procedures guaranteed us a set of representative configurations |
| 686 |
< |
from chemically-relevant systems sampled from appropriate |
| 687 |
< |
ensembles. Force field parameters for the ions and Argon were taken |
| 688 |
< |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 689 |
< |
|
| 690 |
< |
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 691 |
< |
We compared the following alternative summation methods with results |
| 692 |
< |
from the reference method ({\sc spme}): |
| 693 |
< |
|
| 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, |
| 698 |
< |
and 0.3\AA$^{-1}$, |
| 699 |
< |
\item reaction field with an infinite dielectric constant, and |
| 700 |
< |
\item an unmodified cutoff. |
| 701 |
< |
\end{enumerate} |
| 702 |
< |
|
| 703 |
< |
Group-based cutoffs with a fifth-order polynomial switching function |
| 704 |
< |
were utilized for the reaction field simulations. Additionally, we |
| 705 |
< |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 706 |
< |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 707 |
< |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 708 |
< |
were performed using the {\sc oopse} molecular mechanics |
| 709 |
< |
package.\cite{Meineke05} All other portions of the energy calculation |
| 710 |
< |
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 711 |
< |
manner across all systems and configurations. |
| 712 |
< |
|
| 713 |
< |
The alternative methods were also evaluated with three different |
| 714 |
< |
cutoff radii (9, 12, and 15\AA). As noted previously, the |
| 715 |
< |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 716 |
< |
real-space and reciprocal-space portions of the Ewald calculation. |
| 717 |
< |
Typical molecular mechanics packages set this to a value dependent on |
| 718 |
< |
the cutoff radius and a tolerance (typically less than $1 \times |
| 719 |
< |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 720 |
< |
increasing accuracy at the expense of computational time spent on the |
| 721 |
< |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 722 |
< |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 723 |
< |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 724 |
< |
0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\ |
| 725 |
< |
respectively. |
| 726 |
< |
|
| 727 |
< |
\section{Configuration Energy Difference Results}\label{sec:EnergyResults} |
| 728 |
< |
In order to evaluate the performance of the pairwise electrostatic |
| 729 |
< |
summation methods for Monte Carlo (MC) simulations, the energy |
| 730 |
< |
differences between configurations were compared to the values |
| 731 |
< |
obtained when using {\sc spme}. The results for the combined |
| 732 |
< |
regression analysis of all of the systems are shown in figure |
| 733 |
< |
\ref{fig:delE}. |
| 734 |
< |
|
| 735 |
< |
\begin{figure} |
| 736 |
< |
\centering |
| 737 |
< |
\includegraphics[width=4.75in]{./figures/delEplot.pdf} |
| 738 |
< |
\caption{Statistical analysis of the quality of configurational energy |
| 739 |
< |
differences for a given electrostatic method compared with the |
| 740 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 741 |
< |
indicate $\Delta E$ values indistinguishable from those obtained using |
| 742 |
< |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 743 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 744 |
< |
inverted triangles).} |
| 745 |
< |
\label{fig:delE} |
| 746 |
< |
\end{figure} |
| 747 |
< |
|
| 748 |
< |
The most striking feature of this plot is how well the Shifted Force |
| 749 |
< |
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 750 |
< |
differences. For the undamped {\sc sf} method, and the |
| 751 |
< |
moderately-damped {\sc sp} methods, the results are nearly |
| 752 |
< |
indistinguishable from the Ewald results. The other common methods do |
| 753 |
< |
significantly less well. |
| 754 |
< |
|
| 755 |
< |
The unmodified cutoff method is essentially unusable. This is not |
| 756 |
< |
surprising since hard cutoffs give large energy fluctuations as atoms |
| 757 |
< |
or molecules move in and out of the cutoff |
| 758 |
< |
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
| 759 |
< |
some degree by using group based cutoffs with a switching |
| 760 |
< |
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 761 |
< |
significant improvement using the group-switched cutoff because the |
| 762 |
< |
salt and salt solution systems contain non-neutral groups. Section |
| 763 |
< |
\ref{sec:IndividualResults} includes results for systems comprised entirely |
| 764 |
< |
of neutral groups. |
| 765 |
< |
|
| 766 |
< |
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 767 |
< |
the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$ |
| 768 |
< |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 769 |
< |
results, particularly with a cutoff radius greater than 12 |
| 770 |
< |
\AA . Use of a larger damping parameter is more helpful for the |
| 771 |
< |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 772 |
< |
with larger cutoffs. |
| 773 |
< |
|
| 774 |
< |
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 775 |
< |
worse} correlation with Ewald. Overall, the undamped case is the best |
| 776 |
< |
performing set, as the correlation and quality of fits are |
| 777 |
< |
consistently superior regardless of the cutoff distance. The undamped |
| 778 |
< |
case is also less computationally demanding (because no evaluation of |
| 779 |
< |
the complementary error function is required). |
| 780 |
< |
|
| 781 |
< |
The reaction field results illustrates some of that method's |
| 782 |
< |
limitations, primarily that it was developed for use in homogeneous |
| 783 |
< |
systems; although it does provide results that are an improvement over |
| 784 |
< |
those from an unmodified cutoff. |
| 785 |
< |
|
| 786 |
< |
\section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults} |
| 787 |
< |
|
| 788 |
< |
Evaluation of pairwise methods for use in Molecular Dynamics |
| 789 |
< |
simulations requires consideration of effects on the forces and |
| 790 |
< |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
| 791 |
< |
regression results for the force and torque vector magnitudes, |
| 792 |
< |
respectively. The data in these figures was generated from an |
| 793 |
< |
accumulation of the statistics from all of the system types. |
| 794 |
< |
|
| 795 |
< |
\begin{figure} |
| 796 |
< |
\centering |
| 797 |
< |
\includegraphics[width=4.75in]{./figures/frcMagplot.pdf} |
| 798 |
< |
\caption{Statistical analysis of the quality of the force vector |
| 799 |
< |
magnitudes for a given electrostatic method compared with the |
| 800 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 801 |
< |
indicate force magnitude values indistinguishable from those obtained |
| 802 |
< |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 803 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 804 |
< |
inverted triangles).} |
| 805 |
< |
\label{fig:frcMag} |
| 806 |
< |
\end{figure} |
| 807 |
< |
|
| 808 |
< |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 809 |
< |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 810 |
< |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 811 |
< |
This is perhaps expected because this method explicitly incorporates a |
| 812 |
< |
smooth transition in the forces at the cutoff radius as well as the |
| 813 |
< |
neutralizing image charges. |
| 814 |
< |
|
| 815 |
< |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 816 |
< |
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 817 |
< |
poor, but using group based cutoffs and a switching function provides |
| 818 |
< |
an improvement much more significant than what was seen with $\Delta |
| 819 |
< |
E$. |
| 820 |
< |
|
| 821 |
< |
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 822 |
< |
method is generating usable forces. Further increases in damping, |
| 823 |
< |
while beneficial for simulations with a cutoff radius of 9\AA\ , is |
| 824 |
< |
detrimental to simulations with larger cutoff radii. |
| 825 |
< |
|
| 826 |
< |
The reaction field results are surprisingly good, considering the poor |
| 827 |
< |
quality of the fits for the $\Delta E$ results. There is still a |
| 828 |
< |
considerable degree of scatter in the data, but the forces correlate |
| 829 |
< |
well with the Ewald forces in general. We note that the reaction |
| 830 |
< |
field calculations do not include the pure NaCl systems, so these |
| 831 |
< |
results are partly biased towards conditions in which the method |
| 832 |
< |
performs more favorably. |
| 833 |
< |
|
| 834 |
< |
\begin{figure} |
| 835 |
< |
\centering |
| 836 |
< |
\includegraphics[width=4.75in]{./figures/trqMagplot.pdf} |
| 837 |
< |
\caption{Statistical analysis of the quality of the torque vector |
| 838 |
< |
magnitudes for a given electrostatic method compared with the |
| 839 |
< |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 840 |
< |
indicate torque magnitude values indistinguishable from those obtained |
| 841 |
< |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 842 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 843 |
< |
inverted triangles).} |
| 844 |
< |
\label{fig:trqMag} |
| 845 |
< |
\end{figure} |
| 846 |
< |
|
| 847 |
< |
Molecular torques were only available from the systems which contained |
| 848 |
< |
rigid molecules (i.e. the systems containing water). The data in |
| 849 |
< |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 850 |
< |
|
| 851 |
< |
Torques appear to be much more sensitive to charges at a longer |
| 852 |
< |
distance. The striking feature in comparing the new electrostatic |
| 853 |
< |
methods with {\sc spme} is how much the agreement improves with increasing |
| 854 |
< |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 855 |
< |
appears to be reproducing the {\sc spme} torques most accurately. |
| 856 |
< |
|
| 857 |
< |
Water molecules are dipolar, and the reaction field method reproduces |
| 858 |
< |
the effect of the surrounding polarized medium on each of the |
| 859 |
< |
molecular bodies. Therefore it is not surprising that reaction field |
| 860 |
< |
performs best of all of the methods on molecular torques. |
| 861 |
< |
|
| 862 |
< |
\section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults} |
| 863 |
< |
|
| 864 |
< |
It is clearly important that a new electrostatic method can reproduce |
| 865 |
< |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 866 |
< |
sum. However, the {\it directionality} of these vectors will also be |
| 867 |
< |
vital in calculating dynamical quantities accurately. Force and |
| 868 |
< |
torque directionalities were investigated by measuring the angles |
| 869 |
< |
formed between these vectors and the same vectors calculated using |
| 870 |
< |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 871 |
< |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 872 |
< |
distributions of the combined set over all system types. |
| 873 |
< |
|
| 874 |
< |
\begin{figure} |
| 875 |
< |
\centering |
| 876 |
< |
\includegraphics[width=4.75in]{./figures/frcTrqAngplot.pdf} |
| 877 |
< |
\caption{Statistical analysis of the width of the angular distribution |
| 878 |
< |
that the force and torque vectors from a given electrostatic method |
| 879 |
< |
make with their counterparts obtained using the reference Ewald sum. |
| 880 |
< |
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 881 |
< |
indicate force and torque directions indistinguishable from those |
| 882 |
< |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 883 |
< |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 884 |
< |
and 15\AA\ = inverted triangles).} |
| 885 |
< |
\label{fig:frcTrqAng} |
| 886 |
< |
\end{figure} |
| 887 |
< |
|
| 888 |
< |
Both the force and torque $\sigma^2$ results from the analysis of the |
| 889 |
< |
total accumulated system data are tabulated in figure |
| 890 |
< |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 891 |
< |
sp}) method would be essentially unusable for molecular dynamics |
| 892 |
< |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 893 |
< |
method, however, is generating force and torque vectors which are |
| 894 |
< |
within a few degrees of the Ewald results even with weak (or no) |
| 895 |
< |
damping. |
| 896 |
< |
|
| 897 |
< |
All of the sets (aside from the over-damped case) show the improvement |
| 898 |
< |
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 899 |
< |
from 9 to 12\AA\ typically results in a halving of the width of the |
| 900 |
< |
distribution, with a similar improvement when going from 12 to 15 |
| 901 |
< |
\AA . |
| 902 |
< |
|
| 903 |
< |
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 904 |
< |
all do equivalently well at capturing the direction of both the force |
| 905 |
< |
and torque vectors. Using the electrostatic damping improves the |
| 906 |
< |
angular behavior significantly for the {\sc sp} and moderately for the |
| 907 |
< |
{\sc sf} methods. Over-damping is detrimental to both methods. Again |
| 908 |
< |
it is important to recognize that the force vectors cover all |
| 909 |
< |
particles in all seven systems, while torque vectors are only |
| 910 |
< |
available for neutral molecular groups. Damping is more beneficial to |
| 911 |
< |
charged bodies, and this observation is investigated further in |
| 912 |
< |
section \ref{sec:IndividualResults}. |
| 913 |
< |
|
| 914 |
< |
Although not discussed previously, group based cutoffs can be applied |
| 915 |
< |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 916 |
< |
will reintroduce small discontinuities at the cutoff radius, but the |
| 917 |
< |
effects of these can be minimized by utilizing a switching function. |
| 918 |
< |
Though there are no significant benefits or drawbacks observed in |
| 919 |
< |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 920 |
< |
is a measurable improvement in the directionality of the forces and |
| 921 |
< |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 922 |
< |
obtained both without (N) and with (Y) group based cutoffs and a |
| 923 |
< |
switching function. Note that the $\alpha$ values have units of |
| 924 |
< |
\AA$^{-1}$ and the variance values have units of degrees$^2$. The |
| 925 |
< |
{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much |
| 926 |
< |
narrower angular distributions when using group-based cutoffs. The |
| 927 |
< |
{\sc sf} method likewise shows improvement in the undamped and lightly |
| 928 |
< |
damped cases. |
| 929 |
< |
|
| 930 |
< |
\begin{table} |
| 931 |
< |
\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$) |
| 932 |
< |
THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A |
| 933 |
< |
GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING |
| 934 |
< |
THE REFERENCE EWALD SUMMATION} |
| 935 |
< |
|
| 936 |
< |
\footnotesize |
| 937 |
< |
\begin{center} |
| 938 |
< |
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 939 |
< |
\toprule |
| 940 |
< |
\toprule |
| 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 & |
| 944 |
< |
$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & |
| 945 |
< |
$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\ |
| 946 |
< |
|
| 947 |
< |
\midrule |
| 948 |
< |
|
| 949 |
< |
9\AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
| 950 |
< |
& \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} \\ |
| 951 |
< |
12\AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
| 952 |
< |
& \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} \\ |
| 953 |
< |
15\AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
| 954 |
< |
& \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} \\ |
| 955 |
< |
|
| 956 |
< |
\midrule |
| 957 |
< |
|
| 958 |
< |
9\AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
| 959 |
< |
& \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} \\ |
| 960 |
< |
12\AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
| 961 |
< |
& \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} \\ |
| 962 |
< |
15\AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
| 963 |
< |
& \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} \\ |
| 964 |
< |
|
| 965 |
< |
\bottomrule |
| 966 |
< |
\end{tabular} |
| 967 |
< |
\end{center} |
| 968 |
< |
\label{tab:groupAngle} |
| 969 |
< |
\end{table} |
| 970 |
< |
|
| 971 |
< |
One additional trend in table \ref{tab:groupAngle} is that the |
| 972 |
< |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 973 |
< |
increases, something that is more obvious with group-based cutoffs. |
| 974 |
< |
The complimentary error function inserted into the potential weakens |
| 975 |
< |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 976 |
< |
However, at larger values of $\alpha$, it is possible to over-damp the |
| 977 |
< |
electrostatic interaction and to remove it completely. Kast |
| 978 |
< |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 979 |
< |
values for these types of electrostatic summation methods by fitting |
| 980 |
< |
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 981 |
< |
0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\ |
| 982 |
< |
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 983 |
< |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 984 |
< |
these findings, choices this high would introduce error in the |
| 985 |
< |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 986 |
< |
observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial, |
| 987 |
< |
but damping may be unnecessary when using the {\sc sf} method. |
| 988 |
< |
|
| 989 |
< |
\section{Individual System Analysis Results}\label{sec:IndividualResults} |
| 990 |
< |
|
| 991 |
< |
The combined results of the previous sections show how the pairwise |
| 992 |
< |
methods compare to the Ewald summation in the general sense over all |
| 993 |
< |
of the system types. It is also useful to consider each of the |
| 994 |
< |
studied systems in an individual fashion, so that we can identify |
| 995 |
< |
conditions that are particularly difficult for a selected pairwise |
| 996 |
< |
method to address. This allows us to further establish the limitations |
| 997 |
< |
of these pairwise techniques. Below, the energy difference, force |
| 998 |
< |
vector, and torque vector analyses are presented on an individual |
| 999 |
< |
system basis. |
| 1000 |
< |
|
| 1001 |
< |
\subsection{SPC/E Water Results}\label{sec:WaterResults} |
| 1002 |
< |
|
| 1003 |
< |
The first system considered was liquid water at 300K using the SPC/E |
| 1004 |
< |
model of water.\cite{Berendsen87} The results for the energy gap |
| 1005 |
< |
comparisons and the force and torque vector magnitude comparisons are |
| 1006 |
< |
shown in table \ref{tab:spce}. The force and torque vector |
| 1007 |
< |
directionality results are displayed separately in table |
| 1008 |
< |
\ref{tab:spceAng}, where the effect of group-based cutoffs and |
| 1009 |
< |
switching functions on the {\sc sp} and {\sc sf} potentials are also |
| 1010 |
< |
investigated. In all of the individual results table, the method |
| 1011 |
< |
abbreviations are as follows: |
| 1012 |
< |
|
| 1013 |
< |
\begin{itemize}[itemsep=0pt] |
| 1014 |
< |
\item PC = Pure Cutoff, |
| 1015 |
< |
\item SP = Shifted Potential, |
| 1016 |
< |
\item SF = Shifted Force, |
| 1017 |
< |
\item GSC = Group Switched Cutoff, |
| 1018 |
< |
\item RF = Reaction Field (where $\varepsilon \approx\infty$), |
| 1019 |
< |
\item GSSP = Group Switched Shifted Potential, and |
| 1020 |
< |
\item GSSF = Group Switched Shifted Force. |
| 1021 |
< |
\end{itemize} |
| 1022 |
< |
|
| 1023 |
< |
\begin{table}[htbp] |
| 1024 |
< |
\centering |
| 1025 |
< |
\caption{REGRESSION RESULTS OF THE LIQUID WATER SYSTEM FOR THE |
| 1026 |
< |
$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it middle}) |
| 1027 |
< |
AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1028 |
< |
|
| 1029 |
< |
\footnotesize |
| 1030 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1031 |
< |
\toprule |
| 1032 |
< |
\toprule |
| 1033 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1034 |
< |
\cmidrule(lr){3-4} |
| 1035 |
< |
\cmidrule(lr){5-6} |
| 1036 |
< |
\cmidrule(l){7-8} |
| 1037 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1038 |
< |
\midrule |
| 1039 |
< |
PC & & 3.046 & 0.002 & -3.018 & 0.002 & 4.719 & 0.005 \\ |
| 1040 |
< |
SP & 0.0 & 1.035 & 0.218 & 0.908 & 0.313 & 1.037 & 0.470 \\ |
| 1041 |
< |
& 0.1 & 1.021 & 0.387 & 0.965 & 0.752 & 1.006 & 0.947 \\ |
| 1042 |
< |
& 0.2 & 0.997 & 0.962 & 1.001 & 0.994 & 0.994 & 0.996 \\ |
| 1043 |
< |
& 0.3 & 0.984 & 0.980 & 0.997 & 0.985 & 0.982 & 0.987 \\ |
| 1044 |
< |
SF & 0.0 & 0.977 & 0.974 & 0.996 & 0.992 & 0.991 & 0.997 \\ |
| 1045 |
< |
& 0.1 & 0.983 & 0.974 & 1.001 & 0.994 & 0.996 & 0.998 \\ |
| 1046 |
< |
& 0.2 & 0.992 & 0.989 & 1.001 & 0.995 & 0.994 & 0.996 \\ |
| 1047 |
< |
& 0.3 & 0.984 & 0.980 & 0.996 & 0.985 & 0.982 & 0.987 \\ |
| 1048 |
< |
GSC & & 0.918 & 0.862 & 0.852 & 0.756 & 0.801 & 0.700 \\ |
| 1049 |
< |
RF & & 0.971 & 0.958 & 0.975 & 0.987 & 0.959 & 0.983 \\ |
| 1050 |
< |
\midrule |
| 1051 |
< |
PC & & -1.647 & 0.000 & -0.127 & 0.000 & -0.979 & 0.000 \\ |
| 1052 |
< |
SP & 0.0 & 0.735 & 0.368 & 0.813 & 0.537 & 0.865 & 0.659 \\ |
| 1053 |
< |
& 0.1 & 0.850 & 0.612 & 0.956 & 0.887 & 0.992 & 0.979 \\ |
| 1054 |
< |
& 0.2 & 0.996 & 0.989 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1055 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ |
| 1056 |
< |
SF & 0.0 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 0.999 \\ |
| 1057 |
< |
& 0.1 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1058 |
< |
& 0.2 & 0.999 & 0.998 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1059 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ |
| 1060 |
< |
GSC & & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1061 |
< |
RF & & 0.999 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1062 |
< |
\midrule |
| 1063 |
< |
PC & & 2.387 & 0.000 & 0.183 & 0.000 & 1.282 & 0.000 \\ |
| 1064 |
< |
SP & 0.0 & 0.847 & 0.543 & 0.904 & 0.694 & 0.935 & 0.786 \\ |
| 1065 |
< |
& 0.1 & 0.922 & 0.749 & 0.980 & 0.934 & 0.996 & 0.988 \\ |
| 1066 |
< |
& 0.2 & 0.987 & 0.985 & 0.989 & 0.992 & 0.990 & 0.993 \\ |
| 1067 |
< |
& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ |
| 1068 |
< |
SF & 0.0 & 0.978 & 0.990 & 0.988 & 0.997 & 0.993 & 0.999 \\ |
| 1069 |
< |
& 0.1 & 0.983 & 0.991 & 0.993 & 0.997 & 0.997 & 0.999 \\ |
| 1070 |
< |
& 0.2 & 0.986 & 0.989 & 0.989 & 0.992 & 0.990 & 0.993 \\ |
| 1071 |
< |
& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ |
| 1072 |
< |
GSC & & 0.995 & 0.981 & 0.999 & 0.991 & 1.001 & 0.994 \\ |
| 1073 |
< |
RF & & 0.993 & 0.989 & 0.998 & 0.996 & 1.000 & 0.999 \\ |
| 1074 |
< |
\bottomrule |
| 1075 |
< |
\end{tabular} |
| 1076 |
< |
\label{tab:spce} |
| 1077 |
< |
\end{table} |
| 1078 |
< |
|
| 1079 |
< |
\begin{table}[htbp] |
| 1080 |
< |
\centering |
| 1081 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1082 |
< |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE LIQUID WATER |
| 1083 |
< |
SYSTEM} |
| 1084 |
< |
|
| 1085 |
< |
\footnotesize |
| 1086 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1087 |
< |
\toprule |
| 1088 |
< |
\toprule |
| 1089 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1090 |
< |
\cmidrule(lr){3-5} |
| 1091 |
< |
\cmidrule(l){6-8} |
| 1092 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1093 |
< |
\midrule |
| 1094 |
< |
PC & & 783.759 & 481.353 & 332.677 & 248.674 & 144.382 & 98.535 \\ |
| 1095 |
< |
SP & 0.0 & 659.440 & 380.699 & 250.002 & 235.151 & 134.661 & 88.135 \\ |
| 1096 |
< |
& 0.1 & 293.849 & 67.772 & 11.609 & 105.090 & 23.813 & 4.369 \\ |
| 1097 |
< |
& 0.2 & 5.975 & 0.136 & 0.094 & 5.553 & 1.784 & 1.536 \\ |
| 1098 |
< |
& 0.3 & 0.725 & 0.707 & 0.693 & 7.293 & 6.933 & 6.748 \\ |
| 1099 |
< |
SF & 0.0 & 2.238 & 0.713 & 0.292 & 3.290 & 1.090 & 0.416 \\ |
| 1100 |
< |
& 0.1 & 2.238 & 0.524 & 0.115 & 3.184 & 0.945 & 0.326 \\ |
| 1101 |
< |
& 0.2 & 0.374 & 0.102 & 0.094 & 2.598 & 1.755 & 1.537 \\ |
| 1102 |
< |
& 0.3 & 0.721 & 0.707 & 0.693 & 7.322 & 6.933 & 6.748 \\ |
| 1103 |
< |
GSC & & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ |
| 1104 |
< |
RF & & 2.091 & 0.403 & 0.113 & 3.583 & 1.071 & 0.399 \\ |
| 1105 |
< |
\midrule |
| 1106 |
< |
GSSP & 0.0 & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ |
| 1107 |
< |
& 0.1 & 1.879 & 0.291 & 0.057 & 3.983 & 1.117 & 0.370 \\ |
| 1108 |
< |
& 0.2 & 0.443 & 0.103 & 0.093 & 2.821 & 1.794 & 1.532 \\ |
| 1109 |
< |
& 0.3 & 0.728 & 0.694 & 0.692 & 7.387 & 6.942 & 6.748 \\ |
| 1110 |
< |
GSSF & 0.0 & 1.298 & 0.270 & 0.083 & 3.098 & 0.992 & 0.375 \\ |
| 1111 |
< |
& 0.1 & 1.296 & 0.210 & 0.044 & 3.055 & 0.922 & 0.330 \\ |
| 1112 |
< |
& 0.2 & 0.433 & 0.104 & 0.093 & 2.895 & 1.797 & 1.532 \\ |
| 1113 |
< |
& 0.3 & 0.728 & 0.694 & 0.692 & 7.410 & 6.942 & 6.748 \\ |
| 1114 |
< |
\bottomrule |
| 1115 |
< |
\end{tabular} |
| 1116 |
< |
\label{tab:spceAng} |
| 1117 |
< |
\end{table} |
| 1118 |
< |
|
| 1119 |
< |
The water results parallel the combined results seen in sections |
| 1120 |
< |
\ref{sec:EnergyResults} through \ref{sec:FTDirResults}. There is good |
| 1121 |
< |
agreement with {\sc spme} in both energetic and dynamic behavior when |
| 1122 |
< |
using the {\sc sf} method with and without damping. The {\sc sp} |
| 1123 |
< |
method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly |
| 1124 |
< |
with cutoff radii greater than 12\AA. Over-damping the electrostatics |
| 1125 |
< |
reduces the agreement between both these methods and {\sc spme}. |
| 1126 |
< |
|
| 1127 |
< |
The pure cutoff ({\sc pc}) method performs poorly, again mirroring the |
| 1128 |
< |
observations from the combined results. In contrast to these results, however, the use of a switching function and group |
| 1129 |
< |
based cutoffs greatly improves the results for these neutral water |
| 1130 |
< |
molecules. The group switched cutoff ({\sc gsc}) does not mimic the |
| 1131 |
< |
energetics of {\sc spme} as well as the {\sc sp} (with moderate |
| 1132 |
< |
damping) and {\sc sf} methods, but the dynamics are quite good. The |
| 1133 |
< |
switching functions correct discontinuities in the potential and |
| 1134 |
< |
forces, leading to these improved results. Such improvements with the |
| 1135 |
< |
use of a switching function have been recognized in previous |
| 1136 |
< |
studies,\cite{Andrea83,Steinbach94} and this proves to be a useful |
| 1137 |
< |
tactic for stably incorporating local area electrostatic effects. |
| 1138 |
< |
|
| 1139 |
< |
The reaction field ({\sc rf}) method simply extends upon the results |
| 1140 |
< |
observed in the {\sc gsc} case. Both methods are similar in form |
| 1141 |
< |
(i.e. neutral groups, switching function), but {\sc rf} incorporates |
| 1142 |
< |
an added effect from the external dielectric. This similarity |
| 1143 |
< |
translates into the same good dynamic results and improved energetic |
| 1144 |
< |
agreement with {\sc spme}. Though this agreement is not to the level |
| 1145 |
< |
of the moderately damped {\sc sp} and {\sc sf} methods, these results |
| 1146 |
< |
show how incorporating some implicit properties of the surroundings |
| 1147 |
< |
(i.e. $\epsilon_\textrm{S}$) can improve the solvent depiction. |
| 1148 |
< |
|
| 1149 |
< |
As a final note for the liquid water system, use of group cutoffs and a |
| 1150 |
< |
switching function leads to noticeable improvements in the {\sc sp} |
| 1151 |
< |
and {\sc sf} methods, primarily in directionality of the force and |
| 1152 |
< |
torque vectors (table \ref{tab:spceAng}). The {\sc sp} method shows |
| 1153 |
< |
significant narrowing of the angle distribution when using little to |
| 1154 |
< |
no damping and only modest improvement for the recommended conditions |
| 1155 |
< |
($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA). The |
| 1156 |
< |
{\sc sf} method shows modest narrowing across all damping and cutoff |
| 1157 |
< |
ranges of interest. When over-damping these methods, group cutoffs and |
| 1158 |
< |
the switching function do not improve the force and torque |
| 1159 |
< |
directionalities. |
| 1160 |
< |
|
| 1161 |
< |
\subsection{SPC/E Ice I$_\textrm{c}$ Results}\label{sec:IceResults} |
| 1162 |
< |
|
| 1163 |
< |
In addition to the disordered molecular system above, the ordered |
| 1164 |
< |
molecular system of ice I$_\textrm{c}$ was also considered. Ice |
| 1165 |
< |
polymorph could have been used to fit this role; however, ice |
| 1166 |
< |
I$_\textrm{c}$ was chosen because it can form an ideal periodic |
| 1167 |
< |
lattice with the same number of water molecules used in the disordered |
| 1168 |
< |
liquid state case. The results for the energy gap comparisons and the |
| 1169 |
< |
force and torque vector magnitude comparisons are shown in table |
| 1170 |
< |
\ref{tab:ice}. The force and torque vector directionality results are |
| 1171 |
< |
displayed separately in table \ref{tab:iceAng}, where the effect of |
| 1172 |
< |
group-based cutoffs and switching functions on the {\sc sp} and {\sc |
| 1173 |
< |
sf} potentials are also displayed. |
| 1174 |
< |
|
| 1175 |
< |
\begin{table}[htbp] |
| 1176 |
< |
\centering |
| 1177 |
< |
\caption{REGRESSION RESULTS OF THE ICE I$_\textrm{c}$ SYSTEM FOR |
| 1178 |
< |
$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it |
| 1179 |
< |
middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1180 |
< |
|
| 1181 |
< |
\footnotesize |
| 1182 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1183 |
< |
\toprule |
| 1184 |
< |
\toprule |
| 1185 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1186 |
< |
\cmidrule(lr){3-4} |
| 1187 |
< |
\cmidrule(lr){5-6} |
| 1188 |
< |
\cmidrule(l){7-8} |
| 1189 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1190 |
< |
\midrule |
| 1191 |
< |
PC & & 19.897 & 0.047 & -29.214 & 0.048 & -3.771 & 0.001 \\ |
| 1192 |
< |
SP & 0.0 & -0.014 & 0.000 & 2.135 & 0.347 & 0.457 & 0.045 \\ |
| 1193 |
< |
& 0.1 & 0.321 & 0.017 & 1.490 & 0.584 & 0.886 & 0.796 \\ |
| 1194 |
< |
& 0.2 & 0.896 & 0.872 & 1.011 & 0.998 & 0.997 & 0.999 \\ |
| 1195 |
< |
& 0.3 & 0.983 & 0.997 & 0.992 & 0.997 & 0.991 & 0.997 \\ |
| 1196 |
< |
SF & 0.0 & 0.943 & 0.979 & 1.048 & 0.978 & 0.995 & 0.999 \\ |
| 1197 |
< |
& 0.1 & 0.948 & 0.979 & 1.044 & 0.983 & 1.000 & 0.999 \\ |
| 1198 |
< |
& 0.2 & 0.982 & 0.997 & 0.969 & 0.960 & 0.997 & 0.999 \\ |
| 1199 |
< |
& 0.3 & 0.985 & 0.997 & 0.961 & 0.961 & 0.991 & 0.997 \\ |
| 1200 |
< |
GSC & & 0.983 & 0.985 & 0.966 & 0.994 & 1.003 & 0.999 \\ |
| 1201 |
< |
RF & & 0.924 & 0.944 & 0.990 & 0.996 & 0.991 & 0.998 \\ |
| 1202 |
< |
\midrule |
| 1203 |
< |
PC & & -4.375 & 0.000 & 6.781 & 0.000 & -3.369 & 0.000 \\ |
| 1204 |
< |
SP & 0.0 & 0.515 & 0.164 & 0.856 & 0.426 & 0.743 & 0.478 \\ |
| 1205 |
< |
& 0.1 & 0.696 & 0.405 & 0.977 & 0.817 & 0.974 & 0.964 \\ |
| 1206 |
< |
& 0.2 & 0.981 & 0.980 & 1.001 & 1.000 & 1.000 & 1.000 \\ |
| 1207 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.999 & 0.997 & 0.999 \\ |
| 1208 |
< |
SF & 0.0 & 0.991 & 0.995 & 1.003 & 0.998 & 0.999 & 1.000 \\ |
| 1209 |
< |
& 0.1 & 0.992 & 0.995 & 1.003 & 0.998 & 1.000 & 1.000 \\ |
| 1210 |
< |
& 0.2 & 0.998 & 0.998 & 0.981 & 0.962 & 1.000 & 1.000 \\ |
| 1211 |
< |
& 0.3 & 0.996 & 0.998 & 0.976 & 0.957 & 0.997 & 0.999 \\ |
| 1212 |
< |
GSC & & 0.997 & 0.996 & 0.998 & 0.999 & 1.000 & 1.000 \\ |
| 1213 |
< |
RF & & 0.988 & 0.989 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1214 |
< |
\midrule |
| 1215 |
< |
PC & & -6.367 & 0.000 & -3.552 & 0.000 & -3.447 & 0.000 \\ |
| 1216 |
< |
SP & 0.0 & 0.643 & 0.409 & 0.833 & 0.607 & 0.961 & 0.805 \\ |
| 1217 |
< |
& 0.1 & 0.791 & 0.683 & 0.957 & 0.914 & 1.000 & 0.989 \\ |
| 1218 |
< |
& 0.2 & 0.974 & 0.991 & 0.993 & 0.998 & 0.993 & 0.998 \\ |
| 1219 |
< |
& 0.3 & 0.976 & 0.992 & 0.977 & 0.992 & 0.977 & 0.992 \\ |
| 1220 |
< |
SF & 0.0 & 0.979 & 0.997 & 0.992 & 0.999 & 0.994 & 1.000 \\ |
| 1221 |
< |
& 0.1 & 0.984 & 0.997 & 0.996 & 0.999 & 0.998 & 1.000 \\ |
| 1222 |
< |
& 0.2 & 0.991 & 0.997 & 0.974 & 0.958 & 0.993 & 0.998 \\ |
| 1223 |
< |
& 0.3 & 0.977 & 0.992 & 0.956 & 0.948 & 0.977 & 0.992 \\ |
| 1224 |
< |
GSC & & 0.999 & 0.997 & 0.996 & 0.999 & 1.002 & 1.000 \\ |
| 1225 |
< |
RF & & 0.994 & 0.997 & 0.997 & 0.999 & 1.000 & 1.000 \\ |
| 1226 |
< |
\bottomrule |
| 1227 |
< |
\end{tabular} |
| 1228 |
< |
\label{tab:ice} |
| 1229 |
< |
\end{table} |
| 1230 |
< |
|
| 1231 |
< |
\begin{table}[htbp] |
| 1232 |
< |
\centering |
| 1233 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1234 |
< |
OF THE FORCE AND TORQUE VECTORS IN THE ICE I$_\textrm{c}$ SYSTEM} |
| 1235 |
< |
|
| 1236 |
< |
\footnotesize |
| 1237 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1238 |
< |
\toprule |
| 1239 |
< |
\toprule |
| 1240 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque |
| 1241 |
< |
$\sigma^2$} \\ |
| 1242 |
< |
\cmidrule(lr){3-5} |
| 1243 |
< |
\cmidrule(l){6-8} |
| 1244 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1245 |
< |
\midrule |
| 1246 |
< |
PC & & 2128.921 & 603.197 & 715.579 & 329.056 & 221.397 & 81.042 \\ |
| 1247 |
< |
SP & 0.0 & 1429.341 & 470.320 & 447.557 & 301.678 & 197.437 & 73.840 \\ |
| 1248 |
< |
& 0.1 & 590.008 & 107.510 & 18.883 & 118.201 & 32.472 & 3.599 \\ |
| 1249 |
< |
& 0.2 & 10.057 & 0.105 & 0.038 & 2.875 & 0.572 & 0.518 \\ |
| 1250 |
< |
& 0.3 & 0.245 & 0.260 & 0.262 & 2.365 & 2.396 & 2.327 \\ |
| 1251 |
< |
SF & 0.0 & 1.745 & 1.161 & 0.212 & 1.135 & 0.426 & 0.155 \\ |
| 1252 |
< |
& 0.1 & 1.721 & 0.868 & 0.082 & 1.118 & 0.358 & 0.118 \\ |
| 1253 |
< |
& 0.2 & 0.201 & 0.040 & 0.038 & 0.786 & 0.555 & 0.518 \\ |
| 1254 |
< |
& 0.3 & 0.241 & 0.260 & 0.262 & 2.368 & 2.400 & 2.327 \\ |
| 1255 |
< |
GSC & & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ |
| 1256 |
< |
RF & & 2.887 & 0.217 & 0.107 & 1.006 & 0.281 & 0.085 \\ |
| 1257 |
< |
\midrule |
| 1258 |
< |
GSSP & 0.0 & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ |
| 1259 |
< |
& 0.1 & 1.341 & 0.123 & 0.037 & 0.835 & 0.234 & 0.085 \\ |
| 1260 |
< |
& 0.2 & 0.558 & 0.040 & 0.037 & 0.823 & 0.557 & 0.519 \\ |
| 1261 |
< |
& 0.3 & 0.250 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ |
| 1262 |
< |
GSSF & 0.0 & 2.124 & 0.132 & 0.069 & 0.919 & 0.263 & 0.099 \\ |
| 1263 |
< |
& 0.1 & 2.165 & 0.101 & 0.035 & 0.895 & 0.244 & 0.096 \\ |
| 1264 |
< |
& 0.2 & 0.706 & 0.040 & 0.037 & 0.870 & 0.559 & 0.519 \\ |
| 1265 |
< |
& 0.3 & 0.251 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ |
| 1266 |
< |
\bottomrule |
| 1267 |
< |
\end{tabular} |
| 1268 |
< |
\label{tab:iceAng} |
| 1269 |
< |
\end{table} |
| 1270 |
< |
|
| 1271 |
< |
Highly ordered systems are a difficult test for the pairwise methods |
| 1272 |
< |
in that they lack the implicit periodicity of the Ewald summation. As |
| 1273 |
< |
expected, the energy gap agreement with {\sc spme} is reduced for the |
| 1274 |
< |
{\sc sp} and {\sc sf} methods with parameters that were ideal for the |
| 1275 |
< |
disordered liquid system. Moving to higher $R_\textrm{c}$ helps |
| 1276 |
< |
improve the agreement, though at an increase in computational cost. |
| 1277 |
< |
The dynamics of this crystalline system (both in magnitude and |
| 1278 |
< |
direction) are little affected. Both methods still reproduce the Ewald |
| 1279 |
< |
behavior with the same parameter recommendations from the previous |
| 1280 |
< |
section. |
| 1281 |
< |
|
| 1282 |
< |
It is also worth noting that {\sc rf} exhibits improved energy gap |
| 1283 |
< |
results over the liquid water system. One possible explanation is |
| 1284 |
< |
that the ice I$_\textrm{c}$ crystal is ordered such that the net |
| 1285 |
< |
dipole moment of the crystal is zero. With $\epsilon_\textrm{S} = |
| 1286 |
< |
\infty$, the reaction field incorporates this structural organization |
| 1287 |
< |
by actively enforcing a zeroed dipole moment within each cutoff |
| 1288 |
< |
sphere. |
| 1289 |
< |
|
| 1290 |
< |
\subsection{NaCl Melt Results}\label{sec:SaltMeltResults} |
| 1291 |
< |
|
| 1292 |
< |
A high temperature NaCl melt was tested to gauge the accuracy of the |
| 1293 |
< |
pairwise summation methods in a disordered system of charges. The |
| 1294 |
< |
results for the energy gap comparisons and the force vector magnitude |
| 1295 |
< |
comparisons are shown in table \ref{tab:melt}. The force vector |
| 1296 |
< |
directionality results are displayed separately in table |
| 1297 |
< |
\ref{tab:meltAng}. |
| 1298 |
< |
|
| 1299 |
< |
\begin{table}[htbp] |
| 1300 |
< |
\centering |
| 1301 |
< |
\caption{REGRESSION RESULTS OF THE MOLTEN SODIUM CHLORIDE SYSTEM FOR |
| 1302 |
< |
$\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES ({\it |
| 1303 |
< |
lower})} |
| 1304 |
< |
|
| 1305 |
< |
\footnotesize |
| 1306 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1307 |
< |
\toprule |
| 1308 |
< |
\toprule |
| 1309 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1310 |
< |
\cmidrule(lr){3-4} |
| 1311 |
< |
\cmidrule(lr){5-6} |
| 1312 |
< |
\cmidrule(l){7-8} |
| 1313 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1314 |
< |
\midrule |
| 1315 |
< |
PC & & -0.008 & 0.000 & -0.049 & 0.005 & -0.136 & 0.020 \\ |
| 1316 |
< |
SP & 0.0 & 0.928 & 0.996 & 0.931 & 0.998 & 0.950 & 0.999 \\ |
| 1317 |
< |
& 0.1 & 0.977 & 0.998 & 0.998 & 1.000 & 0.997 & 1.000 \\ |
| 1318 |
< |
& 0.2 & 0.960 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ |
| 1319 |
< |
& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ |
| 1320 |
< |
SF & 0.0 & 0.996 & 1.000 & 0.995 & 1.000 & 0.997 & 1.000 \\ |
| 1321 |
< |
& 0.1 & 1.021 & 1.000 & 1.024 & 1.000 & 1.007 & 1.000 \\ |
| 1322 |
< |
& 0.2 & 0.966 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ |
| 1323 |
< |
& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ |
| 1324 |
< |
\midrule |
| 1325 |
< |
PC & & 1.103 & 0.000 & 0.989 & 0.000 & 0.802 & 0.000 \\ |
| 1326 |
< |
SP & 0.0 & 0.973 & 0.981 & 0.975 & 0.988 & 0.979 & 0.992 \\ |
| 1327 |
< |
& 0.1 & 0.987 & 0.992 & 0.993 & 0.998 & 0.997 & 0.999 \\ |
| 1328 |
< |
& 0.2 & 0.993 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ |
| 1329 |
< |
& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ |
| 1330 |
< |
SF & 0.0 & 0.996 & 0.997 & 0.997 & 0.999 & 0.998 & 1.000 \\ |
| 1331 |
< |
& 0.1 & 1.000 & 0.997 & 1.001 & 0.999 & 1.000 & 1.000 \\ |
| 1332 |
< |
& 0.2 & 0.994 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ |
| 1333 |
< |
& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ |
| 1334 |
< |
\bottomrule |
| 1335 |
< |
\end{tabular} |
| 1336 |
< |
\label{tab:melt} |
| 1337 |
< |
\end{table} |
| 1338 |
< |
|
| 1339 |
< |
\begin{table}[htbp] |
| 1340 |
< |
\centering |
| 1341 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1342 |
< |
OF THE FORCE VECTORS IN THE MOLTEN SODIUM CHLORIDE SYSTEM} |
| 1343 |
< |
|
| 1344 |
< |
\footnotesize |
| 1345 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1346 |
< |
\toprule |
| 1347 |
< |
\toprule |
| 1348 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1349 |
< |
\cmidrule(lr){3-5} |
| 1350 |
< |
\cmidrule(l){6-8} |
| 1351 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ |
| 1352 |
< |
\midrule |
| 1353 |
< |
PC & & 13.294 & 8.035 & 5.366 \\ |
| 1354 |
< |
SP & 0.0 & 13.316 & 8.037 & 5.385 \\ |
| 1355 |
< |
& 0.1 & 5.705 & 1.391 & 0.360 \\ |
| 1356 |
< |
& 0.2 & 2.415 & 7.534 & 13.927 \\ |
| 1357 |
< |
& 0.3 & 23.769 & 67.306 & 57.252 \\ |
| 1358 |
< |
SF & 0.0 & 1.693 & 0.603 & 0.256 \\ |
| 1359 |
< |
& 0.1 & 1.687 & 0.653 & 0.272 \\ |
| 1360 |
< |
& 0.2 & 2.598 & 7.523 & 13.930 \\ |
| 1361 |
< |
& 0.3 & 23.734 & 67.305 & 57.252 \\ |
| 1362 |
< |
\bottomrule |
| 1363 |
< |
\end{tabular} |
| 1364 |
< |
\label{tab:meltAng} |
| 1365 |
< |
\end{table} |
| 1366 |
< |
|
| 1367 |
< |
The molten NaCl system shows more sensitivity to the electrostatic |
| 1368 |
< |
damping than the water systems. The most noticeable point is that the |
| 1369 |
< |
undamped {\sc sf} method does very well at replicating the {\sc spme} |
| 1370 |
< |
configurational energy differences and forces. Light damping appears |
| 1371 |
< |
to minimally improve the dynamics, but this comes with a deterioration |
| 1372 |
< |
of the energy gap results. In contrast, this light damping improves |
| 1373 |
< |
the {\sc sp} energy gaps and forces. Moderate and heavy electrostatic |
| 1374 |
< |
damping reduce the agreement with {\sc spme} for both methods. From |
| 1375 |
< |
these observations, the undamped {\sc sf} method is the best choice |
| 1376 |
< |
for disordered systems of charges. |
| 1377 |
< |
|
| 1378 |
< |
\subsection{NaCl Crystal Results}\label{sec:SaltCrystalResults} |
| 1379 |
< |
|
| 1380 |
< |
Similar to the use of ice I$_\textrm{c}$ to investigate the role of |
| 1381 |
< |
order in molecular systems on the effectiveness of the pairwise |
| 1382 |
< |
methods, the 1000K NaCl crystal system was used to investigate the |
| 1383 |
< |
accuracy of the pairwise summation methods in an ordered system of |
| 1384 |
< |
charged particles. The results for the energy gap comparisons and the |
| 1385 |
< |
force vector magnitude comparisons are shown in table \ref{tab:salt}. |
| 1386 |
< |
The force vector directionality results are displayed separately in |
| 1387 |
< |
table \ref{tab:saltAng}. |
| 1388 |
< |
|
| 1389 |
< |
\begin{table}[htbp] |
| 1390 |
< |
\centering |
| 1391 |
< |
\caption{REGRESSION RESULTS OF THE CRYSTALLINE SODIUM CHLORIDE |
| 1392 |
< |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES |
| 1393 |
< |
({\it lower})} |
| 1394 |
< |
|
| 1395 |
< |
\footnotesize |
| 1396 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1397 |
< |
\toprule |
| 1398 |
< |
\toprule |
| 1399 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1400 |
< |
\cmidrule(lr){3-4} |
| 1401 |
< |
\cmidrule(lr){5-6} |
| 1402 |
< |
\cmidrule(l){7-8} |
| 1403 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1404 |
< |
\midrule |
| 1405 |
< |
PC & & -20.241 & 0.228 & -20.248 & 0.229 & -20.239 & 0.228 \\ |
| 1406 |
< |
SP & 0.0 & 1.039 & 0.733 & 2.037 & 0.565 & 1.225 & 0.743 \\ |
| 1407 |
< |
& 0.1 & 1.049 & 0.865 & 1.424 & 0.784 & 1.029 & 0.980 \\ |
| 1408 |
< |
& 0.2 & 0.982 & 0.976 & 0.969 & 0.980 & 0.960 & 0.980 \\ |
| 1409 |
< |
& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.945 \\ |
| 1410 |
< |
SF & 0.0 & 1.041 & 0.967 & 0.994 & 0.989 & 0.957 & 0.993 \\ |
| 1411 |
< |
& 0.1 & 1.050 & 0.968 & 0.996 & 0.991 & 0.972 & 0.995 \\ |
| 1412 |
< |
& 0.2 & 0.982 & 0.975 & 0.959 & 0.980 & 0.960 & 0.980 \\ |
| 1413 |
< |
& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.944 \\ |
| 1414 |
< |
\midrule |
| 1415 |
< |
PC & & 0.795 & 0.000 & 0.792 & 0.000 & 0.793 & 0.000 \\ |
| 1416 |
< |
SP & 0.0 & 0.916 & 0.829 & 1.086 & 0.791 & 1.010 & 0.936 \\ |
| 1417 |
< |
& 0.1 & 0.958 & 0.917 & 1.049 & 0.943 & 1.001 & 0.995 \\ |
| 1418 |
< |
& 0.2 & 0.981 & 0.981 & 0.982 & 0.984 & 0.981 & 0.984 \\ |
| 1419 |
< |
& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ |
| 1420 |
< |
SF & 0.0 & 1.002 & 0.983 & 0.997 & 0.994 & 0.991 & 0.997 \\ |
| 1421 |
< |
& 0.1 & 1.003 & 0.984 & 0.996 & 0.995 & 0.993 & 0.997 \\ |
| 1422 |
< |
& 0.2 & 0.983 & 0.980 & 0.981 & 0.984 & 0.981 & 0.984 \\ |
| 1423 |
< |
& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ |
| 1424 |
< |
\bottomrule |
| 1425 |
< |
\end{tabular} |
| 1426 |
< |
\label{tab:salt} |
| 1427 |
< |
\end{table} |
| 1428 |
< |
|
| 1429 |
< |
\begin{table}[htbp] |
| 1430 |
< |
\centering |
| 1431 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1432 |
< |
DISTRIBUTIONS OF THE FORCE VECTORS IN THE CRYSTALLINE SODIUM CHLORIDE |
| 1433 |
< |
SYSTEM} |
| 1434 |
< |
|
| 1435 |
< |
\footnotesize |
| 1436 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1437 |
< |
\toprule |
| 1438 |
< |
\toprule |
| 1439 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1440 |
< |
\cmidrule(lr){3-5} |
| 1441 |
< |
\cmidrule(l){6-8} |
| 1442 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ |
| 1443 |
< |
\midrule |
| 1444 |
< |
PC & & 111.945 & 111.824 & 111.866 \\ |
| 1445 |
< |
SP & 0.0 & 112.414 & 152.215 & 38.087 \\ |
| 1446 |
< |
& 0.1 & 52.361 & 42.574 & 2.819 \\ |
| 1447 |
< |
& 0.2 & 10.847 & 9.709 & 9.686 \\ |
| 1448 |
< |
& 0.3 & 31.128 & 31.104 & 31.029 \\ |
| 1449 |
< |
SF & 0.0 & 10.025 & 3.555 & 1.648 \\ |
| 1450 |
< |
& 0.1 & 9.462 & 3.303 & 1.721 \\ |
| 1451 |
< |
& 0.2 & 11.454 & 9.813 & 9.701 \\ |
| 1452 |
< |
& 0.3 & 31.120 & 31.105 & 31.029 \\ |
| 1453 |
< |
\bottomrule |
| 1454 |
< |
\end{tabular} |
| 1455 |
< |
\label{tab:saltAng} |
| 1456 |
< |
\end{table} |
| 1457 |
< |
|
| 1458 |
< |
The crystalline NaCl system is the most challenging test case for the |
| 1459 |
< |
pairwise summation methods, as evidenced by the results in tables |
| 1460 |
< |
\ref{tab:salt} and \ref{tab:saltAng}. The undamped and weakly damped |
| 1461 |
< |
{\sc sf} methods seem to be the best choices. These methods match well |
| 1462 |
< |
with {\sc spme} across the energy gap, force magnitude, and force |
| 1463 |
< |
directionality tests. The {\sc sp} method struggles in all cases, |
| 1464 |
< |
with the exception of good dynamics reproduction when using weak |
| 1465 |
< |
electrostatic damping with a large cutoff radius. |
| 1466 |
< |
|
| 1467 |
< |
The moderate electrostatic damping case is not as good as we would |
| 1468 |
< |
expect given the long-time dynamics results observed for this system |
| 1469 |
< |
(see section \ref{sec:LongTimeDynamics}). Since the data tabulated in |
| 1470 |
< |
tables \ref{tab:salt} and \ref{tab:saltAng} are a test of |
| 1471 |
< |
instantaneous dynamics, this indicates that good long-time dynamics |
| 1472 |
< |
comes in part at the expense of short-time dynamics. |
| 1473 |
< |
|
| 1474 |
< |
\subsection{0.11M NaCl Solution Results} |
| 1475 |
< |
|
| 1476 |
< |
In an effort to bridge the charged atomic and neutral molecular |
| 1477 |
< |
systems, Na$^+$ and Cl$^-$ ion charge defects were incorporated into |
| 1478 |
< |
the liquid water system. This low ionic strength system consists of 4 |
| 1479 |
< |
ions in the 1000 SPC/E water solvent ($\approx$0.11 M). The results |
| 1480 |
< |
for the energy gap comparisons and the force and torque vector |
| 1481 |
< |
magnitude comparisons are shown in table \ref{tab:solnWeak}. The |
| 1482 |
< |
force and torque vector directionality results are displayed |
| 1483 |
< |
separately in table \ref{tab:solnWeakAng}, where the effect of |
| 1484 |
< |
group-based cutoffs and switching functions on the {\sc sp} and {\sc |
| 1485 |
< |
sf} potentials are investigated. |
| 1486 |
< |
|
| 1487 |
< |
\begin{table}[htbp] |
| 1488 |
< |
\centering |
| 1489 |
< |
\caption{REGRESSION RESULTS OF THE WEAK SODIUM CHLORIDE SOLUTION |
| 1490 |
< |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES |
| 1491 |
< |
({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1492 |
< |
|
| 1493 |
< |
\footnotesize |
| 1494 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1495 |
< |
\toprule |
| 1496 |
< |
\toprule |
| 1497 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1498 |
< |
\cmidrule(lr){3-4} |
| 1499 |
< |
\cmidrule(lr){5-6} |
| 1500 |
< |
\cmidrule(l){7-8} |
| 1501 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1502 |
< |
\midrule |
| 1503 |
< |
PC & & 0.247 & 0.000 & -1.103 & 0.001 & 5.480 & 0.015 \\ |
| 1504 |
< |
SP & 0.0 & 0.935 & 0.388 & 0.984 & 0.541 & 1.010 & 0.685 \\ |
| 1505 |
< |
& 0.1 & 0.951 & 0.603 & 0.993 & 0.875 & 1.001 & 0.979 \\ |
| 1506 |
< |
& 0.2 & 0.969 & 0.968 & 0.996 & 0.997 & 0.994 & 0.997 \\ |
| 1507 |
< |
& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ |
| 1508 |
< |
SF & 0.0 & 0.963 & 0.971 & 0.989 & 0.996 & 0.991 & 0.998 \\ |
| 1509 |
< |
& 0.1 & 0.970 & 0.971 & 0.995 & 0.997 & 0.997 & 0.999 \\ |
| 1510 |
< |
& 0.2 & 0.972 & 0.975 & 0.996 & 0.997 & 0.994 & 0.997 \\ |
| 1511 |
< |
& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ |
| 1512 |
< |
GSC & & 0.964 & 0.731 & 0.984 & 0.704 & 1.005 & 0.770 \\ |
| 1513 |
< |
RF & & 0.968 & 0.605 & 0.974 & 0.541 & 1.014 & 0.614 \\ |
| 1514 |
< |
\midrule |
| 1515 |
< |
PC & & 1.354 & 0.000 & -1.190 & 0.000 & -0.314 & 0.000 \\ |
| 1516 |
< |
SP & 0.0 & 0.720 & 0.338 & 0.808 & 0.523 & 0.860 & 0.643 \\ |
| 1517 |
< |
& 0.1 & 0.839 & 0.583 & 0.955 & 0.882 & 0.992 & 0.978 \\ |
| 1518 |
< |
& 0.2 & 0.995 & 0.987 & 0.999 & 1.000 & 0.999 & 1.000 \\ |
| 1519 |
< |
& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ |
| 1520 |
< |
SF & 0.0 & 0.998 & 0.994 & 1.000 & 0.998 & 1.000 & 0.999 \\ |
| 1521 |
< |
& 0.1 & 0.997 & 0.994 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1522 |
< |
& 0.2 & 0.999 & 0.998 & 0.999 & 1.000 & 0.999 & 1.000 \\ |
| 1523 |
< |
& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ |
| 1524 |
< |
GSC & & 0.995 & 0.990 & 0.998 & 0.997 & 0.998 & 0.996 \\ |
| 1525 |
< |
RF & & 0.998 & 0.993 & 0.999 & 0.998 & 0.999 & 0.996 \\ |
| 1526 |
< |
\midrule |
| 1527 |
< |
PC & & 2.437 & 0.000 & -1.872 & 0.000 & 2.138 & 0.000 \\ |
| 1528 |
< |
SP & 0.0 & 0.838 & 0.525 & 0.901 & 0.686 & 0.932 & 0.779 \\ |
| 1529 |
< |
& 0.1 & 0.914 & 0.733 & 0.979 & 0.932 & 0.995 & 0.987 \\ |
| 1530 |
< |
& 0.2 & 0.977 & 0.969 & 0.988 & 0.990 & 0.989 & 0.990 \\ |
| 1531 |
< |
& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ |
| 1532 |
< |
SF & 0.0 & 0.969 & 0.977 & 0.987 & 0.996 & 0.993 & 0.998 \\ |
| 1533 |
< |
& 0.1 & 0.975 & 0.978 & 0.993 & 0.996 & 0.997 & 0.998 \\ |
| 1534 |
< |
& 0.2 & 0.976 & 0.973 & 0.988 & 0.990 & 0.989 & 0.990 \\ |
| 1535 |
< |
& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ |
| 1536 |
< |
GSC & & 0.980 & 0.959 & 0.990 & 0.983 & 0.992 & 0.989 \\ |
| 1537 |
< |
RF & & 0.984 & 0.975 & 0.996 & 0.995 & 0.998 & 0.998 \\ |
| 1538 |
< |
\bottomrule |
| 1539 |
< |
\end{tabular} |
| 1540 |
< |
\label{tab:solnWeak} |
| 1541 |
< |
\end{table} |
| 1542 |
< |
|
| 1543 |
< |
\begin{table}[htbp] |
| 1544 |
< |
\centering |
| 1545 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1546 |
< |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE WEAK SODIUM |
| 1547 |
< |
CHLORIDE SOLUTION SYSTEM} |
| 1548 |
< |
|
| 1549 |
< |
\footnotesize |
| 1550 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1551 |
< |
\toprule |
| 1552 |
< |
\toprule |
| 1553 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1554 |
< |
\cmidrule(lr){3-5} |
| 1555 |
< |
\cmidrule(l){6-8} |
| 1556 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1557 |
< |
\midrule |
| 1558 |
< |
PC & & 882.863 & 510.435 & 344.201 & 277.691 & 154.231 & 100.131 \\ |
| 1559 |
< |
SP & 0.0 & 732.569 & 405.704 & 257.756 & 261.445 & 142.245 & 91.497 \\ |
| 1560 |
< |
& 0.1 & 329.031 & 70.746 & 12.014 & 118.496 & 25.218 & 4.711 \\ |
| 1561 |
< |
& 0.2 & 6.772 & 0.153 & 0.118 & 9.780 & 2.101 & 2.102 \\ |
| 1562 |
< |
& 0.3 & 0.951 & 0.774 & 0.784 & 12.108 & 7.673 & 7.851 \\ |
| 1563 |
< |
SF & 0.0 & 2.555 & 0.762 & 0.313 & 6.590 & 1.328 & 0.558 \\ |
| 1564 |
< |
& 0.1 & 2.561 & 0.560 & 0.123 & 6.464 & 1.162 & 0.457 \\ |
| 1565 |
< |
& 0.2 & 0.501 & 0.118 & 0.118 & 5.698 & 2.074 & 2.099 \\ |
| 1566 |
< |
& 0.3 & 0.943 & 0.774 & 0.784 & 12.118 & 7.674 & 7.851 \\ |
| 1567 |
< |
GSC & & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ |
| 1568 |
< |
RF & & 2.415 & 0.452 & 0.130 & 6.915 & 1.423 & 0.507 \\ |
| 1569 |
< |
\midrule |
| 1570 |
< |
GSSP & 0.0 & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ |
| 1571 |
< |
& 0.1 & 2.251 & 0.324 & 0.064 & 7.628 & 1.639 & 0.497 \\ |
| 1572 |
< |
& 0.2 & 0.590 & 0.118 & 0.116 & 6.080 & 2.096 & 2.103 \\ |
| 1573 |
< |
& 0.3 & 0.953 & 0.759 & 0.780 & 12.347 & 7.683 & 7.849 \\ |
| 1574 |
< |
GSSF & 0.0 & 1.541 & 0.301 & 0.096 & 6.407 & 1.316 & 0.496 \\ |
| 1575 |
< |
& 0.1 & 1.541 & 0.237 & 0.050 & 6.356 & 1.202 & 0.457 \\ |
| 1576 |
< |
& 0.2 & 0.568 & 0.118 & 0.116 & 6.166 & 2.105 & 2.105 \\ |
| 1577 |
< |
& 0.3 & 0.954 & 0.759 & 0.780 & 12.337 & 7.684 & 7.849 \\ |
| 1578 |
< |
\bottomrule |
| 1579 |
< |
\end{tabular} |
| 1580 |
< |
\label{tab:solnWeakAng} |
| 1581 |
< |
\end{table} |
| 1582 |
< |
|
| 1583 |
< |
Because this system is a perturbation of the pure liquid water system, |
| 1584 |
< |
comparisons are best drawn between these two sets. The {\sc sp} and |
| 1585 |
< |
{\sc sf} methods are not significantly affected by the inclusion of a |
| 1586 |
< |
few ions. The aspect of cutoff sphere neutralization aids in the |
| 1587 |
< |
smooth incorporation of these ions; thus, all of the observations |
| 1588 |
< |
regarding these methods carry over from section |
| 1589 |
< |
\ref{sec:WaterResults}. The differences between these systems are more |
| 1590 |
< |
visible for the {\sc rf} method. Though good force agreement is still |
| 1591 |
< |
maintained, the energy gaps show a significant increase in the scatter |
| 1592 |
< |
of the data. |
| 1593 |
< |
|
| 1594 |
< |
\subsection{1.1M NaCl Solution Results} |
| 1595 |
< |
|
| 1596 |
< |
The bridging of the charged atomic and neutral molecular systems was |
| 1597 |
< |
further developed by considering a high ionic strength system |
| 1598 |
< |
consisting of 40 ions in the 1000 SPC/E water solvent ($\approx$1.1 |
| 1599 |
< |
M). The results for the energy gap comparisons and the force and |
| 1600 |
< |
torque vector magnitude comparisons are shown in table |
| 1601 |
< |
\ref{tab:solnStr}. The force and torque vector directionality |
| 1602 |
< |
results are displayed separately in table \ref{tab:solnStrAng}, where |
| 1603 |
< |
the effect of group-based cutoffs and switching functions on the {\sc |
| 1604 |
< |
sp} and {\sc sf} potentials are investigated. |
| 1605 |
< |
|
| 1606 |
< |
\begin{table}[htbp] |
| 1607 |
< |
\centering |
| 1608 |
< |
\caption{REGRESSION RESULTS OF THE STRONG SODIUM CHLORIDE SOLUTION |
| 1609 |
< |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES |
| 1610 |
< |
({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1611 |
< |
|
| 1612 |
< |
\footnotesize |
| 1613 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1614 |
< |
\toprule |
| 1615 |
< |
\toprule |
| 1616 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1617 |
< |
\cmidrule(lr){3-4} |
| 1618 |
< |
\cmidrule(lr){5-6} |
| 1619 |
< |
\cmidrule(l){7-8} |
| 1620 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1621 |
< |
\midrule |
| 1622 |
< |
PC & & -0.081 & 0.000 & 0.945 & 0.001 & 0.073 & 0.000 \\ |
| 1623 |
< |
SP & 0.0 & 0.978 & 0.469 & 0.996 & 0.672 & 0.975 & 0.668 \\ |
| 1624 |
< |
& 0.1 & 0.944 & 0.645 & 0.997 & 0.886 & 0.991 & 0.978 \\ |
| 1625 |
< |
& 0.2 & 0.873 & 0.896 & 0.985 & 0.993 & 0.980 & 0.993 \\ |
| 1626 |
< |
& 0.3 & 0.831 & 0.860 & 0.960 & 0.979 & 0.955 & 0.977 \\ |
| 1627 |
< |
SF & 0.0 & 0.858 & 0.905 & 0.985 & 0.970 & 0.990 & 0.998 \\ |
| 1628 |
< |
& 0.1 & 0.865 & 0.907 & 0.992 & 0.974 & 0.994 & 0.999 \\ |
| 1629 |
< |
& 0.2 & 0.862 & 0.894 & 0.985 & 0.993 & 0.980 & 0.993 \\ |
| 1630 |
< |
& 0.3 & 0.831 & 0.859 & 0.960 & 0.979 & 0.955 & 0.977 \\ |
| 1631 |
< |
GSC & & 1.985 & 0.152 & 0.760 & 0.031 & 1.106 & 0.062 \\ |
| 1632 |
< |
RF & & 2.414 & 0.116 & 0.813 & 0.017 & 1.434 & 0.047 \\ |
| 1633 |
< |
\midrule |
| 1634 |
< |
PC & & -7.028 & 0.000 & -9.364 & 0.000 & 0.925 & 0.865 \\ |
| 1635 |
< |
SP & 0.0 & 0.701 & 0.319 & 0.909 & 0.773 & 0.861 & 0.665 \\ |
| 1636 |
< |
& 0.1 & 0.824 & 0.565 & 0.970 & 0.930 & 0.990 & 0.979 \\ |
| 1637 |
< |
& 0.2 & 0.988 & 0.981 & 0.995 & 0.998 & 0.991 & 0.998 \\ |
| 1638 |
< |
& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ |
| 1639 |
< |
SF & 0.0 & 0.993 & 0.988 & 0.992 & 0.984 & 0.998 & 0.999 \\ |
| 1640 |
< |
& 0.1 & 0.993 & 0.989 & 0.993 & 0.986 & 0.998 & 1.000 \\ |
| 1641 |
< |
& 0.2 & 0.993 & 0.992 & 0.995 & 0.998 & 0.991 & 0.998 \\ |
| 1642 |
< |
& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ |
| 1643 |
< |
GSC & & 0.964 & 0.897 & 0.970 & 0.917 & 0.925 & 0.865 \\ |
| 1644 |
< |
RF & & 0.994 & 0.864 & 0.988 & 0.865 & 0.980 & 0.784 \\ |
| 1645 |
< |
\midrule |
| 1646 |
< |
PC & & -2.212 & 0.000 & -0.588 & 0.000 & 0.953 & 0.925 \\ |
| 1647 |
< |
SP & 0.0 & 0.800 & 0.479 & 0.930 & 0.804 & 0.924 & 0.759 \\ |
| 1648 |
< |
& 0.1 & 0.883 & 0.694 & 0.976 & 0.942 & 0.993 & 0.986 \\ |
| 1649 |
< |
& 0.2 & 0.952 & 0.943 & 0.980 & 0.984 & 0.980 & 0.983 \\ |
| 1650 |
< |
& 0.3 & 0.914 & 0.909 & 0.943 & 0.948 & 0.944 & 0.946 \\ |
| 1651 |
< |
SF & 0.0 & 0.945 & 0.953 & 0.980 & 0.984 & 0.991 & 0.998 \\ |
| 1652 |
< |
& 0.1 & 0.951 & 0.954 & 0.987 & 0.986 & 0.995 & 0.998 \\ |
| 1653 |
< |
& 0.2 & 0.951 & 0.946 & 0.980 & 0.984 & 0.980 & 0.983 \\ |
| 1654 |
< |
& 0.3 & 0.914 & 0.908 & 0.943 & 0.948 & 0.944 & 0.946 \\ |
| 1655 |
< |
GSC & & 0.882 & 0.818 & 0.939 & 0.902 & 0.953 & 0.925 \\ |
| 1656 |
< |
RF & & 0.949 & 0.939 & 0.988 & 0.988 & 0.992 & 0.993 \\ |
| 1657 |
< |
\bottomrule |
| 1658 |
< |
\end{tabular} |
| 1659 |
< |
\label{tab:solnStr} |
| 1660 |
< |
\end{table} |
| 1661 |
< |
|
| 1662 |
< |
\begin{table}[htbp] |
| 1663 |
< |
\centering |
| 1664 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1665 |
< |
OF THE FORCE AND TORQUE VECTORS IN THE STRONG SODIUM CHLORIDE SOLUTION |
| 1666 |
< |
SYSTEM} |
| 1667 |
< |
|
| 1668 |
< |
\footnotesize |
| 1669 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1670 |
< |
\toprule |
| 1671 |
< |
\toprule |
| 1672 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1673 |
< |
\cmidrule(lr){3-5} |
| 1674 |
< |
\cmidrule(l){6-8} |
| 1675 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1676 |
< |
\midrule |
| 1677 |
< |
PC & & 957.784 & 513.373 & 2.260 & 340.043 & 179.443 & 13.079 \\ |
| 1678 |
< |
SP & 0.0 & 786.244 & 139.985 & 259.289 & 311.519 & 90.280 & 105.187 \\ |
| 1679 |
< |
& 0.1 & 354.697 & 38.614 & 12.274 & 144.531 & 23.787 & 5.401 \\ |
| 1680 |
< |
& 0.2 & 7.674 & 0.363 & 0.215 & 16.655 & 3.601 & 3.634 \\ |
| 1681 |
< |
& 0.3 & 1.745 & 1.456 & 1.449 & 23.669 & 14.376 & 14.240 \\ |
| 1682 |
< |
SF & 0.0 & 3.282 & 8.567 & 0.369 & 11.904 & 6.589 & 0.717 \\ |
| 1683 |
< |
& 0.1 & 3.263 & 7.479 & 0.142 & 11.634 & 5.750 & 0.591 \\ |
| 1684 |
< |
& 0.2 & 0.686 & 0.324 & 0.215 & 10.809 & 3.580 & 3.635 \\ |
| 1685 |
< |
& 0.3 & 1.749 & 1.456 & 1.449 & 23.635 & 14.375 & 14.240 \\ |
| 1686 |
< |
GSC & & 6.181 & 2.904 & 2.263 & 44.349 & 19.442 & 12.873 \\ |
| 1687 |
< |
RF & & 3.891 & 0.847 & 0.323 & 18.628 & 3.995 & 2.072 \\ |
| 1688 |
< |
\midrule |
| 1689 |
< |
GSSP & 0.0 & 6.197 & 2.929 & 2.290 & 44.441 & 19.442 & 12.873 \\ |
| 1690 |
< |
& 0.1 & 4.688 & 1.064 & 0.260 & 31.208 & 6.967 & 2.303 \\ |
| 1691 |
< |
& 0.2 & 1.021 & 0.218 & 0.213 & 14.425 & 3.629 & 3.649 \\ |
| 1692 |
< |
& 0.3 & 1.752 & 1.454 & 1.451 & 23.540 & 14.390 & 14.245 \\ |
| 1693 |
< |
GSSF & 0.0 & 2.494 & 0.546 & 0.217 & 16.391 & 3.230 & 1.613 \\ |
| 1694 |
< |
& 0.1 & 2.448 & 0.429 & 0.106 & 16.390 & 2.827 & 1.159 \\ |
| 1695 |
< |
& 0.2 & 0.899 & 0.214 & 0.213 & 13.542 & 3.583 & 3.645 \\ |
| 1696 |
< |
& 0.3 & 1.752 & 1.454 & 1.451 & 23.587 & 14.390 & 14.245 \\ |
| 1697 |
< |
\bottomrule |
| 1698 |
< |
\end{tabular} |
| 1699 |
< |
\label{tab:solnStrAng} |
| 1700 |
< |
\end{table} |
| 1701 |
< |
|
| 1702 |
< |
The {\sc rf} method struggles with the jump in ionic strength. The |
| 1703 |
< |
configuration energy differences degrade to unusable levels while the |
| 1704 |
< |
forces and torques show a more modest reduction in the agreement with |
| 1705 |
< |
{\sc spme}. The {\sc rf} method was designed for homogeneous systems, |
| 1706 |
< |
and this attribute is apparent in these results. |
| 1707 |
< |
|
| 1708 |
< |
The {\sc sp} and {\sc sf} methods require larger cutoffs to maintain |
| 1709 |
< |
their agreement with {\sc spme}. With these results, we still |
| 1710 |
< |
recommend undamped to moderate damping for the {\sc sf} method and |
| 1711 |
< |
moderate damping for the {\sc sp} method, both with cutoffs greater |
| 1712 |
< |
than 12\AA. |
| 1713 |
< |
|
| 1714 |
< |
\subsection{6\AA\ Argon Sphere in SPC/E Water Results} |
| 1715 |
< |
|
| 1716 |
< |
The final model system studied was a 6\AA\ sphere of Argon solvated |
| 1717 |
< |
by SPC/E water. This serves as a test case of a specifically sized |
| 1718 |
< |
electrostatic defect in a disordered molecular system. The results for |
| 1719 |
< |
the energy gap comparisons and the force and torque vector magnitude |
| 1720 |
< |
comparisons are shown in table \ref{tab:argon}. The force and torque |
| 1721 |
< |
vector directionality results are displayed separately in table |
| 1722 |
< |
\ref{tab:argonAng}, where the effect of group-based cutoffs and |
| 1723 |
< |
switching functions on the {\sc sp} and {\sc sf} potentials are |
| 1724 |
< |
investigated. |
| 1725 |
< |
|
| 1726 |
< |
\begin{table}[htbp] |
| 1727 |
< |
\centering |
| 1728 |
< |
\caption{REGRESSION RESULTS OF THE 6\AA\ ARGON SPHERE IN LIQUID |
| 1729 |
< |
WATER SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR |
| 1730 |
< |
MAGNITUDES ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1731 |
< |
|
| 1732 |
< |
\footnotesize |
| 1733 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1734 |
< |
\toprule |
| 1735 |
< |
\toprule |
| 1736 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1737 |
< |
\cmidrule(lr){3-4} |
| 1738 |
< |
\cmidrule(lr){5-6} |
| 1739 |
< |
\cmidrule(l){7-8} |
| 1740 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1741 |
< |
\midrule |
| 1742 |
< |
PC & & 2.320 & 0.008 & -0.650 & 0.001 & 3.848 & 0.029 \\ |
| 1743 |
< |
SP & 0.0 & 1.053 & 0.711 & 0.977 & 0.820 & 0.974 & 0.882 \\ |
| 1744 |
< |
& 0.1 & 1.032 & 0.846 & 0.989 & 0.965 & 0.992 & 0.994 \\ |
| 1745 |
< |
& 0.2 & 0.993 & 0.995 & 0.982 & 0.998 & 0.986 & 0.998 \\ |
| 1746 |
< |
& 0.3 & 0.968 & 0.995 & 0.954 & 0.992 & 0.961 & 0.994 \\ |
| 1747 |
< |
SF & 0.0 & 0.982 & 0.996 & 0.992 & 0.999 & 0.993 & 1.000 \\ |
| 1748 |
< |
& 0.1 & 0.987 & 0.996 & 0.996 & 0.999 & 0.997 & 1.000 \\ |
| 1749 |
< |
& 0.2 & 0.989 & 0.998 & 0.984 & 0.998 & 0.989 & 0.998 \\ |
| 1750 |
< |
& 0.3 & 0.971 & 0.995 & 0.957 & 0.992 & 0.965 & 0.994 \\ |
| 1751 |
< |
GSC & & 1.002 & 0.983 & 0.992 & 0.973 & 0.996 & 0.971 \\ |
| 1752 |
< |
RF & & 0.998 & 0.995 & 0.999 & 0.998 & 0.998 & 0.998 \\ |
| 1753 |
< |
\midrule |
| 1754 |
< |
PC & & -36.559 & 0.002 & -44.917 & 0.004 & -52.945 & 0.006 \\ |
| 1755 |
< |
SP & 0.0 & 0.890 & 0.786 & 0.927 & 0.867 & 0.949 & 0.909 \\ |
| 1756 |
< |
& 0.1 & 0.942 & 0.895 & 0.984 & 0.974 & 0.997 & 0.995 \\ |
| 1757 |
< |
& 0.2 & 0.999 & 0.997 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1758 |
< |
& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ |
| 1759 |
< |
SF & 0.0 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1760 |
< |
& 0.1 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1761 |
< |
& 0.2 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1762 |
< |
& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ |
| 1763 |
< |
GSC & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1764 |
< |
RF & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1765 |
< |
\midrule |
| 1766 |
< |
PC & & 1.984 & 0.000 & 0.012 & 0.000 & 1.357 & 0.000 \\ |
| 1767 |
< |
SP & 0.0 & 0.850 & 0.552 & 0.907 & 0.703 & 0.938 & 0.793 \\ |
| 1768 |
< |
& 0.1 & 0.924 & 0.755 & 0.980 & 0.936 & 0.995 & 0.988 \\ |
| 1769 |
< |
& 0.2 & 0.985 & 0.983 & 0.986 & 0.988 & 0.987 & 0.988 \\ |
| 1770 |
< |
& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ |
| 1771 |
< |
SF & 0.0 & 0.977 & 0.989 & 0.987 & 0.995 & 0.992 & 0.998 \\ |
| 1772 |
< |
& 0.1 & 0.982 & 0.989 & 0.992 & 0.996 & 0.997 & 0.998 \\ |
| 1773 |
< |
& 0.2 & 0.984 & 0.987 & 0.986 & 0.987 & 0.987 & 0.988 \\ |
| 1774 |
< |
& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ |
| 1775 |
< |
GSC & & 0.995 & 0.981 & 0.999 & 0.990 & 1.000 & 0.993 \\ |
| 1776 |
< |
RF & & 0.993 & 0.988 & 0.997 & 0.995 & 0.999 & 0.998 \\ |
| 1777 |
< |
\bottomrule |
| 1778 |
< |
\end{tabular} |
| 1779 |
< |
\label{tab:argon} |
| 1780 |
< |
\end{table} |
| 1781 |
< |
|
| 1782 |
< |
\begin{table}[htbp] |
| 1783 |
< |
\centering |
| 1784 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1785 |
< |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6\AA\ SPHERE OF |
| 1786 |
< |
ARGON IN LIQUID WATER SYSTEM} |
| 1787 |
< |
|
| 1788 |
< |
\footnotesize |
| 1789 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1790 |
< |
\toprule |
| 1791 |
< |
\toprule |
| 1792 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1793 |
< |
\cmidrule(lr){3-5} |
| 1794 |
< |
\cmidrule(l){6-8} |
| 1795 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1796 |
< |
\midrule |
| 1797 |
< |
PC & & 568.025 & 265.993 & 195.099 & 246.626 & 138.600 & 91.654 \\ |
| 1798 |
< |
SP & 0.0 & 504.578 & 251.694 & 179.932 & 231.568 & 131.444 & 85.119 \\ |
| 1799 |
< |
& 0.1 & 224.886 & 49.746 & 9.346 & 104.482 & 23.683 & 4.480 \\ |
| 1800 |
< |
& 0.2 & 4.889 & 0.197 & 0.155 & 6.029 & 2.507 & 2.269 \\ |
| 1801 |
< |
& 0.3 & 0.817 & 0.833 & 0.812 & 8.286 & 8.436 & 8.135 \\ |
| 1802 |
< |
SF & 0.0 & 1.924 & 0.675 & 0.304 & 3.658 & 1.448 & 0.600 \\ |
| 1803 |
< |
& 0.1 & 1.937 & 0.515 & 0.143 & 3.565 & 1.308 & 0.546 \\ |
| 1804 |
< |
& 0.2 & 0.407 & 0.166 & 0.156 & 3.086 & 2.501 & 2.274 \\ |
| 1805 |
< |
& 0.3 & 0.815 & 0.833 & 0.812 & 8.330 & 8.437 & 8.135 \\ |
| 1806 |
< |
GSC & & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ |
| 1807 |
< |
RF & & 1.822 & 0.408 & 0.142 & 3.799 & 1.362 & 0.550 \\ |
| 1808 |
< |
\midrule |
| 1809 |
< |
GSSP & 0.0 & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ |
| 1810 |
< |
& 0.1 & 1.652 & 0.309 & 0.087 & 4.197 & 1.401 & 0.590 \\ |
| 1811 |
< |
& 0.2 & 0.465 & 0.165 & 0.153 & 3.323 & 2.529 & 2.273 \\ |
| 1812 |
< |
& 0.3 & 0.813 & 0.825 & 0.816 & 8.316 & 8.447 & 8.132 \\ |
| 1813 |
< |
GSSF & 0.0 & 1.173 & 0.292 & 0.113 & 3.452 & 1.347 & 0.583 \\ |
| 1814 |
< |
& 0.1 & 1.166 & 0.240 & 0.076 & 3.381 & 1.281 & 0.575 \\ |
| 1815 |
< |
& 0.2 & 0.459 & 0.165 & 0.153 & 3.430 & 2.542 & 2.273 \\ |
| 1816 |
< |
& 0.3 & 0.814 & 0.825 & 0.816 & 8.325 & 8.447 & 8.132 \\ |
| 1817 |
< |
\bottomrule |
| 1818 |
< |
\end{tabular} |
| 1819 |
< |
\label{tab:argonAng} |
| 1820 |
< |
\end{table} |
| 1821 |
< |
|
| 1822 |
< |
This system does not appear to show any significant deviations from |
| 1823 |
< |
the previously observed results. The {\sc sp} and {\sc sf} methods |
| 1824 |
< |
have agreements similar to those observed in section |
| 1825 |
< |
\ref{sec:WaterResults}. The only significant difference is the |
| 1826 |
< |
improvement in the configuration energy differences for the {\sc rf} |
| 1827 |
< |
method. This is surprising in that we are introducing an inhomogeneity |
| 1828 |
< |
to the system; however, this inhomogeneity is charge-neutral and does |
| 1829 |
< |
not result in charged cutoff spheres. The charge-neutrality of the |
| 1830 |
< |
cutoff spheres, which the {\sc sp} and {\sc sf} methods explicitly |
| 1831 |
< |
enforce, seems to play a greater role in the stability of the {\sc rf} |
| 1832 |
< |
method than the required homogeneity of the environment. |
| 1833 |
< |
|
| 1834 |
< |
|
| 1835 |
< |
\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics} |
| 1836 |
< |
|
| 1837 |
< |
Zahn {\it et al.} investigated the structure and dynamics of water |
| 1838 |
< |
using eqs. (\ref{eq:ZahnPot}) and |
| 1839 |
< |
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 1840 |
< |
that a method similar (but not identical with) the damped {\sc sf} |
| 1841 |
< |
method resulted in properties very similar to those obtained when |
| 1842 |
< |
using the Ewald summation. The properties they studied (pair |
| 1843 |
< |
distribution functions, diffusion constants, and velocity and |
| 1844 |
< |
orientational correlation functions) may not be particularly sensitive |
| 1845 |
< |
to the long-range and collective behavior that governs the |
| 1846 |
< |
low-frequency behavior in crystalline systems. Additionally, the |
| 1847 |
< |
ionic crystals are the worst case scenario for the pairwise methods |
| 1848 |
< |
because they lack the reciprocal space contribution contained in the |
| 1849 |
< |
Ewald summation. |
| 1850 |
< |
|
| 1851 |
< |
We are using two separate measures to probe the effects of these |
| 1852 |
< |
alternative electrostatic methods on the dynamics in crystalline |
| 1853 |
< |
materials. For short- and intermediate-time dynamics, we are |
| 1854 |
< |
computing the velocity autocorrelation function, and for long-time |
| 1855 |
< |
and large length-scale collective motions, we are looking at the |
| 1856 |
< |
low-frequency portion of the power spectrum. |
| 1857 |
< |
|
| 1858 |
< |
\begin{figure} |
| 1859 |
< |
\centering |
| 1860 |
< |
\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf} |
| 1861 |
< |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1862 |
< |
1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& |
| 1863 |
< |
0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2\AA$^{-1}$). The inset is |
| 1864 |
< |
a magnification of the area around the first minimum. The times to |
| 1865 |
< |
first collision are nearly identical, but differences can be seen in |
| 1866 |
< |
the peaks and troughs, where the undamped and weakly damped methods |
| 1867 |
< |
are stiffer than the moderately damped and {\sc spme} methods.} |
| 1868 |
< |
\label{fig:vCorrPlot} |
| 1869 |
< |
\end{figure} |
| 1870 |
< |
|
| 1871 |
< |
The short-time decay of the velocity autocorrelation function through |
| 1872 |
< |
the first collision are nearly identical in figure |
| 1873 |
< |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1874 |
< |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1875 |
< |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1876 |
< |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1877 |
< |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1878 |
< |
approaches the {\sc spme} results. With $\alpha$ values of 0.2\AA$^{-1}$, |
| 1879 |
< |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1880 |
< |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1881 |
< |
and {\sc sp} potentials become nearly identical with increased |
| 1882 |
< |
damping. However, this appears to indicate that once damping is |
| 1883 |
< |
utilized, the details of the form of the potential (and forces) |
| 1884 |
< |
constructed out of the damped electrostatic interaction are less |
| 1885 |
< |
important. |
| 1886 |
< |
|
| 1887 |
< |
\section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics} |
| 1888 |
< |
|
| 1889 |
< |
To evaluate how the differences between the methods affect the |
| 1890 |
< |
collective long-time motion, we computed power spectra from long-time |
| 1891 |
< |
traces of the velocity autocorrelation function. The power spectra for |
| 1892 |
< |
the best-performing alternative methods are shown in |
| 1893 |
< |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1894 |
< |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1895 |
< |
ringing resulting from data truncation. This procedure had no |
| 1896 |
< |
noticeable effect on peak location or magnitude. |
| 1897 |
< |
|
| 1898 |
< |
\begin{figure} |
| 1899 |
< |
\centering |
| 1900 |
< |
\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf} |
| 1901 |
< |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1902 |
< |
functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf} |
| 1903 |
< |
($\alpha$ = 0, 0.1, \& 0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = |
| 1904 |
< |
0.2\AA$^{-1}$). The inset shows the frequency region below 100 |
| 1905 |
< |
cm$^{-1}$ to highlight where the spectra differ.} |
| 1906 |
< |
\label{fig:methodPS} |
| 1907 |
< |
\end{figure} |
| 1908 |
< |
|
| 1909 |
< |
While the high frequency regions of the power spectra for the |
| 1910 |
< |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1911 |
< |
the low frequency region shows how the summation methods differ. |
| 1912 |
< |
Considering the low-frequency inset (expanded in the upper frame of |
| 1913 |
< |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1914 |
< |
correlated motions are blue-shifted when using undamped or weakly |
| 1915 |
< |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1916 |
< |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1917 |
< |
correlated motion to the Ewald method (which has a convergence |
| 1918 |
< |
parameter of 0.3119\AA$^{-1}$). This weakening of the electrostatic |
| 1919 |
< |
interaction with increased damping explains why the long-ranged |
| 1920 |
< |
correlated motions are at lower frequencies for the moderately damped |
| 1921 |
< |
methods than for undamped or weakly damped methods. |
| 1922 |
< |
|
| 1923 |
< |
To isolate the role of the damping constant, we have computed the |
| 1924 |
< |
spectra for a single method ({\sc sf}) with a range of damping |
| 1925 |
< |
constants and compared this with the {\sc spme} spectrum. |
| 1926 |
< |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1927 |
< |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1928 |
< |
However, even without any electrostatic damping, the {\sc sf} method |
| 1929 |
< |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1930 |
< |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1931 |
< |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1932 |
< |
Most} of the collective behavior in the crystal is accurately captured |
| 1933 |
< |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1934 |
< |
obtained using moderate damping in addition to the shifting at the |
| 1935 |
< |
cutoff distance. |
| 1936 |
< |
|
| 1937 |
< |
\begin{figure} |
| 1938 |
< |
\centering |
| 1939 |
< |
\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf} |
| 1940 |
< |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1941 |
< |
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1942 |
< |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1943 |
< |
electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement |
| 1944 |
< |
with the power spectrum obtained using the Ewald sum. Over-damping can |
| 1945 |
< |
result in underestimates of frequencies of the long-wavelength |
| 1946 |
< |
motions.} |
| 1947 |
< |
\label{fig:dampInc} |
| 1948 |
< |
\end{figure} |
| 1949 |
< |
|
| 1950 |
< |
\section{An Application: TIP5P-E Water}\label{sec:t5peApplied} |
| 1951 |
< |
|
| 1952 |
< |
The above sections focused on the energetics and dynamics of a variety |
| 1953 |
< |
of systems when utilizing the {\sc sp} and {\sc sf} pairwise |
| 1954 |
< |
techniques. A unitary correlation with results obtained using the |
| 1955 |
< |
Ewald summation should result in a successful reproduction of both the |
| 1956 |
< |
static and dynamic properties of any selected system. To test this, |
| 1957 |
< |
we decided to calculate a series of properties for the TIP5P-E water |
| 1958 |
< |
model when using the {\sc sf} technique. |
| 1959 |
< |
|
| 1960 |
< |
The TIP5P-E water model is a variant of Mahoney and Jorgensen's |
| 1961 |
< |
five-point transferable intermolecular potential (TIP5P) model for |
| 1962 |
< |
water.\cite{Mahoney00} TIP5P was developed to reproduce the density |
| 1963 |
< |
maximum anomaly present in liquid water near 4$^\circ$C. As with many |
| 1964 |
< |
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC, |
| 1965 |
< |
and SPC/E), TIP5P was parametrized using a simple cutoff with no |
| 1966 |
< |
long-range electrostatic |
| 1967 |
< |
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
| 1968 |
< |
Without this correction, the pressure term on the central particle |
| 1969 |
< |
from the surroundings is missing. Because they expand to compensate |
| 1970 |
< |
for this added pressure term when this correction is included, systems |
| 1971 |
< |
composed of these particles tend to underpredict the density of water |
| 1972 |
< |
under standard conditions. When using any form of long-range |
| 1973 |
< |
electrostatic correction, it has become common practice to develop or |
| 1974 |
< |
utilize a reparametrized water model that corrects for this |
| 1975 |
< |
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 1976 |
< |
this practice and was optimized specifically for use with the Ewald |
| 1977 |
< |
summation.\cite{Rick04} In his publication, Rick preserved the |
| 1978 |
< |
geometry and point charge magnitudes in TIP5P and focused on altering |
| 1979 |
< |
the Lennard-Jones parameters to correct the density at |
| 1980 |
< |
298K.\cite{Rick04} With the density corrected, he compared common |
| 1981 |
< |
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1982 |
< |
9\AA\ cutoff. |
| 1983 |
< |
|
| 1984 |
< |
In the following sections, we compared these same water properties |
| 1985 |
< |
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1986 |
< |
{\sc sf} technique. In the above evaluation of the pairwise |
| 1987 |
< |
techniques, we observed some flexibility in the choice of parameters. |
| 1988 |
< |
Because of this, the following comparisons include the {\sc sf} |
| 1989 |
< |
technique with a 12\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 1990 |
< |
0.2\AA$^{-1}$, as well as a 9\AA\ cutoff with an $\alpha$ = |
| 1991 |
< |
0.2\AA$^{-1}$. |
| 1992 |
< |
|
| 1993 |
< |
\subsection{Density}\label{sec:t5peDensity} |
| 1994 |
< |
|
| 1995 |
< |
As stated previously, the property that prompted the development of |
| 1996 |
< |
TIP5P-E was the density at 1 atm. The density depends upon the |
| 1997 |
< |
internal pressure of the system in the $NPT$ ensemble, and the |
| 1998 |
< |
calculation of the pressure includes a components from both the |
| 1999 |
< |
kinetic energy and the virial. More specifically, the instantaneous |
| 2000 |
< |
molecular pressure ($p(t)$) is given by |
| 2001 |
< |
\begin{equation} |
| 2002 |
< |
p(t) = \frac{1}{\textrm{d}V}\sum_\mu |
| 2003 |
< |
\left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} |
| 2004 |
< |
+ \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right], |
| 2005 |
< |
\label{eq:MolecularPressure} |
| 2006 |
< |
\end{equation} |
| 2007 |
< |
where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of |
| 2008 |
< |
molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass |
| 2009 |
< |
($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on |
| 2010 |
< |
atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the |
| 2011 |
< |
right term in the brackets of eq. \ref{eq:MolecularPressure}) is |
| 2012 |
< |
directly dependent on the interatomic forces. Since the {\sc sp} |
| 2013 |
< |
method does not modify the forces (see |
| 2014 |
< |
sec. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will |
| 2015 |
< |
be identical to that obtained without an electrostatic correction. |
| 2016 |
< |
The {\sc sf} method does alter the virial component and, by way of the |
| 2017 |
< |
modified pressures, should provide densities more in line with those |
| 2018 |
< |
obtained using the Ewald summation. |
| 2019 |
< |
|
| 2020 |
< |
To compare densities, $NPT$ simulations were performed with the same |
| 2021 |
< |
temperatures as those selected by Rick in his Ewald summation |
| 2022 |
< |
simulations.\cite{Rick04} In order to improve statistics around the |
| 2023 |
< |
density maximum, 3ns trajectories were accumulated at 0, 12.5, and |
| 2024 |
< |
25$^\circ$C, while 2ns trajectories were obtained at all other |
| 2025 |
< |
temperatures. The average densities were calculated from the later |
| 2026 |
< |
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 2027 |
< |
method for accumulating statistics, these sequences were spliced into |
| 2028 |
< |
200 segments to calculate the average density and standard deviation |
| 2029 |
< |
at each temperature.\cite{Mahoney00} |
| 2030 |
< |
|
| 2031 |
< |
\begin{figure} |
| 2032 |
< |
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 2033 |
< |
\caption{Density versus temperature for the TIP5P-E water model when |
| 2034 |
< |
using the Ewald summation (Ref. \cite{Rick04}) and the {\sc sf} method |
| 2035 |
< |
with various parameters. The pressure term from the image-charge shell |
| 2036 |
< |
is larger than that provided by the reciprocal-space portion of the |
| 2037 |
< |
Ewald summation, leading to slightly lower densities. This effect is |
| 2038 |
< |
more visible with the 9\AA\ cutoff, where the image charges exert a |
| 2039 |
< |
greater force on the central particle. The error bars for the {\sc sf} |
| 2040 |
< |
methods show plus or minus the standard deviation of the density |
| 2041 |
< |
measurement at each temperature.} |
| 2042 |
< |
\label{fig:t5peDensities} |
| 2043 |
< |
\end{figure} |
| 2044 |
< |
|
| 2045 |
< |
Figure \ref{fig:t5peDensities} shows the densities calculated for |
| 2046 |
< |
TIP5P-E using differing electrostatic corrections overlaid on the |
| 2047 |
< |
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 2048 |
< |
technique are close to, though typically lower than, those calculated |
| 2049 |
< |
while using the Ewald summation. These slightly reduced densities |
| 2050 |
< |
indicate that the pressure component from the image charges at |
| 2051 |
< |
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 2052 |
< |
portion of the Ewald summation. Bringing the image charges closer to |
| 2053 |
< |
the central particle by choosing a 9\AA\ R$_\textrm{c}$ (rather than |
| 2054 |
< |
the preferred 12\AA\ R$_\textrm{c}$) increases the strength of their |
| 2055 |
< |
interactions, resulting in a further reduction of the densities. |
| 2056 |
< |
|
| 2057 |
< |
Because the strength of the image charge interactions has a noticable |
| 2058 |
< |
effect on the density, we would expect the use of electrostatic |
| 2059 |
< |
damping to also play a role in these calculations. Larger values of |
| 2060 |
< |
$\alpha$ weaken the pair-interactions; and since electrostatic damping |
| 2061 |
< |
is distance-dependent, force components from the image charges will be |
| 2062 |
< |
reduced more than those from particles close the the central |
| 2063 |
< |
charge. This effect is visible in figure \ref{fig:t5peDensities} with |
| 2064 |
< |
the damped {\sc sf} sums showing slightly higher densities; however, |
| 2065 |
< |
it is apparent that the choice of cutoff radius plays a much more |
| 2066 |
< |
important role in the resulting densities. |
| 2067 |
< |
|
| 2068 |
< |
As a final note, all of the above density calculations were performed |
| 2069 |
< |
with systems of 512 water molecules. Rick observed a system sized |
| 2070 |
< |
dependence of the computed densities when using the Ewald summation, |
| 2071 |
< |
most likely due to his tying of the convergence parameter to the box |
| 2072 |
< |
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
| 2073 |
< |
calculated densities were on average 0.002 to 0.003 g/cm$^3$ lower. A |
| 2074 |
< |
system size of 256 molecules would force the use of a shorter |
| 2075 |
< |
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also |
| 2076 |
< |
lower the densities. Moving to larger systems, as long as the |
| 2077 |
< |
R$_\textrm{c}$ remains at a fixed value, we would expect the densities |
| 2078 |
< |
to remain constant. |
| 2079 |
< |
|
| 2080 |
< |
\subsection{Liquid Structure}\label{sec:t5peLiqStructure} |
| 2081 |
< |
|
| 2082 |
< |
A common function considered when developing and comparing water |
| 2083 |
< |
models is the oxygen-oxygen radial distribution function |
| 2084 |
< |
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of |
| 2085 |
< |
finding a pair of oxygen atoms some distance ($r$) apart relative to a |
| 2086 |
< |
random distribution at the same density.\cite{Allen87} It is |
| 2087 |
< |
calculated via |
| 2088 |
< |
\begin{equation} |
| 2089 |
< |
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i} |
| 2090 |
< |
\delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle, |
| 2091 |
< |
\label{eq:GOOofR} |
| 2092 |
< |
\end{equation} |
| 2093 |
< |
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen |
| 2094 |
< |
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or |
| 2095 |
< |
neutron scattering experiments through the oxygen-oxygen structure |
| 2096 |
< |
factor ($S_\textrm{OO}(k)$) by the following relationship: |
| 2097 |
< |
\begin{equation} |
| 2098 |
< |
S_\textrm{OO}(k) = 1 + 4\pi\rho |
| 2099 |
< |
\int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r. |
| 2100 |
< |
\label{eq:SOOofK} |
| 2101 |
< |
\end{equation} |
| 2102 |
< |
Thus, $S_\textrm{OO}(k)$ is related to the Fourier-Laplace transform |
| 2103 |
< |
of $g_\textrm{OO}(r)$. |
| 2104 |
< |
|
| 2105 |
< |
The expermentally determined $g_\textrm{OO}(r)$ for liquid water has |
| 2106 |
< |
been compared in great detail with the various common water models, |
| 2107 |
< |
and TIP5P was found to be in better agreement than other rigid, |
| 2108 |
< |
non-polarizable models.\cite{Sorenson00} This excellent agreement with |
| 2109 |
< |
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To |
| 2110 |
< |
check whether the choice of using the Ewald summation or the {\sc sf} |
| 2111 |
< |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K |
| 2112 |
< |
and 1atm were determined for the systems compared in the previous |
| 2113 |
< |
section. |
| 2114 |
< |
|
| 2115 |
< |
\begin{figure} |
| 2116 |
< |
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 2117 |
< |
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298K and |
| 2118 |
< |
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc |
| 2119 |
< |
sf} technique with varying parameters. Even with the reduced densities |
| 2120 |
< |
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially |
| 2121 |
< |
identical.} |
| 2122 |
< |
\label{fig:t5peGofRs} |
| 2123 |
< |
\end{figure} |
| 2124 |
< |
|
| 2125 |
< |
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc |
| 2126 |
< |
sf} technique with a various parameters are overlaid on the |
| 2127 |
< |
$g_\textrm{OO}(r)$ while using the Ewald summation. The differences in |
| 2128 |
< |
density do not appear to have any effect on the liquid structure as |
| 2129 |
< |
the $g_\textrm{OO}(r)$s are indistinquishable. These results indicate |
| 2130 |
< |
that the $g_\textrm{OO}(r)$ is insensitive to the choice of |
| 2131 |
< |
electrostatic correction. |
| 2132 |
< |
|
| 2133 |
< |
\subsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 2134 |
< |
|
| 2135 |
< |
In addition to the density, there are a variety of thermodynamic |
| 2136 |
< |
quantities that can be calculated for water and compared directly to |
| 2137 |
< |
experimental values. Some of these additional quatities include the |
| 2138 |
< |
latent heat of vaporization ($\Delta H_\textrm{vap}$), the constant |
| 2139 |
< |
pressure heat capacity ($C_p$), the isothermal compressibility |
| 2140 |
< |
($\kappa_T$), the thermal expansivity ($\alpha_p$), and the static |
| 2141 |
< |
dielectric constant ($\epsilon$). All of these properties were |
| 2142 |
< |
calculated for TIP5P-E with the Ewald summation, so they provide a |
| 2143 |
< |
good set for comparisons involving the {\sc sf} technique. |
| 2144 |
< |
|
| 2145 |
< |
The $\Delta H_\textrm{vap}$ is the enthalpy change required to |
| 2146 |
< |
transform one mol of substance from the liquid phase to the gas |
| 2147 |
< |
phase.\cite{Berry00} In molecular simulations, this quantity can be |
| 2148 |
< |
determined via |
| 2149 |
< |
\begin{equation} |
| 2150 |
< |
\begin{split} |
| 2151 |
< |
\Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq.} \\ |
| 2152 |
< |
&= E_\textrm{gas} - E_\textrm{liq.} |
| 2153 |
< |
+ p(V_\textrm{gas} - V_\textrm{liq.}) \\ |
| 2154 |
< |
&\approx -\frac{\langle U_\textrm{liq.}\rangle}{N}+ RT, |
| 2155 |
< |
\end{split} |
| 2156 |
< |
\label{eq:DeltaHVap} |
| 2157 |
< |
\end{equation} |
| 2158 |
< |
where $E$ is the total energy, $U$ is the potential energy, $p$ is the |
| 2159 |
< |
pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is |
| 2160 |
< |
the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As |
| 2161 |
< |
seen in the last line of equation \ref{eq:DeltaHVap}, we can |
| 2162 |
< |
approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas |
| 2163 |
< |
state. This allows us to cancel the kinetic energy terms, leaving only |
| 2164 |
< |
the $U_\textrm{liq.}$ term. Additionally, the $pV$ term for the gas is |
| 2165 |
< |
several orders of magnitude larger than that of the liquid, so we can |
| 2166 |
< |
neglect the liquid $pV$ term. |
| 2167 |
< |
|
| 2168 |
< |
The remaining thermodynamic properties can all be calculated from |
| 2169 |
< |
fluctuations of the enthalpy, volume, and system dipole |
| 2170 |
< |
moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the |
| 2171 |
< |
enthalpy in constant pressure simulations via |
| 2172 |
< |
\begin{equation} |
| 2173 |
< |
\begin{split} |
| 2174 |
< |
C_p = \left(\frac{\partial H}{\partial T}\right)_{N,p} |
| 2175 |
< |
= \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2}, |
| 2176 |
< |
\end{split} |
| 2177 |
< |
\label{eq:Cp} |
| 2178 |
< |
\end{equation} |
| 2179 |
< |
where $k_B$ is Boltzmann's constant. $\kappa_T$ can be calculated via |
| 2180 |
< |
\begin{equation} |
| 2181 |
< |
\begin{split} |
| 2182 |
< |
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T} |
| 2183 |
< |
= \frac{(\langle V^2\rangle_{N,P,T} - \langle V\rangle^2_{N,P,T})} |
| 2184 |
< |
{k_BT\langle V\rangle_{N,P,T}}, |
| 2185 |
< |
\end{split} |
| 2186 |
< |
\label{eq:kappa} |
| 2187 |
< |
\end{equation} |
| 2188 |
< |
and $\alpha_p$ can be calculated via |
| 2189 |
< |
\begin{equation} |
| 2190 |
< |
\begin{split} |
| 2191 |
< |
\alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P} |
| 2192 |
< |
= \frac{(\langle VH\rangle_{N,P,T} |
| 2193 |
< |
- \langle V\rangle_{N,P,T}\langle H\rangle_{N,P,T})} |
| 2194 |
< |
{k_BT^2\langle V\rangle_{N,P,T}}. |
| 2195 |
< |
\end{split} |
| 2196 |
< |
\label{eq:alpha} |
| 2197 |
< |
\end{equation} |
| 2198 |
< |
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can |
| 2199 |
< |
be calculated for systems of non-polarizable substances via |
| 2200 |
< |
\begin{equation} |
| 2201 |
< |
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
| 2202 |
< |
\label{eq:staticDielectric} |
| 2203 |
< |
\end{equation} |
| 2204 |
< |
where $\epsilon_0$ is the permittivity of free space and $\langle |
| 2205 |
< |
M^2\rangle$ is the fluctuation of the system dipole |
| 2206 |
< |
moment.\cite{Allen87} The numerator in the fractional term in equation |
| 2207 |
< |
\ref{eq:staticDielectric} is the fluctuation of the simulation-box |
| 2208 |
< |
dipole moment, identical to the quantity calculated in the |
| 2209 |
< |
finite-system Kirkwood $g$ factor ($G_k$): |
| 2210 |
< |
\begin{equation} |
| 2211 |
< |
G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
| 2212 |
< |
\label{eq:KirkwoodFactor} |
| 2213 |
< |
\end{equation} |
| 2214 |
< |
where $\mu$ is the dipole moment of a single molecule of the |
| 2215 |
< |
homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The |
| 2216 |
< |
fluctuation term in both equation \ref{eq:staticDielectric} and |
| 2217 |
< |
\ref{eq:KirkwoodFactor} is calculated as follows, |
| 2218 |
< |
\begin{equation} |
| 2219 |
< |
\begin{split} |
| 2220 |
< |
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
| 2221 |
< |
- \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
| 2222 |
< |
&= \langle M_x^2+M_y^2+M_z^2\rangle |
| 2223 |
< |
- (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
| 2224 |
< |
+ \langle M_x\rangle^2). |
| 2225 |
< |
\end{split} |
| 2226 |
< |
\label{eq:fluctBoxDipole} |
| 2227 |
< |
\end{equation} |
| 2228 |
< |
This fluctuation term can be accumulated during the simulation; |
| 2229 |
< |
however, it converges rather slowly, thus requiring multi-nanosecond |
| 2230 |
< |
simulation times.\cite{Horn04} In the case of tin-foil boundary |
| 2231 |
< |
conditions, the dielectric/surface term of equation \ref{eq:EwaldSum} |
| 2232 |
< |
is equal to zero. Since the {\sc sf} method also lacks this |
| 2233 |
< |
dielectric/surface term, equation \ref{eq:staticDielectric} is still |
| 2234 |
< |
valid for determining static dielectric constants. |
| 2235 |
< |
|
| 2236 |
< |
All of the above properties were calculated from the same trajectories |
| 2237 |
< |
used to determine the densities in section \ref{sec:t5peDensity} |
| 2238 |
< |
except for the static dielectric constants. The $\epsilon$ values were |
| 2239 |
< |
accumulated from 2ns $NVE$ ensemble trajectories with system densities |
| 2240 |
< |
fixed at the average values from the $NPT$ simulations at each of the |
| 2241 |
< |
temperatures. The resulting values are displayed in figure |
| 2242 |
< |
\ref{fig:t5peThermo}. |
| 2243 |
< |
\begin{figure} |
| 2244 |
< |
\centering |
| 2245 |
< |
\includegraphics[width=5.5in]{./figures/t5peThermo.pdf} |
| 2246 |
< |
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation |
| 2247 |
< |
and the {\sc sf} techniques along with the experimental values. Units |
| 2248 |
< |
for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$, |
| 2249 |
< |
cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$, |
| 2250 |
< |
and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from |
| 2251 |
< |
reference \cite{Rick04}. Experimental values for $\Delta |
| 2252 |
< |
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference |
| 2253 |
< |
\cite{Kell75}. Experimental values for $C_p$ are from reference |
| 2254 |
< |
\cite{Wagner02}. Experimental values for $\epsilon$ are from reference |
| 2255 |
< |
\cite{Malmberg56}.} |
| 2256 |
< |
\label{fig:t5peThermo} |
| 2257 |
< |
\end{figure} |
| 2258 |
< |
|
| 2259 |
< |
As observed for the density in section \ref{sec:t5peDensity}, the |
| 2260 |
< |
property trends with temperature seen when using the Ewald summation |
| 2261 |
< |
are reproduced with the {\sc sf} technique. Differences include the |
| 2262 |
< |
calculated values of $\Delta H_\textrm{vap}$ underpredicting the Ewald |
| 2263 |
< |
values. This is to be expected due to the direct weakening of the |
| 2264 |
< |
electrostatic interaction through forced neutralization in {\sc |
| 2265 |
< |
sf}. This results in an increase of the intermolecular potential |
| 2266 |
< |
producing lower values from equation \ref{eq:DeltaHVap}. The slopes of |
| 2267 |
< |
these values with temperature are similar to that seen using the Ewald |
| 2268 |
< |
summation; however, they are both steeper than the experimental trend, |
| 2269 |
< |
indirectly resulting in the inflated $C_p$ values at all temperatures. |
| 2270 |
< |
|
| 2271 |
< |
Above the supercooled regim\'{e}, $C_p$, $\kappa_T$, and $\alpha_p$ |
| 2272 |
< |
values all overlap within error. As indicated for the $\Delta |
| 2273 |
< |
H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph, |
| 2274 |
< |
the deviations between experiment and simulation in this region are |
| 2275 |
< |
not the fault of the electrostatic summation methods but are due to |
| 2276 |
< |
the TIP5P class model itself. Like most rigid, non-polarizable, |
| 2277 |
< |
point-charge water models, the density decreases with temperature at a |
| 2278 |
< |
much faster rate than experiment (see figure |
| 2279 |
< |
\ref{fig:t5peDensities}). The reduced density leads to the inflated |
| 2280 |
< |
compressibility and expansivity values at higher temperatures seen |
| 2281 |
< |
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability |
| 2282 |
< |
and many-body effects are required in order for simulation to overcome |
| 2283 |
< |
these differences with experiment.\cite{Laasonen93,Donchev06} |
| 2284 |
< |
|
| 2285 |
< |
At temperatures below the freezing point for experimental water, the |
| 2286 |
< |
differences between {\sc sf} and the Ewald summation results are more |
| 2287 |
< |
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region |
| 2288 |
< |
indicate a more pronounced transition in the supercooled regim\'{e}, |
| 2289 |
< |
particularly in the case of {\sc sf} without damping. This points to |
| 2290 |
< |
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 2291 |
< |
temperatures below 250K in these simulations. Because the systems are |
| 2292 |
< |
locked in different regions of phase-space, comparisons between |
| 2293 |
< |
properties at these temperatures are not exactly fair. This |
| 2294 |
< |
observation is explored in more detail in section |
| 2295 |
< |
\ref{sec:t5peDynamics}. |
| 2296 |
< |
|
| 2297 |
< |
The final thermodynamic property displayed in figure |
| 2298 |
< |
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
| 2299 |
< |
between the Ewald summation and the {\sc sf} technique (and experiment |
| 2300 |
< |
for that matter). It is known that the dielectric constant is |
| 2301 |
< |
dependent upon and quite sensitive to the imposed boundary |
| 2302 |
< |
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the |
| 2303 |
< |
converged $\epsilon$ values accumulated for the {\sc sf} |
| 2304 |
< |
simulations. Lack of a damping function results in dielectric |
| 2305 |
< |
constants significantly smaller than that obtained using the Ewald |
| 2306 |
< |
sum. Increasing the damping coefficient to 0.2\AA$^{-1}$ improves the |
| 2307 |
< |
agreement considerably. It should be noted that the choice of the |
| 2308 |
< |
``Ewald coefficient'' value also has a significant effect on the |
| 2309 |
< |
calculated value when using the Ewald summation. In the simulations of |
| 2310 |
< |
TIP5P-E with the Ewald sum, this screening parameter was tethered to |
| 2311 |
< |
the simulation box size (as was the $R_\textrm{c}$).\cite{Rick04} |
| 2312 |
< |
Systems with larger screening parameters reported larger dielectric |
| 2313 |
< |
constant values, the same behavior we see here with {\sc sf}. In |
| 2314 |
< |
section \ref{sec:dampingDielectric}, this connection is further |
| 2315 |
< |
explored as optimal damping coefficients are determined for {\sc |
| 2316 |
< |
sf} for capturing the dielectric behavior. |
| 2317 |
< |
|
| 2318 |
< |
\subsection{Dynamic Properties}\label{sec:t5peDynamics} |
| 2319 |
< |
|
| 2320 |
< |
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 2321 |
< |
method, 200ps $NVE$ simulations were performed for each temperature at |
| 2322 |
< |
the average density reported by the $NPT$ simulations. The |
| 2323 |
< |
self-diffusion constants ($D$) were calculated with the Einstein |
| 2324 |
< |
relation using the mean square displacement (MSD), |
| 2325 |
< |
\begin{equation} |
| 2326 |
< |
D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t}, |
| 2327 |
< |
\label{eq:MSD} |
| 2328 |
< |
\end{equation} |
| 2329 |
< |
where $t$ is time, and $\mathbf{r}_i$ is the position of particle |
| 2330 |
< |
$i$.\cite{Allen87} Figure \ref{fig:ExampleMSD} shows an example MSD |
| 2331 |
< |
plot. As labeled in the figure, MSD plots consist of three distinct |
| 2332 |
< |
regions: |
| 2333 |
< |
|
| 2334 |
< |
\begin{enumerate}[itemsep=0pt] |
| 2335 |
< |
\item parabolic short-time ballistic motion, |
| 2336 |
< |
\item linear diffusive regime, and |
| 2337 |
< |
\item poor statistic region at long-time. |
| 2338 |
< |
\end{enumerate} |
| 2339 |
< |
The slope from the linear region (region 2) is used to calculate $D$. |
| 2340 |
< |
\begin{figure} |
| 2341 |
< |
\centering |
| 2342 |
< |
\includegraphics[width=3.5in]{./figures/ExampleMSD.pdf} |
| 2343 |
< |
\caption{Example plot of mean square displacement verses time. The |
| 2344 |
< |
left red region is the ballistic motion regime, the middle green |
| 2345 |
< |
region is the linear diffusive regime, and the right blue region is |
| 2346 |
< |
the region with poor statistics.} |
| 2347 |
< |
\label{fig:ExampleMSD} |
| 2348 |
< |
\end{figure} |
| 2349 |
< |
|
| 2350 |
< |
\begin{figure} |
| 2351 |
< |
\centering |
| 2352 |
< |
\includegraphics[width=3.5in]{./figures/waterFrame.pdf} |
| 2353 |
< |
\caption{Body-fixed coordinate frame for a water molecule. The |
| 2354 |
< |
respective molecular principle axes point in the direction of the |
| 2355 |
< |
labeled frame axes.} |
| 2356 |
< |
\label{fig:waterFrame} |
| 2357 |
< |
\end{figure} |
| 2358 |
< |
In addition to translational diffusion, reorientational time constants |
| 2359 |
< |
were calculated for comparisons with the Ewald simulations and with |
| 2360 |
< |
experiments. These values were determined from 25ps $NVE$ trajectories |
| 2361 |
< |
through calculation of the orientational time correlation function, |
| 2362 |
< |
\begin{equation} |
| 2363 |
< |
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t) |
| 2364 |
< |
\cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle, |
| 2365 |
< |
\label{eq:OrientCorr} |
| 2366 |
< |
\end{equation} |
| 2367 |
< |
where $P_l$ is the Legendre polynomial of order $l$ and |
| 2368 |
< |
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along |
| 2369 |
< |
principle axis $\alpha$. The principle axis frame for these water |
| 2370 |
< |
molecules is shown in figure \ref{fig:waterFrame}. As an example, |
| 2371 |
< |
$C_l^y$ is calculated from the time evolution of the unit vector |
| 2372 |
< |
connecting the two hydrogen atoms. |
| 2373 |
< |
|
| 2374 |
< |
\begin{figure} |
| 2375 |
< |
\centering |
| 2376 |
< |
\includegraphics[width=3.5in]{./figures/ExampleOrientCorr.pdf} |
| 2377 |
< |
\caption{Example plots of the orientational autocorrelation functions |
| 2378 |
< |
for the first and second Legendre polynomials. These curves show the |
| 2379 |
< |
time decay of the unit vector along the $y$ principle axis.} |
| 2380 |
< |
\label{fig:OrientCorr} |
| 2381 |
< |
\end{figure} |
| 2382 |
< |
From the orientation autocorrelation functions, we can obtain time |
| 2383 |
< |
constants for rotational relaxation. Figure \ref{fig:OrientCorr} shows |
| 2384 |
< |
some example plots of orientational autocorrelation functions for the |
| 2385 |
< |
first and second Legendre polynomials. The relatively short time |
| 2386 |
< |
portions (between 1 and 3ps for water) of these curves can be fit to |
| 2387 |
< |
an exponential decay to obtain these constants, and they are directly |
| 2388 |
< |
comparable to water orientational relaxation times from nuclear |
| 2389 |
< |
magnetic resonance (NMR). The relaxation constant obtained from |
| 2390 |
< |
$C_2^y(t)$ is of particular interest because it is about the principle |
| 2391 |
< |
axis with the minimum moment of inertia and should thereby dominate |
| 2392 |
< |
the orientational relaxation of the molecule.\cite{Impey82} This means |
| 2393 |
< |
that $C_2^y(t)$ should provide the best comparison to the NMR |
| 2394 |
< |
relaxation data. |
| 2395 |
< |
|
| 2396 |
< |
\begin{figure} |
| 2397 |
< |
\centering |
| 2398 |
< |
\includegraphics[width=5.5in]{./figures/t5peDynamics.pdf} |
| 2399 |
< |
\caption{Diffusion constants ({\it upper}) and reorientational time |
| 2400 |
< |
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf} |
| 2401 |
< |
technique compared with experiment. Data at temperatures less that |
| 2402 |
< |
0$^\circ$C were not included in the $\tau_2^y$ plot to allow for |
| 2403 |
< |
easier comparisons in the more relavent temperature regime.} |
| 2404 |
< |
\label{fig:t5peDynamics} |
| 2405 |
< |
\end{figure} |
| 2406 |
< |
Results for the diffusion constants and reorientational time constants |
| 2407 |
< |
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is |
| 2408 |
< |
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using |
| 2409 |
< |
the Ewald sum are reproduced with the {\sc sf} techinque. The enhanced |
| 2410 |
< |
diffusion at high temperatures are again the product of the lower |
| 2411 |
< |
densities in comparison with experiment and do not provide any special |
| 2412 |
< |
insight into differences between the electrostatic summation |
| 2413 |
< |
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to |
| 2414 |
< |
diffuse a little faster than with the Ewald sum; however, use of light |
| 2415 |
< |
to moderate damping results in indistiguishable $D$ values. Though not |
| 2416 |
< |
apparent in this figure, {\sc sf} values at the lowest temperature are |
| 2417 |
< |
approximately an order of magnitude lower than with Ewald. These |
| 2418 |
< |
values support the observation from section \ref{sec:t5peThermo} that |
| 2419 |
< |
there appeared to be a change to a more glassy-like phase with the |
| 2420 |
< |
{\sc sf} technique at these lower temperatures. |
| 2421 |
< |
|
| 2422 |
< |
The $\tau_2^y$ results in the lower frame of figure |
| 2423 |
< |
\ref{fig:t5peDynamics} show a much greater difference between the {\sc |
| 2424 |
< |
sf} results and the Ewald results. At all temperatures shown, TIP5P-E |
| 2425 |
< |
relaxes faster than experiment with the Ewald sum while tracking |
| 2426 |
< |
experiment fairly well when using the {\sc sf} technique, independent |
| 2427 |
< |
of the choice of damping constant. Their are several possible reasons |
| 2428 |
< |
for this deviation between techniques. The Ewald results were taken |
| 2429 |
< |
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick |
| 2430 |
< |
calculation from a 10ps trajectory with {\sc sf} with an $\alpha$ of |
| 2431 |
< |
0.2\AA$^-1$ at 25$^\circ$C showed a 0.4ps drop in $\tau_2^y$, placing |
| 2432 |
< |
the result more in line with that obtained using the Ewald sum. These |
| 2433 |
< |
results support this explanation; however, recomputing the results to |
| 2434 |
< |
meet a poorer statistical standard is counter-productive. Assuming the |
| 2435 |
< |
Ewald results are not the product of poor statistics, differences in |
| 2436 |
< |
techniques to integrate the orientational motion could also play a |
| 2437 |
< |
role. {\sc shake} is the most commonly used technique for |
| 2438 |
< |
approximating rigid-body orientational motion,\cite{Ryckaert77} where |
| 2439 |
< |
as in {\sc oopse}, we maintain and integrate the entire rotation |
| 2440 |
< |
matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake} |
| 2441 |
< |
is an iterative constraint technique, if the convergence tolerances |
| 2442 |
< |
are raised for increased performance, error will accumulate in the |
| 2443 |
< |
orientational motion. Finally, the Ewald results were calculated using |
| 2444 |
< |
the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf} |
| 2445 |
< |
calculations. The additional mode of motion due to the thermostat will |
| 2446 |
< |
alter the dynamics, resulting in differences between $NVT$ and $NVE$ |
| 2447 |
< |
results. These differences are increasingly noticable as the |
| 2448 |
< |
thermostat time constant decreases. |
| 2449 |
< |
|
| 2450 |
< |
\section{Damping of Point Multipoles}\label{sec:dampingMultipoles} |
| 2451 |
< |
|
| 2452 |
< |
|
| 2453 |
< |
|
| 2454 |
< |
\section{Damping and Dielectric Constants}\label{sec:dampingDielectric} |
| 2455 |
< |
|
| 2456 |
< |
\section{Conclusions}\label{sec:PairwiseConclusions} |
| 2457 |
< |
|
| 2458 |
< |
The above investigation of pairwise electrostatic summation techniques |
| 2459 |
< |
shows that there are viable and computationally efficient alternatives |
| 2460 |
< |
to the Ewald summation. These methods are derived from the damped and |
| 2461 |
< |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 2462 |
< |
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
| 2463 |
< |
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
| 2464 |
< |
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 2465 |
< |
energetic and dynamic characteristics exhibited by simulations |
| 2466 |
< |
employing lattice summation techniques. The cumulative energy |
| 2467 |
< |
difference results showed the undamped {\sc sf} and moderately damped |
| 2468 |
< |
{\sc sp} methods produced results nearly identical to {\sc spme}. |
| 2469 |
< |
Similarly for the dynamic features, the undamped or moderately damped |
| 2470 |
< |
{\sc sf} and moderately damped {\sc sp} methods produce force and |
| 2471 |
< |
torque vector magnitude and directions very similar to the expected |
| 2472 |
< |
values. These results translate into long-time dynamic behavior |
| 2473 |
< |
equivalent to that produced in simulations using {\sc spme}. |
| 2474 |
< |
|
| 2475 |
< |
As in all purely-pairwise cutoff methods, these methods are expected |
| 2476 |
< |
to scale approximately {\it linearly} with system size, and they are |
| 2477 |
< |
easily parallelizable. This should result in substantial reductions |
| 2478 |
< |
in the computational cost of performing large simulations. |
| 2479 |
< |
|
| 2480 |
< |
Aside from the computational cost benefit, these techniques have |
| 2481 |
< |
applicability in situations where the use of the Ewald sum can prove |
| 2482 |
< |
problematic. Of greatest interest is their potential use in |
| 2483 |
< |
interfacial systems, where the unmodified lattice sum techniques |
| 2484 |
< |
artificially accentuate the periodicity of the system in an |
| 2485 |
< |
undesirable manner. There have been alterations to the standard Ewald |
| 2486 |
< |
techniques, via corrections and reformulations, to compensate for |
| 2487 |
< |
these systems; but the pairwise techniques discussed here require no |
| 2488 |
< |
modifications, making them natural tools to tackle these problems. |
| 2489 |
< |
Additionally, this transferability gives them benefits over other |
| 2490 |
< |
pairwise methods, like reaction field, because estimations of physical |
| 2491 |
< |
properties (e.g. the dielectric constant) are unnecessary. |
| 2492 |
< |
|
| 2493 |
< |
If a researcher is using Monte Carlo simulations of large chemical |
| 2494 |
< |
systems containing point charges, most structural features will be |
| 2495 |
< |
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 2496 |
< |
method with an electrostatic damping of 0.2\AA$^{-1}$. These methods |
| 2497 |
< |
would also be appropriate for molecular dynamics simulations where the |
| 2498 |
< |
data of interest is either structural or short-time dynamical |
| 2499 |
< |
quantities. For long-time dynamics and collective motions, the safest |
| 2500 |
< |
pairwise method we have evaluated is the {\sc sf} method with an |
| 2501 |
< |
electrostatic damping between 0.2 and 0.25\AA$^{-1}$. |
| 2502 |
< |
|
| 2503 |
< |
We are not suggesting that there is any flaw with the Ewald sum; in |
| 2504 |
< |
fact, it is the standard by which these simple pairwise sums have been |
| 2505 |
< |
judged. However, these results do suggest that in the typical |
| 2506 |
< |
simulations performed today, the Ewald summation may no longer be |
| 2507 |
< |
required to obtain the level of accuracy most researchers have come to |
| 2508 |
< |
expect. |
| 46 |
> |
\input{electrostaticsChapter} |
| 47 |
|
|
| 48 |
|
\chapter{\label{chap:water}SIMPLE MODELS FOR WATER} |
| 49 |
|
|
