| 1 |
< |
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION |
| 2 |
< |
TECHNIQUES} |
| 1 |
> |
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES} |
| 2 |
|
|
| 3 |
< |
In molecular simulations, proper accumulation of the electrostatic |
| 4 |
< |
interactions is essential and is one of the most |
| 3 |
> |
In molecular simulations, proper accumulation of electrostatic |
| 4 |
> |
interactions is essential and one of the most |
| 5 |
|
computationally-demanding tasks. The common molecular mechanics force |
| 6 |
|
fields represent atomic sites with full or partial charges protected |
| 7 |
< |
by repulsive Lennard-Jones interactions. This means that nearly |
| 8 |
< |
every pair interaction involves a calculation of charge-charge forces. |
| 7 |
> |
by repulsive Lennard-Jones interactions. This means that nearly every |
| 8 |
> |
pair interaction involves a calculation of charge-charge forces. |
| 9 |
|
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
| 10 |
|
interactions quickly become the most expensive part of molecular |
| 11 |
|
simulations. Historically, the electrostatic pair interaction would |
| 32 |
|
|
| 33 |
|
In this chapter, we focus on a new set of pairwise methods devised by |
| 34 |
|
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 35 |
< |
methods along with a few other mixed methods (i.e. reaction field) are |
| 36 |
< |
compared with the smooth particle mesh Ewald |
| 35 |
> |
methods, along with a few other mixed methods (i.e. reaction field), |
| 36 |
> |
are compared with the smooth particle mesh Ewald |
| 37 |
|
sum,\cite{Onsager36,Essmann99} which is our reference method for |
| 38 |
|
handling long-range electrostatic interactions. The new methods for |
| 39 |
|
handling electrostatics have the potential to scale linearly with |
| 40 |
< |
increasing system size since they involve only a simple modification |
| 40 |
> |
increasing system size, since they involve only a simple modification |
| 41 |
|
to the direct pairwise sum. They also lack the added periodicity of |
| 42 |
|
the Ewald sum, so they can be used for systems which are non-periodic |
| 43 |
|
or which have one- or two-dimensional periodicity. Below, these |
| 44 |
< |
methods are evaluated using a variety of model systems to |
| 45 |
< |
establish their usability in molecular simulations. |
| 44 |
> |
methods are evaluated using a variety of model systems to establish |
| 45 |
> |
their usability in molecular simulations. |
| 46 |
|
|
| 47 |
|
\section{The Ewald Sum} |
| 48 |
|
|
| 92 |
|
where $\alpha$ is the damping or convergence parameter with units of |
| 93 |
|
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 94 |
|
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 95 |
< |
constant of the surrounding medium. The final two terms of |
| 96 |
< |
equation (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 97 |
< |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 98 |
< |
dipolar term was neglected in early applications in molecular |
| 99 |
< |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 100 |
< |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 101 |
< |
dipole moment which is magnified through replication of the periodic |
| 102 |
< |
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
| 103 |
< |
system is said to be using conducting (or ``tin-foil'') boundary |
| 104 |
< |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 106 |
< |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 107 |
< |
time. Initially, due to the small system sizes that could be |
| 108 |
< |
simulated feasibly, the entire simulation box was replicated to |
| 109 |
< |
convergence. In more modern simulations, the systems have grown large |
| 110 |
< |
enough that a real-space cutoff could potentially give convergent |
| 111 |
< |
behavior. Indeed, it has been observed that with the choice of a |
| 112 |
< |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 113 |
< |
rapidly convergent and small relative to the real-space |
| 114 |
< |
portion.\cite{Karasawa89,Kolafa92} |
| 95 |
> |
constant of the surrounding medium. The final two terms of equation |
| 96 |
> |
(\ref{eq:EwaldSum}) are a particle-self term and a dipolar term for |
| 97 |
> |
interacting with a surrounding dielectric.\cite{Allen87} This dipolar |
| 98 |
> |
term was neglected in early applications of this technique in |
| 99 |
> |
molecular simulations,\cite{Brush66,Woodcock71} until it was |
| 100 |
> |
introduced by de Leeuw {\it et al.} to address situations where the |
| 101 |
> |
unit cell has a dipole moment which is magnified through replication |
| 102 |
> |
of the periodic images.\cite{deLeeuw80,Smith81} If this term is taken |
| 103 |
> |
to be zero, the system is said to be using conducting (or |
| 104 |
> |
``tin-foil'') boundary conditions, $\epsilon_{\rm S} = \infty$. |
| 105 |
|
|
| 106 |
|
\begin{figure} |
| 107 |
|
\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf} |
| 112 |
|
convergence for the larger systems of charges that are common today.} |
| 113 |
|
\label{fig:ewaldTime} |
| 114 |
|
\end{figure} |
| 115 |
+ |
Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied |
| 116 |
+ |
over time. Initially, due to the small system sizes that could be |
| 117 |
+ |
simulated feasibly, the entire simulation box was replicated to |
| 118 |
+ |
convergence. In more modern simulations, the systems have grown large |
| 119 |
+ |
enough that a real-space cutoff could potentially give convergent |
| 120 |
+ |
behavior. Indeed, it has been observed that with the choice of a |
| 121 |
+ |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 122 |
+ |
rapidly convergent and small relative to the real-space |
| 123 |
+ |
portion.\cite{Karasawa89,Kolafa92} |
| 124 |
|
|
| 125 |
< |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 125 |
> |
The original Ewald summation is an $\mathcal{O}(N^2)$ algorithm. The |
| 126 |
|
convergence parameter $(\alpha)$ plays an important role in balancing |
| 127 |
|
the computational cost between the direct and reciprocal-space |
| 128 |
|
portions of the summation. The choice of this value allows one to |
| 129 |
|
select whether the real-space or reciprocal space portion of the |
| 130 |
< |
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
| 131 |
< |
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 130 |
> |
summation is an $\mathcal{O}(N^2)$ calculation (with the other being |
| 131 |
> |
$\mathcal{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 132 |
|
$\alpha$ and thoughtful algorithm development, this cost can be |
| 133 |
< |
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 133 |
> |
reduced to $\mathcal{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 134 |
|
taken to reduce the cost of the Ewald summation even further is to set |
| 135 |
|
$\alpha$ such that the real-space interactions decay rapidly, allowing |
| 136 |
|
for a short spherical cutoff. Then the reciprocal space summation is |
| 139 |
|
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 140 |
|
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 141 |
|
methods, the cost of the reciprocal-space portion of the Ewald |
| 142 |
< |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
| 142 |
> |
summation is reduced from $\mathcal{O}(N^2)$ down to $\mathcal{O}(N |
| 143 |
|
\log N)$. |
| 144 |
|
|
| 145 |
|
These developments and optimizations have made the use of the Ewald |
| 155 |
|
bringing them more in line with the cost of the full 3-D summation. |
| 156 |
|
|
| 157 |
|
Several studies have recognized that the inherent periodicity in the |
| 158 |
< |
Ewald sum can also have an effect on three-dimensional |
| 159 |
< |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 160 |
< |
Solvated proteins are essentially kept at high concentration due to |
| 161 |
< |
the periodicity of the electrostatic summation method. In these |
| 162 |
< |
systems, the more compact folded states of a protein can be |
| 163 |
< |
artificially stabilized by the periodic replicas introduced by the |
| 164 |
< |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 158 |
> |
Ewald sum can have an effect not just on reduced dimensionality |
| 159 |
> |
system, but on three-dimensional systems as |
| 160 |
> |
well.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 161 |
> |
As an example, solvated proteins are essentially kept at high |
| 162 |
> |
concentration due to the periodicity of the electrostatic summation |
| 163 |
> |
method. In these systems, the more compact folded states of a protein |
| 164 |
> |
can be artificially stabilized by the periodic replicas introduced by |
| 165 |
> |
the Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 166 |
|
considering the use of the Ewald summation where the assumed |
| 167 |
< |
periodicity would introduce spurious effects in the system dynamics. |
| 167 |
> |
periodicity would introduce spurious effects. |
| 168 |
|
|
| 169 |
|
|
| 170 |
|
\section{The Wolf and Zahn Methods} |
| 175 |
|
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
| 176 |
|
observed that the electrostatic interaction is effectively |
| 177 |
|
short-ranged in condensed phase systems and that neutralization of the |
| 178 |
< |
charge contained within the cutoff radius is crucial for potential |
| 178 |
> |
charges contained within the cutoff radius is crucial for potential |
| 179 |
|
stability. They devised a pairwise summation method that ensures |
| 180 |
|
charge neutrality and gives results similar to those obtained with the |
| 181 |
|
Ewald summation. The resulting shifted Coulomb potential includes |
| 182 |
|
image-charges subtracted out through placement on the cutoff sphere |
| 183 |
|
and a distance-dependent damping function (identical to that seen in |
| 184 |
< |
the real-space portion of the Ewald sum) to aid convergence |
| 184 |
> |
the real-space portion of the Ewald sum) to aid convergence: |
| 185 |
|
\begin{equation} |
| 186 |
|
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}} |
| 187 |
|
- \lim_{r_{ij}\rightarrow R_\textrm{c}} |
| 241 |
|
\label{eq:ZahnPot} |
| 242 |
|
\end{equation} |
| 243 |
|
and showed that this potential does fairly well at capturing the |
| 244 |
< |
structural and dynamic properties of water compared the same |
| 244 |
> |
structural and dynamic properties of water compared with the same |
| 245 |
|
properties obtained using the Ewald sum. |
| 246 |
|
|
| 247 |
|
\section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation} |
| 312 |
|
v(r) = \frac{q_i q_j}{r}, |
| 313 |
|
\label{eq:Coulomb} |
| 314 |
|
\end{equation} |
| 315 |
< |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 316 |
< |
al.}'s undamped prescription: |
| 315 |
> |
then the {\sc sp} form will give Wolf {\it et al.}'s undamped |
| 316 |
> |
prescription: |
| 317 |
|
\begin{equation} |
| 318 |
|
V_\textrm{SP}(r) = |
| 319 |
|
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 335 |
|
forces at the cutoff radius which results in energy drift during MD |
| 336 |
|
simulations. |
| 337 |
|
|
| 338 |
< |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 339 |
< |
will give, |
| 338 |
> |
The {\sc sf} form using the normal Coulomb potential will give, |
| 339 |
|
\begin{equation} |
| 340 |
|
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}} |
| 341 |
|
+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] |
| 349 |
|
\label{eq:SFForces} |
| 350 |
|
\end{equation} |
| 351 |
|
This formulation has the benefits that there are no discontinuities at |
| 352 |
< |
the cutoff radius, while the neutralizing image charges are present in |
| 352 |
> |
the cutoff radius and the neutralizing image charges are present in |
| 353 |
|
both the energy and force expressions. It would be simple to add the |
| 354 |
|
self-neutralizing term back when computing the total energy of the |
| 355 |
|
system, thereby maintaining the agreement with the Madelung energies. |
| 361 |
|
Wolf \textit{et al.} originally discussed the energetics of the |
| 362 |
|
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
| 363 |
|
insufficient for accurate determination of the energy with reasonable |
| 364 |
< |
cutoff distances. The calculated Madelung energies fluctuated around |
| 365 |
< |
the expected value as the cutoff radius was increased, but the |
| 364 |
> |
cutoff distances. The calculated Madelung energies fluctuated wildly |
| 365 |
> |
around the expected value, but as the cutoff radius was increased, the |
| 366 |
|
oscillations converged toward the correct value.\cite{Wolf99} A |
| 367 |
< |
damping function was incorporated to accelerate the convergence; and |
| 367 |
> |
damping function was incorporated to accelerate this convergence; and |
| 368 |
|
though alternative forms for the damping function could be |
| 369 |
|
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 370 |
|
chosen to mirror the effective screening used in the Ewald summation. |
| 374 |
|
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 375 |
|
\label{eq:dampCoulomb} |
| 376 |
|
\end{equation} |
| 377 |
< |
the shifted potential (Eq. (\ref{eq:SPPot})) becomes |
| 377 |
> |
the {\sc sp} potential function (Eq. (\ref{eq:SPPot})) becomes |
| 378 |
|
\begin{equation} |
| 379 |
|
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} |
| 380 |
|
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) |
| 389 |
|
\quad r\leqslant R_\textrm{c}. |
| 390 |
|
\label{eq:DSPForces} |
| 391 |
|
\end{equation} |
| 392 |
< |
Again, this damped shifted potential suffers from a |
| 393 |
< |
force-discontinuity at the cutoff radius, and the image charges play |
| 394 |
< |
no role in the forces. To remedy these concerns, one may derive a |
| 395 |
< |
{\sc sf} variant by including the derivative term in |
| 396 |
< |
equation (\ref{eq:shiftingForm}), |
| 392 |
> |
Again, this damped shifted potential suffers from a discontinuity in |
| 393 |
> |
the forces at the cutoff radius, and the image charges play no role in |
| 394 |
> |
the forces. To remedy these concerns, one may derive a {\sc sf} |
| 395 |
> |
variant by including the derivative term present in |
| 396 |
> |
equation~(\ref{eq:shiftingForm}), |
| 397 |
|
\begin{equation} |
| 398 |
|
\begin{split} |
| 399 |
|
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 423 |
|
\end{split} |
| 424 |
|
\end{equation} |
| 425 |
|
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
| 426 |
< |
equations (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
| 427 |
< |
recovered from equations (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
| 426 |
> |
equations (\ref{eq:SPPot}) through (\ref{eq:SFForces}) are correctly |
| 427 |
> |
recovered from equations (\ref{eq:DSPPot}) through (\ref{eq:DSFForces}). |
| 428 |
|
|
| 429 |
< |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 429 |
> |
This new {\sc sf} potential is similar to equation (\ref{eq:ZahnPot}) |
| 430 |
|
derived by Zahn \textit{et al.}; however, there are two important |
| 431 |
|
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation |
| 432 |
|
(\ref{eq:shiftingForm}) is equal to equation (\ref{eq:dampCoulomb}) |
| 436 |
|
portion is different. The missing $v_\textrm{c}$ term would not |
| 437 |
|
affect molecular dynamics simulations (although the computed energy |
| 438 |
|
would be expected to have sudden jumps as particle distances crossed |
| 439 |
< |
$R_c$). The sign problem is a potential source of errors, however. |
| 440 |
< |
In fact, it introduces a discontinuity in the forces at the cutoff, |
| 441 |
< |
because the force function is shifted in the wrong direction and |
| 442 |
< |
doesn't cross zero at $R_\textrm{c}$. |
| 439 |
> |
$R_c$); however, the sign problem is a potential source of errors. In |
| 440 |
> |
fact, equation~(\ref{eq:ZahnPot}) introduces a discontinuity in the |
| 441 |
> |
forces at the cutoff, because the force function is shifted in the |
| 442 |
> |
wrong direction and does not cross zero at $R_\textrm{c}$. |
| 443 |
|
|
| 444 |
|
Equations (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 445 |
|
electrostatic summation method in which the potential and forces are |
| 446 |
|
continuous at the cutoff radius and which incorporates the damping |
| 447 |
|
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
| 448 |
< |
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
| 449 |
< |
{\sc sf}, damping) are at reproducing the correct electrostatic |
| 448 |
> |
this chapter, we will evaluate exactly how good these methods ({\sc |
| 449 |
> |
sp}, {\sc sf}, damping) are at reproducing the correct electrostatic |
| 450 |
|
summation performed by the Ewald sum. |
| 451 |
|
|
| 452 |
|
|
| 453 |
|
\section{Evaluating Pairwise Summation Techniques} |
| 454 |
|
|
| 455 |
< |
In classical molecular mechanics simulations, there are two primary |
| 456 |
< |
techniques utilized to obtain information about the system of |
| 457 |
< |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 458 |
< |
techniques utilize pairwise summations of interactions between |
| 459 |
< |
particle sites, but they use these summations in different ways. |
| 455 |
> |
As mentioned in the introduction, there are two primary techniques |
| 456 |
> |
utilized to obtain information about the system of interest in |
| 457 |
> |
classical molecular mechanics simulations: Monte Carlo (MC) and |
| 458 |
> |
molecular dynamics (MD). Both of these techniques utilize pairwise |
| 459 |
> |
summations of interactions between particle sites, but they use these |
| 460 |
> |
summations in different ways. |
| 461 |
|
|
| 462 |
|
In MC, the potential energy difference between configurations dictates |
| 463 |
|
the progression of MC sampling. Going back to the origins of this |
| 476 |
|
electrostatic summation techniques, the dynamics in the short term |
| 477 |
|
will be indistinguishable. Because error in MD calculations is |
| 478 |
|
cumulative, one should expect greater deviation at longer times, |
| 479 |
< |
although methods which have large differences in the force and torque |
| 479 |
> |
and methods which have large differences in the force and torque |
| 480 |
|
vectors will diverge from each other more rapidly. |
| 481 |
|
|
| 482 |
|
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 483 |
|
|
| 484 |
+ |
\begin{figure} |
| 485 |
+ |
\centering |
| 486 |
+ |
\includegraphics[width = 3.5in]{./figures/dualLinear.pdf} |
| 487 |
+ |
\caption{Example least squares regressions of the configuration energy |
| 488 |
+ |
differences for SPC/E water systems. The upper plot shows a data set |
| 489 |
+ |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 490 |
+ |
shows a data set with a good correlation coefficient.} |
| 491 |
+ |
\label{fig:linearFit} |
| 492 |
+ |
\end{figure} |
| 493 |
|
The pairwise summation techniques (outlined in section |
| 494 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 495 |
|
studying the energy differences between conformations. We took the |
| 503 |
|
correlation (slope) and correlation coefficient for these regressions |
| 504 |
|
indicate perfect agreement between the alternative method and {\sc spme}. |
| 505 |
|
Sample correlation plots for two alternate methods are shown in |
| 506 |
< |
Fig. \ref{fig:linearFit}. |
| 506 |
> |
figure \ref{fig:linearFit}. |
| 507 |
|
|
| 499 |
– |
\begin{figure} |
| 500 |
– |
\centering |
| 501 |
– |
\includegraphics[width = 3.5in]{./figures/dualLinear.pdf} |
| 502 |
– |
\caption{Example least squares regressions of the configuration energy |
| 503 |
– |
differences for SPC/E water systems. The upper plot shows a data set |
| 504 |
– |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 505 |
– |
shows a data set with a good correlation coefficient.} |
| 506 |
– |
\label{fig:linearFit} |
| 507 |
– |
\end{figure} |
| 508 |
– |
|
| 508 |
|
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 509 |
|
were represented using 500 independent configurations. Thus, each of |
| 510 |
|
the alternative (non-Ewald) electrostatic summation methods was |
| 511 |
|
evaluated using an accumulated 873,250 configurational energy |
| 512 |
< |
differences. |
| 512 |
> |
differences. Results for and discussions regarding the individual |
| 513 |
> |
analysis of each of the system types appear in appendix |
| 514 |
> |
\ref{app:IndividualResults}, while the cumulative results over all the |
| 515 |
> |
investigated systems appear below in section~\ref{sec:EnergyResults}. |
| 516 |
|
|
| 515 |
– |
Results and discussion for the individual analysis of each of the |
| 516 |
– |
system types appear in sections \ref{sec:IndividualResults}, while the |
| 517 |
– |
cumulative results over all the investigated systems appear below in |
| 518 |
– |
sections \ref{sec:EnergyResults}. |
| 519 |
– |
|
| 517 |
|
\subsection{Molecular Dynamics and the Force and Torque |
| 518 |
|
Vectors}\label{sec:MDMethods} We evaluated the pairwise methods |
| 519 |
|
(outlined in section \ref{sec:ESMethods}) for use in MD simulations by |
| 525 |
|
comparing $\Delta E$ values. Instead of a single energy difference |
| 526 |
|
between two system configurations, we compared the magnitudes of the |
| 527 |
|
forces (and torques) on each molecule in each configuration. For a |
| 528 |
< |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 529 |
< |
vectors and 1000 torque vectors. With 500 configurations, this |
| 530 |
< |
results in 520,000 force and 500,000 torque vector comparisons. |
| 531 |
< |
Additionally, data from seven different system types was aggregated |
| 532 |
< |
before the comparison was made. |
| 528 |
> |
system of 1000 water molecules and 40 ions, there are 1040 force and |
| 529 |
> |
1000 torque vectors. With 500 configurations, this results in 520,000 |
| 530 |
> |
force and 500,000 torque vector comparisons. Additionally, data from |
| 531 |
> |
seven different system types was aggregated before comparisons were |
| 532 |
> |
made. |
| 533 |
|
|
| 534 |
|
The {\it directionality} of the force and torque vectors was |
| 535 |
|
investigated through measurement of the angle ($\theta$) formed |
| 544 |
|
unit sphere. Since this distribution is a measure of angular error |
| 545 |
|
between two different electrostatic summation methods, there is no |
| 546 |
|
{\it a priori} reason for the profile to adhere to any specific |
| 547 |
< |
shape. Thus, gaussian fits were used to measure the width of the |
| 547 |
> |
shape. Thus, Gaussian fits were used to measure the width of the |
| 548 |
|
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 549 |
|
each of these fits and was used to compare distribution widths. |
| 550 |
|
Values of $\sigma^2$ near zero indicate vector directions |
| 607 |
|
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 608 |
|
\item a 6~\AA\ radius sphere of Argon in water. |
| 609 |
|
\end{enumerate} |
| 613 |
– |
|
| 610 |
|
By utilizing the pairwise techniques (outlined in section |
| 611 |
|
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 612 |
|
charged particles, and mixtures of the two, we hope to discern under |
| 695 |
|
inverted triangles).} |
| 696 |
|
\label{fig:delE} |
| 697 |
|
\end{figure} |
| 702 |
– |
|
| 698 |
|
The most striking feature of this plot is how well the Shifted Force |
| 699 |
|
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 700 |
|
differences. For the undamped {\sc sf} method, and the |
| 709 |
|
some degree by using group based cutoffs with a switching |
| 710 |
|
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 711 |
|
significant improvement using the group-switched cutoff because the |
| 712 |
< |
salt and salt solution systems contain non-neutral groups. Section |
| 713 |
< |
\ref{sec:IndividualResults} includes results for systems comprised entirely |
| 714 |
< |
of neutral groups. |
| 712 |
> |
salt and salt solution systems contain non-neutral groups. Appendix |
| 713 |
> |
\ref{app:IndividualResults} includes results for systems comprised |
| 714 |
> |
entirely of neutral groups. |
| 715 |
|
|
| 716 |
|
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 717 |
|
the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$ |
| 718 |
|
shows an excellent correlation and quality of fit with the {\sc spme} |
| 719 |
< |
results, particularly with a cutoff radius greater than 12~\AA\. Use |
| 719 |
> |
results, particularly with a cutoff radius greater than 12~\AA . Use |
| 720 |
|
of a larger damping parameter is more helpful for the shortest cutoff |
| 721 |
|
shown, but it has a detrimental effect on simulations with larger |
| 722 |
|
cutoffs. |
| 730 |
|
|
| 731 |
|
The reaction field results illustrates some of that method's |
| 732 |
|
limitations, primarily that it was developed for use in homogeneous |
| 733 |
< |
systems; although it does provide results that are an improvement over |
| 734 |
< |
those from an unmodified cutoff. |
| 733 |
> |
systems. It does, however, provide results that are an improvement |
| 734 |
> |
over those from an unmodified cutoff. |
| 735 |
|
|
| 736 |
|
\section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults} |
| 737 |
|
|
| 754 |
|
inverted triangles).} |
| 755 |
|
\label{fig:frcMag} |
| 756 |
|
\end{figure} |
| 757 |
< |
|
| 758 |
< |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 759 |
< |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 760 |
< |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 761 |
< |
This is perhaps expected because this method explicitly incorporates a |
| 767 |
< |
smooth transition in the forces at the cutoff radius as well as the |
| 757 |
> |
Again, it is striking how well the {\sc sp} and {\sc sf} methods |
| 758 |
> |
reproduce the {\sc spme} forces. The undamped and weakly-damped {\sc |
| 759 |
> |
sf} method gives the best agreement with Ewald. This is perhaps |
| 760 |
> |
expected because this method explicitly incorporates a smooth |
| 761 |
> |
transition in the forces at the cutoff radius as well as the |
| 762 |
|
neutralizing image charges. |
| 763 |
|
|
| 764 |
|
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 769 |
|
|
| 770 |
|
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 771 |
|
method is generating usable forces. Further increases in damping, |
| 772 |
< |
while beneficial for simulations with a cutoff radius of 9~\AA\ , is |
| 772 |
> |
while beneficial for simulations with a cutoff radius of 9~\AA\ , are |
| 773 |
|
detrimental to simulations with larger cutoff radii. |
| 774 |
|
|
| 775 |
|
The reaction field results are surprisingly good, considering the poor |
| 776 |
|
quality of the fits for the $\Delta E$ results. There is still a |
| 777 |
< |
considerable degree of scatter in the data, but the forces correlate |
| 778 |
< |
well with the Ewald forces in general. We note that the reaction |
| 779 |
< |
field calculations do not include the pure NaCl systems, so these |
| 780 |
< |
results are partly biased towards conditions in which the method |
| 781 |
< |
performs more favorably. |
| 777 |
> |
considerable degree of scatter in the data, but in general, the forces |
| 778 |
> |
correlate well with the Ewald forces. We note that the pure NaCl |
| 779 |
> |
systems were not included in the system set used in the reaction field |
| 780 |
> |
calculations, so these results are partly biased towards conditions in |
| 781 |
> |
which the method performs more favorably. |
| 782 |
|
|
| 783 |
|
\begin{figure} |
| 784 |
|
\centering |
| 792 |
|
inverted triangles).} |
| 793 |
|
\label{fig:trqMag} |
| 794 |
|
\end{figure} |
| 801 |
– |
|
| 795 |
|
Molecular torques were only available from the systems which contained |
| 796 |
|
rigid molecules (i.e. the systems containing water). The data in |
| 797 |
< |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 797 |
> |
figure \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 798 |
|
|
| 799 |
< |
Torques appear to be much more sensitive to charges at a longer |
| 800 |
< |
distance. The striking feature in comparing the new electrostatic |
| 801 |
< |
methods with {\sc spme} is how much the agreement improves with increasing |
| 802 |
< |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 803 |
< |
appears to be reproducing the {\sc spme} torques most accurately. |
| 799 |
> |
Torques appear to be much more sensitive to charge interactions at |
| 800 |
> |
longer distances. The most noticeable feature in comparing the new |
| 801 |
> |
electrostatic methods with {\sc spme} is how much the agreement |
| 802 |
> |
improves with increasing cutoff radius. Again, the weakly damped and |
| 803 |
> |
undamped {\sc sf} method appears to reproduce the {\sc spme} torques |
| 804 |
> |
most accurately. |
| 805 |
|
|
| 806 |
|
Water molecules are dipolar, and the reaction field method reproduces |
| 807 |
|
the effect of the surrounding polarized medium on each of the |
| 810 |
|
|
| 811 |
|
\section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults} |
| 812 |
|
|
| 813 |
< |
It is clearly important that a new electrostatic method can reproduce |
| 814 |
< |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 815 |
< |
sum. However, the {\it directionality} of these vectors will also be |
| 816 |
< |
vital in calculating dynamical quantities accurately. Force and |
| 817 |
< |
torque directionalities were investigated by measuring the angles |
| 818 |
< |
formed between these vectors and the same vectors calculated using |
| 819 |
< |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 820 |
< |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 821 |
< |
distributions of the combined set over all system types. |
| 813 |
> |
It is clearly important that a new electrostatic method should be able |
| 814 |
> |
to reproduce the magnitudes of the force and torque vectors obtained |
| 815 |
> |
via the Ewald sum. However, the {\it directionality} of these vectors |
| 816 |
> |
will also be vital in calculating dynamical quantities accurately. |
| 817 |
> |
Force and torque directionalities were investigated by measuring the |
| 818 |
> |
angles formed between these vectors and the same vectors calculated |
| 819 |
> |
using {\sc spme}. The results (figure \ref{fig:frcTrqAng}) are compared |
| 820 |
> |
through the variance ($\sigma^2$) of the Gaussian fits of the angle |
| 821 |
> |
error distributions of the combined set over all system types. |
| 822 |
|
|
| 823 |
|
\begin{figure} |
| 824 |
|
\centering |
| 833 |
|
and 15~\AA\ = inverted triangles).} |
| 834 |
|
\label{fig:frcTrqAng} |
| 835 |
|
\end{figure} |
| 842 |
– |
|
| 836 |
|
Both the force and torque $\sigma^2$ results from the analysis of the |
| 837 |
|
total accumulated system data are tabulated in figure |
| 838 |
< |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 839 |
< |
sp}) method would be essentially unusable for molecular dynamics |
| 840 |
< |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 841 |
< |
method, however, is generating force and torque vectors which are |
| 842 |
< |
within a few degrees of the Ewald results even with weak (or no) |
| 850 |
< |
damping. |
| 838 |
> |
\ref{fig:frcTrqAng}. Here it is clear that the {\sc sp} method would |
| 839 |
> |
be essentially unusable for molecular dynamics unless the damping |
| 840 |
> |
function is added. The {\sc sf} method, however, is generating force |
| 841 |
> |
and torque vectors which are within a few degrees of the Ewald results |
| 842 |
> |
even with weak (or no) damping. |
| 843 |
|
|
| 844 |
|
All of the sets (aside from the over-damped case) show the improvement |
| 845 |
|
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 856 |
|
particles in all seven systems, while torque vectors are only |
| 857 |
|
available for neutral molecular groups. Damping is more beneficial to |
| 858 |
|
charged bodies, and this observation is investigated further in |
| 859 |
< |
section \ref{sec:IndividualResults}. |
| 859 |
> |
appendix \ref{app:IndividualResults}. |
| 860 |
|
|
| 861 |
|
Although not discussed previously, group based cutoffs can be applied |
| 862 |
|
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 921 |
|
The complimentary error function inserted into the potential weakens |
| 922 |
|
the electrostatic interaction as the value of $\alpha$ is increased. |
| 923 |
|
However, at larger values of $\alpha$, it is possible to over-damp the |
| 924 |
< |
electrostatic interaction and to remove it completely. Kast |
| 924 |
> |
electrostatic interaction and remove it completely. Kast |
| 925 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 926 |
|
values for these types of electrostatic summation methods by fitting |
| 927 |
|
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 928 |
|
0.25, and 0.16~\AA$^{-1}$ for cutoff values of 9, 12, and 15~\AA\ |
| 929 |
|
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 930 |
< |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 930 |
> |
obtain proper MC behavior (figure \ref{fig:delE}); however, based on |
| 931 |
|
these findings, choices this high would introduce error in the |
| 932 |
< |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 933 |
< |
observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial, |
| 934 |
< |
but damping may be unnecessary when using the {\sc sf} method. |
| 935 |
< |
|
| 944 |
< |
\section{Individual System Analysis Results}\label{sec:IndividualResults} |
| 945 |
< |
|
| 946 |
< |
The combined results of the previous sections show how the pairwise |
| 947 |
< |
methods compare to the Ewald summation in the general sense over all |
| 948 |
< |
of the system types. It is also useful to consider each of the |
| 949 |
< |
studied systems in an individual fashion, so that we can identify |
| 950 |
< |
conditions that are particularly difficult for a selected pairwise |
| 951 |
< |
method to address. This allows us to further establish the limitations |
| 952 |
< |
of these pairwise techniques. Below, the energy difference, force |
| 953 |
< |
vector, and torque vector analyses are presented on an individual |
| 954 |
< |
system basis. |
| 955 |
< |
|
| 956 |
< |
\subsection{SPC/E Water Results}\label{sec:WaterResults} |
| 957 |
< |
|
| 958 |
< |
The first system considered was liquid water at 300~K using the SPC/E |
| 959 |
< |
model of water.\cite{Berendsen87} The results for the energy gap |
| 960 |
< |
comparisons and the force and torque vector magnitude comparisons are |
| 961 |
< |
shown in table \ref{tab:spce}. The force and torque vector |
| 962 |
< |
directionality results are displayed separately in table |
| 963 |
< |
\ref{tab:spceAng}, where the effect of group-based cutoffs and |
| 964 |
< |
switching functions on the {\sc sp} and {\sc sf} potentials are also |
| 965 |
< |
investigated. In all of the individual results table, the method |
| 966 |
< |
abbreviations are as follows: |
| 967 |
< |
|
| 968 |
< |
\begin{itemize}[itemsep=0pt] |
| 969 |
< |
\item PC = Pure Cutoff, |
| 970 |
< |
\item SP = Shifted Potential, |
| 971 |
< |
\item SF = Shifted Force, |
| 972 |
< |
\item GSC = Group Switched Cutoff, |
| 973 |
< |
\item RF = Reaction Field (where $\varepsilon \approx\infty$), |
| 974 |
< |
\item GSSP = Group Switched Shifted Potential, and |
| 975 |
< |
\item GSSF = Group Switched Shifted Force. |
| 976 |
< |
\end{itemize} |
| 977 |
< |
|
| 978 |
< |
\begin{table}[htbp] |
| 979 |
< |
\centering |
| 980 |
< |
\caption{REGRESSION RESULTS OF THE LIQUID WATER SYSTEM FOR THE |
| 981 |
< |
$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it middle}) |
| 982 |
< |
AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 983 |
< |
|
| 984 |
< |
\footnotesize |
| 985 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 986 |
< |
\toprule |
| 987 |
< |
\toprule |
| 988 |
< |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 989 |
< |
\cmidrule(lr){3-4} |
| 990 |
< |
\cmidrule(lr){5-6} |
| 991 |
< |
\cmidrule(l){7-8} |
| 992 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 993 |
< |
\midrule |
| 994 |
< |
PC & & 3.046 & 0.002 & -3.018 & 0.002 & 4.719 & 0.005 \\ |
| 995 |
< |
SP & 0.0 & 1.035 & 0.218 & 0.908 & 0.313 & 1.037 & 0.470 \\ |
| 996 |
< |
& 0.1 & 1.021 & 0.387 & 0.965 & 0.752 & 1.006 & 0.947 \\ |
| 997 |
< |
& 0.2 & 0.997 & 0.962 & 1.001 & 0.994 & 0.994 & 0.996 \\ |
| 998 |
< |
& 0.3 & 0.984 & 0.980 & 0.997 & 0.985 & 0.982 & 0.987 \\ |
| 999 |
< |
SF & 0.0 & 0.977 & 0.974 & 0.996 & 0.992 & 0.991 & 0.997 \\ |
| 1000 |
< |
& 0.1 & 0.983 & 0.974 & 1.001 & 0.994 & 0.996 & 0.998 \\ |
| 1001 |
< |
& 0.2 & 0.992 & 0.989 & 1.001 & 0.995 & 0.994 & 0.996 \\ |
| 1002 |
< |
& 0.3 & 0.984 & 0.980 & 0.996 & 0.985 & 0.982 & 0.987 \\ |
| 1003 |
< |
GSC & & 0.918 & 0.862 & 0.852 & 0.756 & 0.801 & 0.700 \\ |
| 1004 |
< |
RF & & 0.971 & 0.958 & 0.975 & 0.987 & 0.959 & 0.983 \\ |
| 1005 |
< |
\midrule |
| 1006 |
< |
PC & & -1.647 & 0.000 & -0.127 & 0.000 & -0.979 & 0.000 \\ |
| 1007 |
< |
SP & 0.0 & 0.735 & 0.368 & 0.813 & 0.537 & 0.865 & 0.659 \\ |
| 1008 |
< |
& 0.1 & 0.850 & 0.612 & 0.956 & 0.887 & 0.992 & 0.979 \\ |
| 1009 |
< |
& 0.2 & 0.996 & 0.989 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1010 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ |
| 1011 |
< |
SF & 0.0 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 0.999 \\ |
| 1012 |
< |
& 0.1 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1013 |
< |
& 0.2 & 0.999 & 0.998 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1014 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ |
| 1015 |
< |
GSC & & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1016 |
< |
RF & & 0.999 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1017 |
< |
\midrule |
| 1018 |
< |
PC & & 2.387 & 0.000 & 0.183 & 0.000 & 1.282 & 0.000 \\ |
| 1019 |
< |
SP & 0.0 & 0.847 & 0.543 & 0.904 & 0.694 & 0.935 & 0.786 \\ |
| 1020 |
< |
& 0.1 & 0.922 & 0.749 & 0.980 & 0.934 & 0.996 & 0.988 \\ |
| 1021 |
< |
& 0.2 & 0.987 & 0.985 & 0.989 & 0.992 & 0.990 & 0.993 \\ |
| 1022 |
< |
& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ |
| 1023 |
< |
SF & 0.0 & 0.978 & 0.990 & 0.988 & 0.997 & 0.993 & 0.999 \\ |
| 1024 |
< |
& 0.1 & 0.983 & 0.991 & 0.993 & 0.997 & 0.997 & 0.999 \\ |
| 1025 |
< |
& 0.2 & 0.986 & 0.989 & 0.989 & 0.992 & 0.990 & 0.993 \\ |
| 1026 |
< |
& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ |
| 1027 |
< |
GSC & & 0.995 & 0.981 & 0.999 & 0.991 & 1.001 & 0.994 \\ |
| 1028 |
< |
RF & & 0.993 & 0.989 & 0.998 & 0.996 & 1.000 & 0.999 \\ |
| 1029 |
< |
\bottomrule |
| 1030 |
< |
\end{tabular} |
| 1031 |
< |
\label{tab:spce} |
| 1032 |
< |
\end{table} |
| 1033 |
< |
|
| 1034 |
< |
\begin{table}[htbp] |
| 1035 |
< |
\centering |
| 1036 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1037 |
< |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE LIQUID WATER |
| 1038 |
< |
SYSTEM} |
| 1039 |
< |
|
| 1040 |
< |
\footnotesize |
| 1041 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1042 |
< |
\toprule |
| 1043 |
< |
\toprule |
| 1044 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1045 |
< |
\cmidrule(lr){3-5} |
| 1046 |
< |
\cmidrule(l){6-8} |
| 1047 |
< |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1048 |
< |
\midrule |
| 1049 |
< |
PC & & 783.759 & 481.353 & 332.677 & 248.674 & 144.382 & 98.535 \\ |
| 1050 |
< |
SP & 0.0 & 659.440 & 380.699 & 250.002 & 235.151 & 134.661 & 88.135 \\ |
| 1051 |
< |
& 0.1 & 293.849 & 67.772 & 11.609 & 105.090 & 23.813 & 4.369 \\ |
| 1052 |
< |
& 0.2 & 5.975 & 0.136 & 0.094 & 5.553 & 1.784 & 1.536 \\ |
| 1053 |
< |
& 0.3 & 0.725 & 0.707 & 0.693 & 7.293 & 6.933 & 6.748 \\ |
| 1054 |
< |
SF & 0.0 & 2.238 & 0.713 & 0.292 & 3.290 & 1.090 & 0.416 \\ |
| 1055 |
< |
& 0.1 & 2.238 & 0.524 & 0.115 & 3.184 & 0.945 & 0.326 \\ |
| 1056 |
< |
& 0.2 & 0.374 & 0.102 & 0.094 & 2.598 & 1.755 & 1.537 \\ |
| 1057 |
< |
& 0.3 & 0.721 & 0.707 & 0.693 & 7.322 & 6.933 & 6.748 \\ |
| 1058 |
< |
GSC & & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ |
| 1059 |
< |
RF & & 2.091 & 0.403 & 0.113 & 3.583 & 1.071 & 0.399 \\ |
| 1060 |
< |
\midrule |
| 1061 |
< |
GSSP & 0.0 & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ |
| 1062 |
< |
& 0.1 & 1.879 & 0.291 & 0.057 & 3.983 & 1.117 & 0.370 \\ |
| 1063 |
< |
& 0.2 & 0.443 & 0.103 & 0.093 & 2.821 & 1.794 & 1.532 \\ |
| 1064 |
< |
& 0.3 & 0.728 & 0.694 & 0.692 & 7.387 & 6.942 & 6.748 \\ |
| 1065 |
< |
GSSF & 0.0 & 1.298 & 0.270 & 0.083 & 3.098 & 0.992 & 0.375 \\ |
| 1066 |
< |
& 0.1 & 1.296 & 0.210 & 0.044 & 3.055 & 0.922 & 0.330 \\ |
| 1067 |
< |
& 0.2 & 0.433 & 0.104 & 0.093 & 2.895 & 1.797 & 1.532 \\ |
| 1068 |
< |
& 0.3 & 0.728 & 0.694 & 0.692 & 7.410 & 6.942 & 6.748 \\ |
| 1069 |
< |
\bottomrule |
| 1070 |
< |
\end{tabular} |
| 1071 |
< |
\label{tab:spceAng} |
| 1072 |
< |
\end{table} |
| 1073 |
< |
|
| 1074 |
< |
The water results parallel the combined results seen in sections |
| 1075 |
< |
\ref{sec:EnergyResults} through \ref{sec:FTDirResults}. There is good |
| 1076 |
< |
agreement with {\sc spme} in both energetic and dynamic behavior when |
| 1077 |
< |
using the {\sc sf} method with and without damping. The {\sc sp} |
| 1078 |
< |
method does well with an $\alpha$ around 0.2~\AA$^{-1}$, particularly |
| 1079 |
< |
with cutoff radii greater than 12~\AA. Over-damping the electrostatics |
| 1080 |
< |
reduces the agreement between both these methods and {\sc spme}. |
| 1081 |
< |
|
| 1082 |
< |
The pure cutoff ({\sc pc}) method performs poorly, again mirroring the |
| 1083 |
< |
observations from the combined results. In contrast to these results, however, the use of a switching function and group |
| 1084 |
< |
based cutoffs greatly improves the results for these neutral water |
| 1085 |
< |
molecules. The group switched cutoff ({\sc gsc}) does not mimic the |
| 1086 |
< |
energetics of {\sc spme} as well as the {\sc sp} (with moderate |
| 1087 |
< |
damping) and {\sc sf} methods, but the dynamics are quite good. The |
| 1088 |
< |
switching functions correct discontinuities in the potential and |
| 1089 |
< |
forces, leading to these improved results. Such improvements with the |
| 1090 |
< |
use of a switching function have been recognized in previous |
| 1091 |
< |
studies,\cite{Andrea83,Steinbach94} and this proves to be a useful |
| 1092 |
< |
tactic for stably incorporating local area electrostatic effects. |
| 1093 |
< |
|
| 1094 |
< |
The reaction field ({\sc rf}) method simply extends upon the results |
| 1095 |
< |
observed in the {\sc gsc} case. Both methods are similar in form |
| 1096 |
< |
(i.e. neutral groups, switching function), but {\sc rf} incorporates |
| 1097 |
< |
an added effect from the external dielectric. This similarity |
| 1098 |
< |
translates into the same good dynamic results and improved energetic |
| 1099 |
< |
agreement with {\sc spme}. Though this agreement is not to the level |
| 1100 |
< |
of the moderately damped {\sc sp} and {\sc sf} methods, these results |
| 1101 |
< |
show how incorporating some implicit properties of the surroundings |
| 1102 |
< |
(i.e. $\epsilon_\textrm{S}$) can improve the solvent depiction. |
| 1103 |
< |
|
| 1104 |
< |
As a final note for the liquid water system, use of group cutoffs and a |
| 1105 |
< |
switching function leads to noticeable improvements in the {\sc sp} |
| 1106 |
< |
and {\sc sf} methods, primarily in directionality of the force and |
| 1107 |
< |
torque vectors (table \ref{tab:spceAng}). The {\sc sp} method shows |
| 1108 |
< |
significant narrowing of the angle distribution when using little to |
| 1109 |
< |
no damping and only modest improvement for the recommended conditions |
| 1110 |
< |
($\alpha = 0.2$~\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$~\AA). The |
| 1111 |
< |
{\sc sf} method shows modest narrowing across all damping and cutoff |
| 1112 |
< |
ranges of interest. When over-damping these methods, group cutoffs and |
| 1113 |
< |
the switching function do not improve the force and torque |
| 1114 |
< |
directionalities. |
| 1115 |
< |
|
| 1116 |
< |
\subsection{SPC/E Ice I$_\textrm{c}$ Results}\label{sec:IceResults} |
| 1117 |
< |
|
| 1118 |
< |
In addition to the disordered molecular system above, the ordered |
| 1119 |
< |
molecular system of ice I$_\textrm{c}$ was also considered. Ice |
| 1120 |
< |
polymorph could have been used to fit this role; however, ice |
| 1121 |
< |
I$_\textrm{c}$ was chosen because it can form an ideal periodic |
| 1122 |
< |
lattice with the same number of water molecules used in the disordered |
| 1123 |
< |
liquid state case. The results for the energy gap comparisons and the |
| 1124 |
< |
force and torque vector magnitude comparisons are shown in table |
| 1125 |
< |
\ref{tab:ice}. The force and torque vector directionality results are |
| 1126 |
< |
displayed separately in table \ref{tab:iceAng}, where the effect of |
| 1127 |
< |
group-based cutoffs and switching functions on the {\sc sp} and {\sc |
| 1128 |
< |
sf} potentials are also displayed. |
| 1129 |
< |
|
| 1130 |
< |
\begin{table}[htbp] |
| 1131 |
< |
\centering |
| 1132 |
< |
\caption{REGRESSION RESULTS OF THE ICE I$_\textrm{c}$ SYSTEM FOR |
| 1133 |
< |
$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it |
| 1134 |
< |
middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1135 |
< |
|
| 1136 |
< |
\footnotesize |
| 1137 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1138 |
< |
\toprule |
| 1139 |
< |
\toprule |
| 1140 |
< |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1141 |
< |
\cmidrule(lr){3-4} |
| 1142 |
< |
\cmidrule(lr){5-6} |
| 1143 |
< |
\cmidrule(l){7-8} |
| 1144 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1145 |
< |
\midrule |
| 1146 |
< |
PC & & 19.897 & 0.047 & -29.214 & 0.048 & -3.771 & 0.001 \\ |
| 1147 |
< |
SP & 0.0 & -0.014 & 0.000 & 2.135 & 0.347 & 0.457 & 0.045 \\ |
| 1148 |
< |
& 0.1 & 0.321 & 0.017 & 1.490 & 0.584 & 0.886 & 0.796 \\ |
| 1149 |
< |
& 0.2 & 0.896 & 0.872 & 1.011 & 0.998 & 0.997 & 0.999 \\ |
| 1150 |
< |
& 0.3 & 0.983 & 0.997 & 0.992 & 0.997 & 0.991 & 0.997 \\ |
| 1151 |
< |
SF & 0.0 & 0.943 & 0.979 & 1.048 & 0.978 & 0.995 & 0.999 \\ |
| 1152 |
< |
& 0.1 & 0.948 & 0.979 & 1.044 & 0.983 & 1.000 & 0.999 \\ |
| 1153 |
< |
& 0.2 & 0.982 & 0.997 & 0.969 & 0.960 & 0.997 & 0.999 \\ |
| 1154 |
< |
& 0.3 & 0.985 & 0.997 & 0.961 & 0.961 & 0.991 & 0.997 \\ |
| 1155 |
< |
GSC & & 0.983 & 0.985 & 0.966 & 0.994 & 1.003 & 0.999 \\ |
| 1156 |
< |
RF & & 0.924 & 0.944 & 0.990 & 0.996 & 0.991 & 0.998 \\ |
| 1157 |
< |
\midrule |
| 1158 |
< |
PC & & -4.375 & 0.000 & 6.781 & 0.000 & -3.369 & 0.000 \\ |
| 1159 |
< |
SP & 0.0 & 0.515 & 0.164 & 0.856 & 0.426 & 0.743 & 0.478 \\ |
| 1160 |
< |
& 0.1 & 0.696 & 0.405 & 0.977 & 0.817 & 0.974 & 0.964 \\ |
| 1161 |
< |
& 0.2 & 0.981 & 0.980 & 1.001 & 1.000 & 1.000 & 1.000 \\ |
| 1162 |
< |
& 0.3 & 0.996 & 0.998 & 0.997 & 0.999 & 0.997 & 0.999 \\ |
| 1163 |
< |
SF & 0.0 & 0.991 & 0.995 & 1.003 & 0.998 & 0.999 & 1.000 \\ |
| 1164 |
< |
& 0.1 & 0.992 & 0.995 & 1.003 & 0.998 & 1.000 & 1.000 \\ |
| 1165 |
< |
& 0.2 & 0.998 & 0.998 & 0.981 & 0.962 & 1.000 & 1.000 \\ |
| 1166 |
< |
& 0.3 & 0.996 & 0.998 & 0.976 & 0.957 & 0.997 & 0.999 \\ |
| 1167 |
< |
GSC & & 0.997 & 0.996 & 0.998 & 0.999 & 1.000 & 1.000 \\ |
| 1168 |
< |
RF & & 0.988 & 0.989 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1169 |
< |
\midrule |
| 1170 |
< |
PC & & -6.367 & 0.000 & -3.552 & 0.000 & -3.447 & 0.000 \\ |
| 1171 |
< |
SP & 0.0 & 0.643 & 0.409 & 0.833 & 0.607 & 0.961 & 0.805 \\ |
| 1172 |
< |
& 0.1 & 0.791 & 0.683 & 0.957 & 0.914 & 1.000 & 0.989 \\ |
| 1173 |
< |
& 0.2 & 0.974 & 0.991 & 0.993 & 0.998 & 0.993 & 0.998 \\ |
| 1174 |
< |
& 0.3 & 0.976 & 0.992 & 0.977 & 0.992 & 0.977 & 0.992 \\ |
| 1175 |
< |
SF & 0.0 & 0.979 & 0.997 & 0.992 & 0.999 & 0.994 & 1.000 \\ |
| 1176 |
< |
& 0.1 & 0.984 & 0.997 & 0.996 & 0.999 & 0.998 & 1.000 \\ |
| 1177 |
< |
& 0.2 & 0.991 & 0.997 & 0.974 & 0.958 & 0.993 & 0.998 \\ |
| 1178 |
< |
& 0.3 & 0.977 & 0.992 & 0.956 & 0.948 & 0.977 & 0.992 \\ |
| 1179 |
< |
GSC & & 0.999 & 0.997 & 0.996 & 0.999 & 1.002 & 1.000 \\ |
| 1180 |
< |
RF & & 0.994 & 0.997 & 0.997 & 0.999 & 1.000 & 1.000 \\ |
| 1181 |
< |
\bottomrule |
| 1182 |
< |
\end{tabular} |
| 1183 |
< |
\label{tab:ice} |
| 1184 |
< |
\end{table} |
| 1185 |
< |
|
| 1186 |
< |
\begin{table}[htbp] |
| 1187 |
< |
\centering |
| 1188 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1189 |
< |
OF THE FORCE AND TORQUE VECTORS IN THE ICE I$_\textrm{c}$ SYSTEM} |
| 1190 |
< |
|
| 1191 |
< |
\footnotesize |
| 1192 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1193 |
< |
\toprule |
| 1194 |
< |
\toprule |
| 1195 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque |
| 1196 |
< |
$\sigma^2$} \\ |
| 1197 |
< |
\cmidrule(lr){3-5} |
| 1198 |
< |
\cmidrule(l){6-8} |
| 1199 |
< |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1200 |
< |
\midrule |
| 1201 |
< |
PC & & 2128.921 & 603.197 & 715.579 & 329.056 & 221.397 & 81.042 \\ |
| 1202 |
< |
SP & 0.0 & 1429.341 & 470.320 & 447.557 & 301.678 & 197.437 & 73.840 \\ |
| 1203 |
< |
& 0.1 & 590.008 & 107.510 & 18.883 & 118.201 & 32.472 & 3.599 \\ |
| 1204 |
< |
& 0.2 & 10.057 & 0.105 & 0.038 & 2.875 & 0.572 & 0.518 \\ |
| 1205 |
< |
& 0.3 & 0.245 & 0.260 & 0.262 & 2.365 & 2.396 & 2.327 \\ |
| 1206 |
< |
SF & 0.0 & 1.745 & 1.161 & 0.212 & 1.135 & 0.426 & 0.155 \\ |
| 1207 |
< |
& 0.1 & 1.721 & 0.868 & 0.082 & 1.118 & 0.358 & 0.118 \\ |
| 1208 |
< |
& 0.2 & 0.201 & 0.040 & 0.038 & 0.786 & 0.555 & 0.518 \\ |
| 1209 |
< |
& 0.3 & 0.241 & 0.260 & 0.262 & 2.368 & 2.400 & 2.327 \\ |
| 1210 |
< |
GSC & & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ |
| 1211 |
< |
RF & & 2.887 & 0.217 & 0.107 & 1.006 & 0.281 & 0.085 \\ |
| 1212 |
< |
\midrule |
| 1213 |
< |
GSSP & 0.0 & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ |
| 1214 |
< |
& 0.1 & 1.341 & 0.123 & 0.037 & 0.835 & 0.234 & 0.085 \\ |
| 1215 |
< |
& 0.2 & 0.558 & 0.040 & 0.037 & 0.823 & 0.557 & 0.519 \\ |
| 1216 |
< |
& 0.3 & 0.250 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ |
| 1217 |
< |
GSSF & 0.0 & 2.124 & 0.132 & 0.069 & 0.919 & 0.263 & 0.099 \\ |
| 1218 |
< |
& 0.1 & 2.165 & 0.101 & 0.035 & 0.895 & 0.244 & 0.096 \\ |
| 1219 |
< |
& 0.2 & 0.706 & 0.040 & 0.037 & 0.870 & 0.559 & 0.519 \\ |
| 1220 |
< |
& 0.3 & 0.251 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ |
| 1221 |
< |
\bottomrule |
| 1222 |
< |
\end{tabular} |
| 1223 |
< |
\label{tab:iceAng} |
| 1224 |
< |
\end{table} |
| 1225 |
< |
|
| 1226 |
< |
Highly ordered systems are a difficult test for the pairwise methods |
| 1227 |
< |
in that they lack the implicit periodicity of the Ewald summation. As |
| 1228 |
< |
expected, the energy gap agreement with {\sc spme} is reduced for the |
| 1229 |
< |
{\sc sp} and {\sc sf} methods with parameters that were ideal for the |
| 1230 |
< |
disordered liquid system. Moving to higher $R_\textrm{c}$ helps |
| 1231 |
< |
improve the agreement, though at an increase in computational cost. |
| 1232 |
< |
The dynamics of this crystalline system (both in magnitude and |
| 1233 |
< |
direction) are little affected. Both methods still reproduce the Ewald |
| 1234 |
< |
behavior with the same parameter recommendations from the previous |
| 1235 |
< |
section. |
| 1236 |
< |
|
| 1237 |
< |
It is also worth noting that {\sc rf} exhibits improved energy gap |
| 1238 |
< |
results over the liquid water system. One possible explanation is |
| 1239 |
< |
that the ice I$_\textrm{c}$ crystal is ordered such that the net |
| 1240 |
< |
dipole moment of the crystal is zero. With $\epsilon_\textrm{S} = |
| 1241 |
< |
\infty$, the reaction field incorporates this structural organization |
| 1242 |
< |
by actively enforcing a zeroed dipole moment within each cutoff |
| 1243 |
< |
sphere. |
| 1244 |
< |
|
| 1245 |
< |
\subsection{NaCl Melt Results}\label{sec:SaltMeltResults} |
| 1246 |
< |
|
| 1247 |
< |
A high temperature NaCl melt was tested to gauge the accuracy of the |
| 1248 |
< |
pairwise summation methods in a disordered system of charges. The |
| 1249 |
< |
results for the energy gap comparisons and the force vector magnitude |
| 1250 |
< |
comparisons are shown in table \ref{tab:melt}. The force vector |
| 1251 |
< |
directionality results are displayed separately in table |
| 1252 |
< |
\ref{tab:meltAng}. |
| 1253 |
< |
|
| 1254 |
< |
\begin{table}[htbp] |
| 1255 |
< |
\centering |
| 1256 |
< |
\caption{REGRESSION RESULTS OF THE MOLTEN SODIUM CHLORIDE SYSTEM FOR |
| 1257 |
< |
$\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES ({\it |
| 1258 |
< |
lower})} |
| 1259 |
< |
|
| 1260 |
< |
\footnotesize |
| 1261 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1262 |
< |
\toprule |
| 1263 |
< |
\toprule |
| 1264 |
< |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1265 |
< |
\cmidrule(lr){3-4} |
| 1266 |
< |
\cmidrule(lr){5-6} |
| 1267 |
< |
\cmidrule(l){7-8} |
| 1268 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1269 |
< |
\midrule |
| 1270 |
< |
PC & & -0.008 & 0.000 & -0.049 & 0.005 & -0.136 & 0.020 \\ |
| 1271 |
< |
SP & 0.0 & 0.928 & 0.996 & 0.931 & 0.998 & 0.950 & 0.999 \\ |
| 1272 |
< |
& 0.1 & 0.977 & 0.998 & 0.998 & 1.000 & 0.997 & 1.000 \\ |
| 1273 |
< |
& 0.2 & 0.960 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ |
| 1274 |
< |
& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ |
| 1275 |
< |
SF & 0.0 & 0.996 & 1.000 & 0.995 & 1.000 & 0.997 & 1.000 \\ |
| 1276 |
< |
& 0.1 & 1.021 & 1.000 & 1.024 & 1.000 & 1.007 & 1.000 \\ |
| 1277 |
< |
& 0.2 & 0.966 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ |
| 1278 |
< |
& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ |
| 1279 |
< |
\midrule |
| 1280 |
< |
PC & & 1.103 & 0.000 & 0.989 & 0.000 & 0.802 & 0.000 \\ |
| 1281 |
< |
SP & 0.0 & 0.973 & 0.981 & 0.975 & 0.988 & 0.979 & 0.992 \\ |
| 1282 |
< |
& 0.1 & 0.987 & 0.992 & 0.993 & 0.998 & 0.997 & 0.999 \\ |
| 1283 |
< |
& 0.2 & 0.993 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ |
| 1284 |
< |
& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ |
| 1285 |
< |
SF & 0.0 & 0.996 & 0.997 & 0.997 & 0.999 & 0.998 & 1.000 \\ |
| 1286 |
< |
& 0.1 & 1.000 & 0.997 & 1.001 & 0.999 & 1.000 & 1.000 \\ |
| 1287 |
< |
& 0.2 & 0.994 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ |
| 1288 |
< |
& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ |
| 1289 |
< |
\bottomrule |
| 1290 |
< |
\end{tabular} |
| 1291 |
< |
\label{tab:melt} |
| 1292 |
< |
\end{table} |
| 1293 |
< |
|
| 1294 |
< |
\begin{table}[htbp] |
| 1295 |
< |
\centering |
| 1296 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1297 |
< |
OF THE FORCE VECTORS IN THE MOLTEN SODIUM CHLORIDE SYSTEM} |
| 1298 |
< |
|
| 1299 |
< |
\footnotesize |
| 1300 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1301 |
< |
\toprule |
| 1302 |
< |
\toprule |
| 1303 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1304 |
< |
\cmidrule(lr){3-5} |
| 1305 |
< |
\cmidrule(l){6-8} |
| 1306 |
< |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA \\ |
| 1307 |
< |
\midrule |
| 1308 |
< |
PC & & 13.294 & 8.035 & 5.366 \\ |
| 1309 |
< |
SP & 0.0 & 13.316 & 8.037 & 5.385 \\ |
| 1310 |
< |
& 0.1 & 5.705 & 1.391 & 0.360 \\ |
| 1311 |
< |
& 0.2 & 2.415 & 7.534 & 13.927 \\ |
| 1312 |
< |
& 0.3 & 23.769 & 67.306 & 57.252 \\ |
| 1313 |
< |
SF & 0.0 & 1.693 & 0.603 & 0.256 \\ |
| 1314 |
< |
& 0.1 & 1.687 & 0.653 & 0.272 \\ |
| 1315 |
< |
& 0.2 & 2.598 & 7.523 & 13.930 \\ |
| 1316 |
< |
& 0.3 & 23.734 & 67.305 & 57.252 \\ |
| 1317 |
< |
\bottomrule |
| 1318 |
< |
\end{tabular} |
| 1319 |
< |
\label{tab:meltAng} |
| 1320 |
< |
\end{table} |
| 1321 |
< |
|
| 1322 |
< |
The molten NaCl system shows more sensitivity to the electrostatic |
| 1323 |
< |
damping than the water systems. The most noticeable point is that the |
| 1324 |
< |
undamped {\sc sf} method does very well at replicating the {\sc spme} |
| 1325 |
< |
configurational energy differences and forces. Light damping appears |
| 1326 |
< |
to minimally improve the dynamics, but this comes with a deterioration |
| 1327 |
< |
of the energy gap results. In contrast, this light damping improves |
| 1328 |
< |
the {\sc sp} energy gaps and forces. Moderate and heavy electrostatic |
| 1329 |
< |
damping reduce the agreement with {\sc spme} for both methods. From |
| 1330 |
< |
these observations, the undamped {\sc sf} method is the best choice |
| 1331 |
< |
for disordered systems of charges. |
| 1332 |
< |
|
| 1333 |
< |
\subsection{NaCl Crystal Results}\label{sec:SaltCrystalResults} |
| 1334 |
< |
|
| 1335 |
< |
Similar to the use of ice I$_\textrm{c}$ to investigate the role of |
| 1336 |
< |
order in molecular systems on the effectiveness of the pairwise |
| 1337 |
< |
methods, the 1000~K NaCl crystal system was used to investigate the |
| 1338 |
< |
accuracy of the pairwise summation methods in an ordered system of |
| 1339 |
< |
charged particles. The results for the energy gap comparisons and the |
| 1340 |
< |
force vector magnitude comparisons are shown in table \ref{tab:salt}. |
| 1341 |
< |
The force vector directionality results are displayed separately in |
| 1342 |
< |
table \ref{tab:saltAng}. |
| 1343 |
< |
|
| 1344 |
< |
\begin{table}[htbp] |
| 1345 |
< |
\centering |
| 1346 |
< |
\caption{REGRESSION RESULTS OF THE CRYSTALLINE SODIUM CHLORIDE |
| 1347 |
< |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES |
| 1348 |
< |
({\it lower})} |
| 1349 |
< |
|
| 1350 |
< |
\footnotesize |
| 1351 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1352 |
< |
\toprule |
| 1353 |
< |
\toprule |
| 1354 |
< |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1355 |
< |
\cmidrule(lr){3-4} |
| 1356 |
< |
\cmidrule(lr){5-6} |
| 1357 |
< |
\cmidrule(l){7-8} |
| 1358 |
< |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1359 |
< |
\midrule |
| 1360 |
< |
PC & & -20.241 & 0.228 & -20.248 & 0.229 & -20.239 & 0.228 \\ |
| 1361 |
< |
SP & 0.0 & 1.039 & 0.733 & 2.037 & 0.565 & 1.225 & 0.743 \\ |
| 1362 |
< |
& 0.1 & 1.049 & 0.865 & 1.424 & 0.784 & 1.029 & 0.980 \\ |
| 1363 |
< |
& 0.2 & 0.982 & 0.976 & 0.969 & 0.980 & 0.960 & 0.980 \\ |
| 1364 |
< |
& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.945 \\ |
| 1365 |
< |
SF & 0.0 & 1.041 & 0.967 & 0.994 & 0.989 & 0.957 & 0.993 \\ |
| 1366 |
< |
& 0.1 & 1.050 & 0.968 & 0.996 & 0.991 & 0.972 & 0.995 \\ |
| 1367 |
< |
& 0.2 & 0.982 & 0.975 & 0.959 & 0.980 & 0.960 & 0.980 \\ |
| 1368 |
< |
& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.944 \\ |
| 1369 |
< |
\midrule |
| 1370 |
< |
PC & & 0.795 & 0.000 & 0.792 & 0.000 & 0.793 & 0.000 \\ |
| 1371 |
< |
SP & 0.0 & 0.916 & 0.829 & 1.086 & 0.791 & 1.010 & 0.936 \\ |
| 1372 |
< |
& 0.1 & 0.958 & 0.917 & 1.049 & 0.943 & 1.001 & 0.995 \\ |
| 1373 |
< |
& 0.2 & 0.981 & 0.981 & 0.982 & 0.984 & 0.981 & 0.984 \\ |
| 1374 |
< |
& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ |
| 1375 |
< |
SF & 0.0 & 1.002 & 0.983 & 0.997 & 0.994 & 0.991 & 0.997 \\ |
| 1376 |
< |
& 0.1 & 1.003 & 0.984 & 0.996 & 0.995 & 0.993 & 0.997 \\ |
| 1377 |
< |
& 0.2 & 0.983 & 0.980 & 0.981 & 0.984 & 0.981 & 0.984 \\ |
| 1378 |
< |
& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ |
| 1379 |
< |
\bottomrule |
| 1380 |
< |
\end{tabular} |
| 1381 |
< |
\label{tab:salt} |
| 1382 |
< |
\end{table} |
| 1383 |
< |
|
| 1384 |
< |
\begin{table}[htbp] |
| 1385 |
< |
\centering |
| 1386 |
< |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1387 |
< |
DISTRIBUTIONS OF THE FORCE VECTORS IN THE CRYSTALLINE SODIUM CHLORIDE |
| 1388 |
< |
SYSTEM} |
| 1389 |
< |
|
| 1390 |
< |
\footnotesize |
| 1391 |
< |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1392 |
< |
\toprule |
| 1393 |
< |
\toprule |
| 1394 |
< |
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1395 |
< |
\cmidrule(lr){3-5} |
| 1396 |
< |
\cmidrule(l){6-8} |
| 1397 |
< |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA \\ |
| 1398 |
< |
\midrule |
| 1399 |
< |
PC & & 111.945 & 111.824 & 111.866 \\ |
| 1400 |
< |
SP & 0.0 & 112.414 & 152.215 & 38.087 \\ |
| 1401 |
< |
& 0.1 & 52.361 & 42.574 & 2.819 \\ |
| 1402 |
< |
& 0.2 & 10.847 & 9.709 & 9.686 \\ |
| 1403 |
< |
& 0.3 & 31.128 & 31.104 & 31.029 \\ |
| 1404 |
< |
SF & 0.0 & 10.025 & 3.555 & 1.648 \\ |
| 1405 |
< |
& 0.1 & 9.462 & 3.303 & 1.721 \\ |
| 1406 |
< |
& 0.2 & 11.454 & 9.813 & 9.701 \\ |
| 1407 |
< |
& 0.3 & 31.120 & 31.105 & 31.029 \\ |
| 1408 |
< |
\bottomrule |
| 1409 |
< |
\end{tabular} |
| 1410 |
< |
\label{tab:saltAng} |
| 1411 |
< |
\end{table} |
| 1412 |
< |
|
| 1413 |
< |
The crystalline NaCl system is the most challenging test case for the |
| 1414 |
< |
pairwise summation methods, as evidenced by the results in tables |
| 1415 |
< |
\ref{tab:salt} and \ref{tab:saltAng}. The undamped and weakly damped |
| 1416 |
< |
{\sc sf} methods seem to be the best choices. These methods match well |
| 1417 |
< |
with {\sc spme} across the energy gap, force magnitude, and force |
| 1418 |
< |
directionality tests. The {\sc sp} method struggles in all cases, |
| 1419 |
< |
with the exception of good dynamics reproduction when using weak |
| 1420 |
< |
electrostatic damping with a large cutoff radius. |
| 1421 |
< |
|
| 1422 |
< |
The moderate electrostatic damping case is not as good as we would |
| 1423 |
< |
expect given the long-time dynamics results observed for this system |
| 1424 |
< |
(see section \ref{sec:LongTimeDynamics}). Since the data tabulated in |
| 1425 |
< |
tables \ref{tab:salt} and \ref{tab:saltAng} are a test of |
| 1426 |
< |
instantaneous dynamics, this indicates that good long-time dynamics |
| 1427 |
< |
comes in part at the expense of short-time dynamics. |
| 932 |
> |
molecular torques, particularly for the shorter cutoffs. Based on the |
| 933 |
> |
above observations, empirical damping up to 0.2~\AA$^{-1}$ is |
| 934 |
> |
beneficial, but damping may be unnecessary when using the {\sc sf} |
| 935 |
> |
method. |
| 936 |
|
|
| 1429 |
– |
\subsection{0.11M NaCl Solution Results} |
| 937 |
|
|
| 1431 |
– |
In an effort to bridge the charged atomic and neutral molecular |
| 1432 |
– |
systems, Na$^+$ and Cl$^-$ ion charge defects were incorporated into |
| 1433 |
– |
the liquid water system. This low ionic strength system consists of 4 |
| 1434 |
– |
ions in the 1000 SPC/E water solvent ($\approx$0.11 M). The results |
| 1435 |
– |
for the energy gap comparisons and the force and torque vector |
| 1436 |
– |
magnitude comparisons are shown in table \ref{tab:solnWeak}. The |
| 1437 |
– |
force and torque vector directionality results are displayed |
| 1438 |
– |
separately in table \ref{tab:solnWeakAng}, where the effect of |
| 1439 |
– |
group-based cutoffs and switching functions on the {\sc sp} and {\sc |
| 1440 |
– |
sf} potentials are investigated. |
| 1441 |
– |
|
| 1442 |
– |
\begin{table}[htbp] |
| 1443 |
– |
\centering |
| 1444 |
– |
\caption{REGRESSION RESULTS OF THE WEAK SODIUM CHLORIDE SOLUTION |
| 1445 |
– |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES |
| 1446 |
– |
({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1447 |
– |
|
| 1448 |
– |
\footnotesize |
| 1449 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1450 |
– |
\toprule |
| 1451 |
– |
\toprule |
| 1452 |
– |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1453 |
– |
\cmidrule(lr){3-4} |
| 1454 |
– |
\cmidrule(lr){5-6} |
| 1455 |
– |
\cmidrule(l){7-8} |
| 1456 |
– |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1457 |
– |
\midrule |
| 1458 |
– |
PC & & 0.247 & 0.000 & -1.103 & 0.001 & 5.480 & 0.015 \\ |
| 1459 |
– |
SP & 0.0 & 0.935 & 0.388 & 0.984 & 0.541 & 1.010 & 0.685 \\ |
| 1460 |
– |
& 0.1 & 0.951 & 0.603 & 0.993 & 0.875 & 1.001 & 0.979 \\ |
| 1461 |
– |
& 0.2 & 0.969 & 0.968 & 0.996 & 0.997 & 0.994 & 0.997 \\ |
| 1462 |
– |
& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ |
| 1463 |
– |
SF & 0.0 & 0.963 & 0.971 & 0.989 & 0.996 & 0.991 & 0.998 \\ |
| 1464 |
– |
& 0.1 & 0.970 & 0.971 & 0.995 & 0.997 & 0.997 & 0.999 \\ |
| 1465 |
– |
& 0.2 & 0.972 & 0.975 & 0.996 & 0.997 & 0.994 & 0.997 \\ |
| 1466 |
– |
& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ |
| 1467 |
– |
GSC & & 0.964 & 0.731 & 0.984 & 0.704 & 1.005 & 0.770 \\ |
| 1468 |
– |
RF & & 0.968 & 0.605 & 0.974 & 0.541 & 1.014 & 0.614 \\ |
| 1469 |
– |
\midrule |
| 1470 |
– |
PC & & 1.354 & 0.000 & -1.190 & 0.000 & -0.314 & 0.000 \\ |
| 1471 |
– |
SP & 0.0 & 0.720 & 0.338 & 0.808 & 0.523 & 0.860 & 0.643 \\ |
| 1472 |
– |
& 0.1 & 0.839 & 0.583 & 0.955 & 0.882 & 0.992 & 0.978 \\ |
| 1473 |
– |
& 0.2 & 0.995 & 0.987 & 0.999 & 1.000 & 0.999 & 1.000 \\ |
| 1474 |
– |
& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ |
| 1475 |
– |
SF & 0.0 & 0.998 & 0.994 & 1.000 & 0.998 & 1.000 & 0.999 \\ |
| 1476 |
– |
& 0.1 & 0.997 & 0.994 & 1.000 & 0.999 & 1.000 & 1.000 \\ |
| 1477 |
– |
& 0.2 & 0.999 & 0.998 & 0.999 & 1.000 & 0.999 & 1.000 \\ |
| 1478 |
– |
& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ |
| 1479 |
– |
GSC & & 0.995 & 0.990 & 0.998 & 0.997 & 0.998 & 0.996 \\ |
| 1480 |
– |
RF & & 0.998 & 0.993 & 0.999 & 0.998 & 0.999 & 0.996 \\ |
| 1481 |
– |
\midrule |
| 1482 |
– |
PC & & 2.437 & 0.000 & -1.872 & 0.000 & 2.138 & 0.000 \\ |
| 1483 |
– |
SP & 0.0 & 0.838 & 0.525 & 0.901 & 0.686 & 0.932 & 0.779 \\ |
| 1484 |
– |
& 0.1 & 0.914 & 0.733 & 0.979 & 0.932 & 0.995 & 0.987 \\ |
| 1485 |
– |
& 0.2 & 0.977 & 0.969 & 0.988 & 0.990 & 0.989 & 0.990 \\ |
| 1486 |
– |
& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ |
| 1487 |
– |
SF & 0.0 & 0.969 & 0.977 & 0.987 & 0.996 & 0.993 & 0.998 \\ |
| 1488 |
– |
& 0.1 & 0.975 & 0.978 & 0.993 & 0.996 & 0.997 & 0.998 \\ |
| 1489 |
– |
& 0.2 & 0.976 & 0.973 & 0.988 & 0.990 & 0.989 & 0.990 \\ |
| 1490 |
– |
& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ |
| 1491 |
– |
GSC & & 0.980 & 0.959 & 0.990 & 0.983 & 0.992 & 0.989 \\ |
| 1492 |
– |
RF & & 0.984 & 0.975 & 0.996 & 0.995 & 0.998 & 0.998 \\ |
| 1493 |
– |
\bottomrule |
| 1494 |
– |
\end{tabular} |
| 1495 |
– |
\label{tab:solnWeak} |
| 1496 |
– |
\end{table} |
| 1497 |
– |
|
| 1498 |
– |
\begin{table}[htbp] |
| 1499 |
– |
\centering |
| 1500 |
– |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1501 |
– |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE WEAK SODIUM |
| 1502 |
– |
CHLORIDE SOLUTION SYSTEM} |
| 1503 |
– |
|
| 1504 |
– |
\footnotesize |
| 1505 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1506 |
– |
\toprule |
| 1507 |
– |
\toprule |
| 1508 |
– |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1509 |
– |
\cmidrule(lr){3-5} |
| 1510 |
– |
\cmidrule(l){6-8} |
| 1511 |
– |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1512 |
– |
\midrule |
| 1513 |
– |
PC & & 882.863 & 510.435 & 344.201 & 277.691 & 154.231 & 100.131 \\ |
| 1514 |
– |
SP & 0.0 & 732.569 & 405.704 & 257.756 & 261.445 & 142.245 & 91.497 \\ |
| 1515 |
– |
& 0.1 & 329.031 & 70.746 & 12.014 & 118.496 & 25.218 & 4.711 \\ |
| 1516 |
– |
& 0.2 & 6.772 & 0.153 & 0.118 & 9.780 & 2.101 & 2.102 \\ |
| 1517 |
– |
& 0.3 & 0.951 & 0.774 & 0.784 & 12.108 & 7.673 & 7.851 \\ |
| 1518 |
– |
SF & 0.0 & 2.555 & 0.762 & 0.313 & 6.590 & 1.328 & 0.558 \\ |
| 1519 |
– |
& 0.1 & 2.561 & 0.560 & 0.123 & 6.464 & 1.162 & 0.457 \\ |
| 1520 |
– |
& 0.2 & 0.501 & 0.118 & 0.118 & 5.698 & 2.074 & 2.099 \\ |
| 1521 |
– |
& 0.3 & 0.943 & 0.774 & 0.784 & 12.118 & 7.674 & 7.851 \\ |
| 1522 |
– |
GSC & & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ |
| 1523 |
– |
RF & & 2.415 & 0.452 & 0.130 & 6.915 & 1.423 & 0.507 \\ |
| 1524 |
– |
\midrule |
| 1525 |
– |
GSSP & 0.0 & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ |
| 1526 |
– |
& 0.1 & 2.251 & 0.324 & 0.064 & 7.628 & 1.639 & 0.497 \\ |
| 1527 |
– |
& 0.2 & 0.590 & 0.118 & 0.116 & 6.080 & 2.096 & 2.103 \\ |
| 1528 |
– |
& 0.3 & 0.953 & 0.759 & 0.780 & 12.347 & 7.683 & 7.849 \\ |
| 1529 |
– |
GSSF & 0.0 & 1.541 & 0.301 & 0.096 & 6.407 & 1.316 & 0.496 \\ |
| 1530 |
– |
& 0.1 & 1.541 & 0.237 & 0.050 & 6.356 & 1.202 & 0.457 \\ |
| 1531 |
– |
& 0.2 & 0.568 & 0.118 & 0.116 & 6.166 & 2.105 & 2.105 \\ |
| 1532 |
– |
& 0.3 & 0.954 & 0.759 & 0.780 & 12.337 & 7.684 & 7.849 \\ |
| 1533 |
– |
\bottomrule |
| 1534 |
– |
\end{tabular} |
| 1535 |
– |
\label{tab:solnWeakAng} |
| 1536 |
– |
\end{table} |
| 1537 |
– |
|
| 1538 |
– |
Because this system is a perturbation of the pure liquid water system, |
| 1539 |
– |
comparisons are best drawn between these two sets. The {\sc sp} and |
| 1540 |
– |
{\sc sf} methods are not significantly affected by the inclusion of a |
| 1541 |
– |
few ions. The aspect of cutoff sphere neutralization aids in the |
| 1542 |
– |
smooth incorporation of these ions; thus, all of the observations |
| 1543 |
– |
regarding these methods carry over from section |
| 1544 |
– |
\ref{sec:WaterResults}. The differences between these systems are more |
| 1545 |
– |
visible for the {\sc rf} method. Though good force agreement is still |
| 1546 |
– |
maintained, the energy gaps show a significant increase in the scatter |
| 1547 |
– |
of the data. |
| 1548 |
– |
|
| 1549 |
– |
\subsection{1.1M NaCl Solution Results} |
| 1550 |
– |
|
| 1551 |
– |
The bridging of the charged atomic and neutral molecular systems was |
| 1552 |
– |
further developed by considering a high ionic strength system |
| 1553 |
– |
consisting of 40 ions in the 1000 SPC/E water solvent ($\approx$1.1 |
| 1554 |
– |
M). The results for the energy gap comparisons and the force and |
| 1555 |
– |
torque vector magnitude comparisons are shown in table |
| 1556 |
– |
\ref{tab:solnStr}. The force and torque vector directionality |
| 1557 |
– |
results are displayed separately in table \ref{tab:solnStrAng}, where |
| 1558 |
– |
the effect of group-based cutoffs and switching functions on the {\sc |
| 1559 |
– |
sp} and {\sc sf} potentials are investigated. |
| 1560 |
– |
|
| 1561 |
– |
\begin{table}[htbp] |
| 1562 |
– |
\centering |
| 1563 |
– |
\caption{REGRESSION RESULTS OF THE STRONG SODIUM CHLORIDE SOLUTION |
| 1564 |
– |
SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES |
| 1565 |
– |
({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1566 |
– |
|
| 1567 |
– |
\footnotesize |
| 1568 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1569 |
– |
\toprule |
| 1570 |
– |
\toprule |
| 1571 |
– |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1572 |
– |
\cmidrule(lr){3-4} |
| 1573 |
– |
\cmidrule(lr){5-6} |
| 1574 |
– |
\cmidrule(l){7-8} |
| 1575 |
– |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1576 |
– |
\midrule |
| 1577 |
– |
PC & & -0.081 & 0.000 & 0.945 & 0.001 & 0.073 & 0.000 \\ |
| 1578 |
– |
SP & 0.0 & 0.978 & 0.469 & 0.996 & 0.672 & 0.975 & 0.668 \\ |
| 1579 |
– |
& 0.1 & 0.944 & 0.645 & 0.997 & 0.886 & 0.991 & 0.978 \\ |
| 1580 |
– |
& 0.2 & 0.873 & 0.896 & 0.985 & 0.993 & 0.980 & 0.993 \\ |
| 1581 |
– |
& 0.3 & 0.831 & 0.860 & 0.960 & 0.979 & 0.955 & 0.977 \\ |
| 1582 |
– |
SF & 0.0 & 0.858 & 0.905 & 0.985 & 0.970 & 0.990 & 0.998 \\ |
| 1583 |
– |
& 0.1 & 0.865 & 0.907 & 0.992 & 0.974 & 0.994 & 0.999 \\ |
| 1584 |
– |
& 0.2 & 0.862 & 0.894 & 0.985 & 0.993 & 0.980 & 0.993 \\ |
| 1585 |
– |
& 0.3 & 0.831 & 0.859 & 0.960 & 0.979 & 0.955 & 0.977 \\ |
| 1586 |
– |
GSC & & 1.985 & 0.152 & 0.760 & 0.031 & 1.106 & 0.062 \\ |
| 1587 |
– |
RF & & 2.414 & 0.116 & 0.813 & 0.017 & 1.434 & 0.047 \\ |
| 1588 |
– |
\midrule |
| 1589 |
– |
PC & & -7.028 & 0.000 & -9.364 & 0.000 & 0.925 & 0.865 \\ |
| 1590 |
– |
SP & 0.0 & 0.701 & 0.319 & 0.909 & 0.773 & 0.861 & 0.665 \\ |
| 1591 |
– |
& 0.1 & 0.824 & 0.565 & 0.970 & 0.930 & 0.990 & 0.979 \\ |
| 1592 |
– |
& 0.2 & 0.988 & 0.981 & 0.995 & 0.998 & 0.991 & 0.998 \\ |
| 1593 |
– |
& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ |
| 1594 |
– |
SF & 0.0 & 0.993 & 0.988 & 0.992 & 0.984 & 0.998 & 0.999 \\ |
| 1595 |
– |
& 0.1 & 0.993 & 0.989 & 0.993 & 0.986 & 0.998 & 1.000 \\ |
| 1596 |
– |
& 0.2 & 0.993 & 0.992 & 0.995 & 0.998 & 0.991 & 0.998 \\ |
| 1597 |
– |
& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ |
| 1598 |
– |
GSC & & 0.964 & 0.897 & 0.970 & 0.917 & 0.925 & 0.865 \\ |
| 1599 |
– |
RF & & 0.994 & 0.864 & 0.988 & 0.865 & 0.980 & 0.784 \\ |
| 1600 |
– |
\midrule |
| 1601 |
– |
PC & & -2.212 & 0.000 & -0.588 & 0.000 & 0.953 & 0.925 \\ |
| 1602 |
– |
SP & 0.0 & 0.800 & 0.479 & 0.930 & 0.804 & 0.924 & 0.759 \\ |
| 1603 |
– |
& 0.1 & 0.883 & 0.694 & 0.976 & 0.942 & 0.993 & 0.986 \\ |
| 1604 |
– |
& 0.2 & 0.952 & 0.943 & 0.980 & 0.984 & 0.980 & 0.983 \\ |
| 1605 |
– |
& 0.3 & 0.914 & 0.909 & 0.943 & 0.948 & 0.944 & 0.946 \\ |
| 1606 |
– |
SF & 0.0 & 0.945 & 0.953 & 0.980 & 0.984 & 0.991 & 0.998 \\ |
| 1607 |
– |
& 0.1 & 0.951 & 0.954 & 0.987 & 0.986 & 0.995 & 0.998 \\ |
| 1608 |
– |
& 0.2 & 0.951 & 0.946 & 0.980 & 0.984 & 0.980 & 0.983 \\ |
| 1609 |
– |
& 0.3 & 0.914 & 0.908 & 0.943 & 0.948 & 0.944 & 0.946 \\ |
| 1610 |
– |
GSC & & 0.882 & 0.818 & 0.939 & 0.902 & 0.953 & 0.925 \\ |
| 1611 |
– |
RF & & 0.949 & 0.939 & 0.988 & 0.988 & 0.992 & 0.993 \\ |
| 1612 |
– |
\bottomrule |
| 1613 |
– |
\end{tabular} |
| 1614 |
– |
\label{tab:solnStr} |
| 1615 |
– |
\end{table} |
| 1616 |
– |
|
| 1617 |
– |
\begin{table}[htbp] |
| 1618 |
– |
\centering |
| 1619 |
– |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS |
| 1620 |
– |
OF THE FORCE AND TORQUE VECTORS IN THE STRONG SODIUM CHLORIDE SOLUTION |
| 1621 |
– |
SYSTEM} |
| 1622 |
– |
|
| 1623 |
– |
\footnotesize |
| 1624 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1625 |
– |
\toprule |
| 1626 |
– |
\toprule |
| 1627 |
– |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1628 |
– |
\cmidrule(lr){3-5} |
| 1629 |
– |
\cmidrule(l){6-8} |
| 1630 |
– |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1631 |
– |
\midrule |
| 1632 |
– |
PC & & 957.784 & 513.373 & 2.260 & 340.043 & 179.443 & 13.079 \\ |
| 1633 |
– |
SP & 0.0 & 786.244 & 139.985 & 259.289 & 311.519 & 90.280 & 105.187 \\ |
| 1634 |
– |
& 0.1 & 354.697 & 38.614 & 12.274 & 144.531 & 23.787 & 5.401 \\ |
| 1635 |
– |
& 0.2 & 7.674 & 0.363 & 0.215 & 16.655 & 3.601 & 3.634 \\ |
| 1636 |
– |
& 0.3 & 1.745 & 1.456 & 1.449 & 23.669 & 14.376 & 14.240 \\ |
| 1637 |
– |
SF & 0.0 & 3.282 & 8.567 & 0.369 & 11.904 & 6.589 & 0.717 \\ |
| 1638 |
– |
& 0.1 & 3.263 & 7.479 & 0.142 & 11.634 & 5.750 & 0.591 \\ |
| 1639 |
– |
& 0.2 & 0.686 & 0.324 & 0.215 & 10.809 & 3.580 & 3.635 \\ |
| 1640 |
– |
& 0.3 & 1.749 & 1.456 & 1.449 & 23.635 & 14.375 & 14.240 \\ |
| 1641 |
– |
GSC & & 6.181 & 2.904 & 2.263 & 44.349 & 19.442 & 12.873 \\ |
| 1642 |
– |
RF & & 3.891 & 0.847 & 0.323 & 18.628 & 3.995 & 2.072 \\ |
| 1643 |
– |
\midrule |
| 1644 |
– |
GSSP & 0.0 & 6.197 & 2.929 & 2.290 & 44.441 & 19.442 & 12.873 \\ |
| 1645 |
– |
& 0.1 & 4.688 & 1.064 & 0.260 & 31.208 & 6.967 & 2.303 \\ |
| 1646 |
– |
& 0.2 & 1.021 & 0.218 & 0.213 & 14.425 & 3.629 & 3.649 \\ |
| 1647 |
– |
& 0.3 & 1.752 & 1.454 & 1.451 & 23.540 & 14.390 & 14.245 \\ |
| 1648 |
– |
GSSF & 0.0 & 2.494 & 0.546 & 0.217 & 16.391 & 3.230 & 1.613 \\ |
| 1649 |
– |
& 0.1 & 2.448 & 0.429 & 0.106 & 16.390 & 2.827 & 1.159 \\ |
| 1650 |
– |
& 0.2 & 0.899 & 0.214 & 0.213 & 13.542 & 3.583 & 3.645 \\ |
| 1651 |
– |
& 0.3 & 1.752 & 1.454 & 1.451 & 23.587 & 14.390 & 14.245 \\ |
| 1652 |
– |
\bottomrule |
| 1653 |
– |
\end{tabular} |
| 1654 |
– |
\label{tab:solnStrAng} |
| 1655 |
– |
\end{table} |
| 1656 |
– |
|
| 1657 |
– |
The {\sc rf} method struggles with the jump in ionic strength. The |
| 1658 |
– |
configuration energy differences degrade to unusable levels while the |
| 1659 |
– |
forces and torques show a more modest reduction in the agreement with |
| 1660 |
– |
{\sc spme}. The {\sc rf} method was designed for homogeneous systems, |
| 1661 |
– |
and this attribute is apparent in these results. |
| 1662 |
– |
|
| 1663 |
– |
The {\sc sp} and {\sc sf} methods require larger cutoffs to maintain |
| 1664 |
– |
their agreement with {\sc spme}. With these results, we still |
| 1665 |
– |
recommend undamped to moderate damping for the {\sc sf} method and |
| 1666 |
– |
moderate damping for the {\sc sp} method, both with cutoffs greater |
| 1667 |
– |
than 12~\AA. |
| 1668 |
– |
|
| 1669 |
– |
\subsection{6~\AA\ Argon Sphere in SPC/E Water Results} |
| 1670 |
– |
|
| 1671 |
– |
The final model system studied was a 6~\AA\ sphere of Argon solvated |
| 1672 |
– |
by SPC/E water. This serves as a test case of a specifically sized |
| 1673 |
– |
electrostatic defect in a disordered molecular system. The results for |
| 1674 |
– |
the energy gap comparisons and the force and torque vector magnitude |
| 1675 |
– |
comparisons are shown in table \ref{tab:argon}. The force and torque |
| 1676 |
– |
vector directionality results are displayed separately in table |
| 1677 |
– |
\ref{tab:argonAng}, where the effect of group-based cutoffs and |
| 1678 |
– |
switching functions on the {\sc sp} and {\sc sf} potentials are |
| 1679 |
– |
investigated. |
| 1680 |
– |
|
| 1681 |
– |
\begin{table}[htbp] |
| 1682 |
– |
\centering |
| 1683 |
– |
\caption{REGRESSION RESULTS OF THE 6~\AA\ ARGON SPHERE IN LIQUID |
| 1684 |
– |
WATER SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR |
| 1685 |
– |
MAGNITUDES ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1686 |
– |
|
| 1687 |
– |
\footnotesize |
| 1688 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1689 |
– |
\toprule |
| 1690 |
– |
\toprule |
| 1691 |
– |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1692 |
– |
\cmidrule(lr){3-4} |
| 1693 |
– |
\cmidrule(lr){5-6} |
| 1694 |
– |
\cmidrule(l){7-8} |
| 1695 |
– |
Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ |
| 1696 |
– |
\midrule |
| 1697 |
– |
PC & & 2.320 & 0.008 & -0.650 & 0.001 & 3.848 & 0.029 \\ |
| 1698 |
– |
SP & 0.0 & 1.053 & 0.711 & 0.977 & 0.820 & 0.974 & 0.882 \\ |
| 1699 |
– |
& 0.1 & 1.032 & 0.846 & 0.989 & 0.965 & 0.992 & 0.994 \\ |
| 1700 |
– |
& 0.2 & 0.993 & 0.995 & 0.982 & 0.998 & 0.986 & 0.998 \\ |
| 1701 |
– |
& 0.3 & 0.968 & 0.995 & 0.954 & 0.992 & 0.961 & 0.994 \\ |
| 1702 |
– |
SF & 0.0 & 0.982 & 0.996 & 0.992 & 0.999 & 0.993 & 1.000 \\ |
| 1703 |
– |
& 0.1 & 0.987 & 0.996 & 0.996 & 0.999 & 0.997 & 1.000 \\ |
| 1704 |
– |
& 0.2 & 0.989 & 0.998 & 0.984 & 0.998 & 0.989 & 0.998 \\ |
| 1705 |
– |
& 0.3 & 0.971 & 0.995 & 0.957 & 0.992 & 0.965 & 0.994 \\ |
| 1706 |
– |
GSC & & 1.002 & 0.983 & 0.992 & 0.973 & 0.996 & 0.971 \\ |
| 1707 |
– |
RF & & 0.998 & 0.995 & 0.999 & 0.998 & 0.998 & 0.998 \\ |
| 1708 |
– |
\midrule |
| 1709 |
– |
PC & & -36.559 & 0.002 & -44.917 & 0.004 & -52.945 & 0.006 \\ |
| 1710 |
– |
SP & 0.0 & 0.890 & 0.786 & 0.927 & 0.867 & 0.949 & 0.909 \\ |
| 1711 |
– |
& 0.1 & 0.942 & 0.895 & 0.984 & 0.974 & 0.997 & 0.995 \\ |
| 1712 |
– |
& 0.2 & 0.999 & 0.997 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1713 |
– |
& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ |
| 1714 |
– |
SF & 0.0 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1715 |
– |
& 0.1 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1716 |
– |
& 0.2 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1717 |
– |
& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ |
| 1718 |
– |
GSC & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1719 |
– |
RF & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ |
| 1720 |
– |
\midrule |
| 1721 |
– |
PC & & 1.984 & 0.000 & 0.012 & 0.000 & 1.357 & 0.000 \\ |
| 1722 |
– |
SP & 0.0 & 0.850 & 0.552 & 0.907 & 0.703 & 0.938 & 0.793 \\ |
| 1723 |
– |
& 0.1 & 0.924 & 0.755 & 0.980 & 0.936 & 0.995 & 0.988 \\ |
| 1724 |
– |
& 0.2 & 0.985 & 0.983 & 0.986 & 0.988 & 0.987 & 0.988 \\ |
| 1725 |
– |
& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ |
| 1726 |
– |
SF & 0.0 & 0.977 & 0.989 & 0.987 & 0.995 & 0.992 & 0.998 \\ |
| 1727 |
– |
& 0.1 & 0.982 & 0.989 & 0.992 & 0.996 & 0.997 & 0.998 \\ |
| 1728 |
– |
& 0.2 & 0.984 & 0.987 & 0.986 & 0.987 & 0.987 & 0.988 \\ |
| 1729 |
– |
& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ |
| 1730 |
– |
GSC & & 0.995 & 0.981 & 0.999 & 0.990 & 1.000 & 0.993 \\ |
| 1731 |
– |
RF & & 0.993 & 0.988 & 0.997 & 0.995 & 0.999 & 0.998 \\ |
| 1732 |
– |
\bottomrule |
| 1733 |
– |
\end{tabular} |
| 1734 |
– |
\label{tab:argon} |
| 1735 |
– |
\end{table} |
| 1736 |
– |
|
| 1737 |
– |
\begin{table}[htbp] |
| 1738 |
– |
\centering |
| 1739 |
– |
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1740 |
– |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6~\AA\ SPHERE OF |
| 1741 |
– |
ARGON IN LIQUID WATER SYSTEM} |
| 1742 |
– |
|
| 1743 |
– |
\footnotesize |
| 1744 |
– |
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1745 |
– |
\toprule |
| 1746 |
– |
\toprule |
| 1747 |
– |
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1748 |
– |
\cmidrule(lr){3-5} |
| 1749 |
– |
\cmidrule(l){6-8} |
| 1750 |
– |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1751 |
– |
\midrule |
| 1752 |
– |
PC & & 568.025 & 265.993 & 195.099 & 246.626 & 138.600 & 91.654 \\ |
| 1753 |
– |
SP & 0.0 & 504.578 & 251.694 & 179.932 & 231.568 & 131.444 & 85.119 \\ |
| 1754 |
– |
& 0.1 & 224.886 & 49.746 & 9.346 & 104.482 & 23.683 & 4.480 \\ |
| 1755 |
– |
& 0.2 & 4.889 & 0.197 & 0.155 & 6.029 & 2.507 & 2.269 \\ |
| 1756 |
– |
& 0.3 & 0.817 & 0.833 & 0.812 & 8.286 & 8.436 & 8.135 \\ |
| 1757 |
– |
SF & 0.0 & 1.924 & 0.675 & 0.304 & 3.658 & 1.448 & 0.600 \\ |
| 1758 |
– |
& 0.1 & 1.937 & 0.515 & 0.143 & 3.565 & 1.308 & 0.546 \\ |
| 1759 |
– |
& 0.2 & 0.407 & 0.166 & 0.156 & 3.086 & 2.501 & 2.274 \\ |
| 1760 |
– |
& 0.3 & 0.815 & 0.833 & 0.812 & 8.330 & 8.437 & 8.135 \\ |
| 1761 |
– |
GSC & & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ |
| 1762 |
– |
RF & & 1.822 & 0.408 & 0.142 & 3.799 & 1.362 & 0.550 \\ |
| 1763 |
– |
\midrule |
| 1764 |
– |
GSSP & 0.0 & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ |
| 1765 |
– |
& 0.1 & 1.652 & 0.309 & 0.087 & 4.197 & 1.401 & 0.590 \\ |
| 1766 |
– |
& 0.2 & 0.465 & 0.165 & 0.153 & 3.323 & 2.529 & 2.273 \\ |
| 1767 |
– |
& 0.3 & 0.813 & 0.825 & 0.816 & 8.316 & 8.447 & 8.132 \\ |
| 1768 |
– |
GSSF & 0.0 & 1.173 & 0.292 & 0.113 & 3.452 & 1.347 & 0.583 \\ |
| 1769 |
– |
& 0.1 & 1.166 & 0.240 & 0.076 & 3.381 & 1.281 & 0.575 \\ |
| 1770 |
– |
& 0.2 & 0.459 & 0.165 & 0.153 & 3.430 & 2.542 & 2.273 \\ |
| 1771 |
– |
& 0.3 & 0.814 & 0.825 & 0.816 & 8.325 & 8.447 & 8.132 \\ |
| 1772 |
– |
\bottomrule |
| 1773 |
– |
\end{tabular} |
| 1774 |
– |
\label{tab:argonAng} |
| 1775 |
– |
\end{table} |
| 1776 |
– |
|
| 1777 |
– |
This system does not appear to show any significant deviations from |
| 1778 |
– |
the previously observed results. The {\sc sp} and {\sc sf} methods |
| 1779 |
– |
have agreements similar to those observed in section |
| 1780 |
– |
\ref{sec:WaterResults}. The only significant difference is the |
| 1781 |
– |
improvement in the configuration energy differences for the {\sc rf} |
| 1782 |
– |
method. This is surprising in that we are introducing an inhomogeneity |
| 1783 |
– |
to the system; however, this inhomogeneity is charge-neutral and does |
| 1784 |
– |
not result in charged cutoff spheres. The charge-neutrality of the |
| 1785 |
– |
cutoff spheres, which the {\sc sp} and {\sc sf} methods explicitly |
| 1786 |
– |
enforce, seems to play a greater role in the stability of the {\sc rf} |
| 1787 |
– |
method than the required homogeneity of the environment. |
| 1788 |
– |
|
| 1789 |
– |
|
| 938 |
|
\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics} |
| 939 |
|
|
| 940 |
|
Zahn {\it et al.} investigated the structure and dynamics of water |
| 941 |
|
using equations (\ref{eq:ZahnPot}) and |
| 942 |
|
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 943 |
|
that a method similar (but not identical with) the damped {\sc sf} |
| 944 |
< |
method resulted in properties very similar to those obtained when |
| 944 |
> |
method resulted in properties very close to those obtained when |
| 945 |
|
using the Ewald summation. The properties they studied (pair |
| 946 |
|
distribution functions, diffusion constants, and velocity and |
| 947 |
|
orientational correlation functions) may not be particularly sensitive |
| 948 |
|
to the long-range and collective behavior that governs the |
| 949 |
|
low-frequency behavior in crystalline systems. Additionally, the |
| 950 |
< |
ionic crystals are the worst case scenario for the pairwise methods |
| 950 |
> |
ionic crystals are a worst case scenario for the pairwise methods |
| 951 |
|
because they lack the reciprocal space contribution contained in the |
| 952 |
|
Ewald summation. |
| 953 |
|
|
| 954 |
< |
We are using two separate measures to probe the effects of these |
| 954 |
> |
We used two separate measures to probe the effects of these |
| 955 |
|
alternative electrostatic methods on the dynamics in crystalline |
| 956 |
< |
materials. For short- and intermediate-time dynamics, we are |
| 957 |
< |
computing the velocity autocorrelation function, and for long-time |
| 958 |
< |
and large length-scale collective motions, we are looking at the |
| 959 |
< |
low-frequency portion of the power spectrum. |
| 956 |
> |
materials. For short- and intermediate-time dynamics, we computed the |
| 957 |
> |
velocity autocorrelation function, and for long-time and large |
| 958 |
> |
length-scale collective motions, we looked at the low-frequency |
| 959 |
> |
portion of the power spectrum. |
| 960 |
|
|
| 961 |
|
\begin{figure} |
| 962 |
|
\centering |
| 970 |
|
are stiffer than the moderately damped and {\sc spme} methods.} |
| 971 |
|
\label{fig:vCorrPlot} |
| 972 |
|
\end{figure} |
| 973 |
< |
|
| 1826 |
< |
The short-time decay of the velocity autocorrelation function through |
| 973 |
> |
The short-time decay of the velocity autocorrelation functions through |
| 974 |
|
the first collision are nearly identical in figure |
| 975 |
|
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 976 |
|
how the methods differ. The undamped {\sc sf} method has deeper |
| 977 |
< |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 977 |
> |
troughs (see inset in figure \ref{fig:vCorrPlot}) and higher peaks than |
| 978 |
|
any of the other methods. As the damping parameter ($\alpha$) is |
| 979 |
|
increased, these peaks are smoothed out, and the {\sc sf} method |
| 980 |
|
approaches the {\sc spme} results. With $\alpha$ values of 0.2~\AA$^{-1}$, |
| 988 |
|
|
| 989 |
|
\section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics} |
| 990 |
|
|
| 1844 |
– |
To evaluate how the differences between the methods affect the |
| 1845 |
– |
collective long-time motion, we computed power spectra from long-time |
| 1846 |
– |
traces of the velocity autocorrelation function. The power spectra for |
| 1847 |
– |
the best-performing alternative methods are shown in |
| 1848 |
– |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1849 |
– |
a cubic switching function between 40 and 50~ps was used to reduce the |
| 1850 |
– |
ringing resulting from data truncation. This procedure had no |
| 1851 |
– |
noticeable effect on peak location or magnitude. |
| 1852 |
– |
|
| 991 |
|
\begin{figure} |
| 992 |
|
\centering |
| 993 |
|
\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf} |
| 998 |
|
100~cm$^{-1}$ to highlight where the spectra differ.} |
| 999 |
|
\label{fig:methodPS} |
| 1000 |
|
\end{figure} |
| 1001 |
+ |
To evaluate how the differences between the methods affect the |
| 1002 |
+ |
collective long-time motion, we computed power spectra from long-time |
| 1003 |
+ |
traces of the velocity autocorrelation function. The power spectra for |
| 1004 |
+ |
the best-performing alternative methods are shown in |
| 1005 |
+ |
figure \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1006 |
+ |
a cubic switching function between 40 and 50~ps was used to reduce the |
| 1007 |
+ |
ringing resulting from data truncation. This procedure had no |
| 1008 |
+ |
noticeable effect on peak location or magnitude. |
| 1009 |
|
|
| 1010 |
|
While the high frequency regions of the power spectra for the |
| 1011 |
|
alternative methods are quantitatively identical with Ewald spectrum, |
| 1012 |
|
the low frequency region shows how the summation methods differ. |
| 1013 |
|
Considering the low-frequency inset (expanded in the upper frame of |
| 1014 |
< |
figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the |
| 1014 |
> |
figure \ref{fig:methodPS}), at frequencies below 100~cm$^{-1}$, the |
| 1015 |
|
correlated motions are blue-shifted when using undamped or weakly |
| 1016 |
|
damped {\sc sf}. When using moderate damping ($\alpha = |
| 1017 |
< |
0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly |
| 1018 |
< |
identical correlated motion to the Ewald method (which has a |
| 1017 |
> |
0.2$~\AA$^{-1}$), both the {\sc sf} and {\sc sp} methods produce |
| 1018 |
> |
correlated motions nearly identical to the Ewald method (which has a |
| 1019 |
|
convergence parameter of 0.3119~\AA$^{-1}$). This weakening of the |
| 1020 |
|
electrostatic interaction with increased damping explains why the |
| 1021 |
|
long-ranged correlated motions are at lower frequencies for the |
| 1022 |
|
moderately damped methods than for undamped or weakly damped methods. |
| 1023 |
|
|
| 1878 |
– |
To isolate the role of the damping constant, we have computed the |
| 1879 |
– |
spectra for a single method ({\sc sf}) with a range of damping |
| 1880 |
– |
constants and compared this with the {\sc spme} spectrum. |
| 1881 |
– |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1882 |
– |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1883 |
– |
However, even without any electrostatic damping, the {\sc sf} method |
| 1884 |
– |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1885 |
– |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1886 |
– |
would predict the lowest frequency peak near 325~cm$^{-1}$. {\it |
| 1887 |
– |
Most} of the collective behavior in the crystal is accurately captured |
| 1888 |
– |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1889 |
– |
obtained using moderate damping in addition to the shifting at the |
| 1890 |
– |
cutoff distance. |
| 1891 |
– |
|
| 1024 |
|
\begin{figure} |
| 1025 |
|
\centering |
| 1026 |
|
\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf} |
| 1033 |
|
motions.} |
| 1034 |
|
\label{fig:dampInc} |
| 1035 |
|
\end{figure} |
| 1036 |
+ |
To isolate the role of the damping constant, we have computed the |
| 1037 |
+ |
spectra for a single method ({\sc sf}) with a range of damping |
| 1038 |
+ |
constants and compared this with the {\sc spme} spectrum. Figure |
| 1039 |
+ |
\ref{fig:dampInc} shows more clearly that increasing the electrostatic |
| 1040 |
+ |
damping red-shifts the lowest frequency phonon modes. However, even |
| 1041 |
+ |
without any electrostatic damping, the {\sc sf} method has at most a |
| 1042 |
+ |
10 cm$^{-1}$ error in the lowest frequency phonon mode. Without the |
| 1043 |
+ |
{\sc sf} modifications, an undamped (pure cutoff) method would predict |
| 1044 |
+ |
the lowest frequency peak near 325~cm$^{-1}$, an error significantly |
| 1045 |
+ |
larger than that of the undamped {\sc sf} technique. This indicates |
| 1046 |
+ |
that {\it most} of the collective behavior in the crystal is |
| 1047 |
+ |
accurately captured using the {\sc sf} method. Quantitative agreement |
| 1048 |
+ |
with Ewald can be obtained using moderate damping in addition to the |
| 1049 |
+ |
shifting at the cutoff distance. |
| 1050 |
|
|
| 1051 |
|
\section{An Application: TIP5P-E Water}\label{sec:t5peApplied} |
| 1052 |
|
|
| 1067 |
|
long-range electrostatic |
| 1068 |
|
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
| 1069 |
|
Without this correction, the pressure term on the central particle |
| 1070 |
< |
from the surroundings is missing. Because they expand to compensate |
| 1071 |
< |
for this added pressure term when this correction is included, systems |
| 1072 |
< |
composed of these particles tend to under-predict the density of water |
| 1073 |
< |
under standard conditions. When using any form of long-range |
| 1074 |
< |
electrostatic correction, it has become common practice to develop or |
| 1075 |
< |
utilize a reparametrized water model that corrects for this |
| 1070 |
> |
from the surroundings is missing. When this correction is included, |
| 1071 |
> |
systems of these particles expand to compensate for this added |
| 1072 |
> |
pressure term and under-predict the density of water under standard |
| 1073 |
> |
conditions. When using any form of long-range electrostatic |
| 1074 |
> |
correction, it has become common practice to develop or utilize a |
| 1075 |
> |
reparametrized water model that corrects for this |
| 1076 |
|
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 1077 |
< |
this practice and was optimized specifically for use with the Ewald |
| 1077 |
> |
this practice and was optimized for use with the Ewald |
| 1078 |
|
summation.\cite{Rick04} In his publication, Rick preserved the |
| 1079 |
|
geometry and point charge magnitudes in TIP5P and focused on altering |
| 1080 |
< |
the Lennard-Jones parameters to correct the density at |
| 1081 |
< |
298K.\cite{Rick04} With the density corrected, he compared common |
| 1082 |
< |
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1937 |
< |
9~\AA\ cutoff. |
| 1080 |
> |
the Lennard-Jones parameters to correct the density at 298~K. With the |
| 1081 |
> |
density corrected, he compared common water properties for TIP5P-E |
| 1082 |
> |
using the Ewald sum with TIP5P using a 9~\AA\ cutoff. |
| 1083 |
|
|
| 1084 |
|
In the following sections, we compared these same water properties |
| 1085 |
|
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1104 |
|
+ \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right], |
| 1105 |
|
\label{eq:MolecularPressure} |
| 1106 |
|
\end{equation} |
| 1107 |
< |
where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of |
| 1108 |
< |
molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass |
| 1109 |
< |
($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on |
| 1110 |
< |
atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the |
| 1111 |
< |
right term in the brackets of equation \ref{eq:MolecularPressure}) is |
| 1112 |
< |
directly dependent on the interatomic forces. Since the {\sc sp} |
| 1113 |
< |
method does not modify the forces (see |
| 1114 |
< |
section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will |
| 1115 |
< |
be identical to that obtained without an electrostatic correction. |
| 1116 |
< |
The {\sc sf} method does alter the virial component and, by way of the |
| 1117 |
< |
modified pressures, should provide densities more in line with those |
| 1118 |
< |
obtained using the Ewald summation. |
| 1107 |
> |
where d is the dimensionality of the system, $V$ is the volume, |
| 1108 |
> |
$\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$ |
| 1109 |
> |
is the position of the center of mass ($M_\mu$) of molecule $\mu$, and |
| 1110 |
> |
$\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule |
| 1111 |
> |
$\mu$.\cite{Melchionna93} The virial term (the right term in the |
| 1112 |
> |
brackets of equation |
| 1113 |
> |
\ref{eq:MolecularPressure}) is directly dependent on the interatomic |
| 1114 |
> |
forces. Since the {\sc sp} method does not modify the forces (see |
| 1115 |
> |
section \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} |
| 1116 |
> |
will be identical to that obtained without an electrostatic |
| 1117 |
> |
correction. The {\sc sf} method does alter the virial component and, |
| 1118 |
> |
by way of the modified pressures, should provide densities more in |
| 1119 |
> |
line with those obtained using the Ewald summation. |
| 1120 |
|
|
| 1121 |
|
To compare densities, $NPT$ simulations were performed with the same |
| 1122 |
|
temperatures as those selected by Rick in his Ewald summation |
| 1126 |
|
temperatures. The average densities were calculated from the later |
| 1127 |
|
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 1128 |
|
method for accumulating statistics, these sequences were spliced into |
| 1129 |
< |
200 segments to calculate the average density and standard deviation |
| 1130 |
< |
at each temperature.\cite{Mahoney00} |
| 1129 |
> |
200 segments, each providing an average density. These 200 density |
| 1130 |
> |
values were used to calculate the average and standard deviation of |
| 1131 |
> |
the density at each temperature.\cite{Mahoney00} |
| 1132 |
|
|
| 1133 |
|
\begin{figure} |
| 1134 |
|
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 1139 |
|
Ewald summation, leading to slightly lower densities. This effect is |
| 1140 |
|
more visible with the 9~\AA\ cutoff, where the image charges exert a |
| 1141 |
|
greater force on the central particle. The error bars for the {\sc sf} |
| 1142 |
< |
methods show plus or minus the standard deviation of the density |
| 1143 |
< |
measurement at each temperature.} |
| 1142 |
> |
methods show the average one-sigma uncertainty of the density |
| 1143 |
> |
measurement, and this uncertainty is the same for all the {\sc sf} |
| 1144 |
> |
curves.} |
| 1145 |
|
\label{fig:t5peDensities} |
| 1146 |
|
\end{figure} |
| 1999 |
– |
|
| 1147 |
|
Figure \ref{fig:t5peDensities} shows the densities calculated for |
| 1148 |
|
TIP5P-E using differing electrostatic corrections overlaid on the |
| 1149 |
|
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 1150 |
|
technique are close to, though typically lower than, those calculated |
| 1151 |
< |
while using the Ewald summation. These slightly reduced densities |
| 1152 |
< |
indicate that the pressure component from the image charges at |
| 1153 |
< |
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 1154 |
< |
portion of the Ewald summation. Bringing the image charges closer to |
| 1155 |
< |
the central particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than |
| 1156 |
< |
the preferred 12~\AA\ R$_\textrm{c}$) increases the strength of their |
| 1157 |
< |
interactions, resulting in a further reduction of the densities. |
| 1151 |
> |
using the Ewald summation. These slightly reduced densities indicate |
| 1152 |
> |
that the pressure component from the image charges at R$_\textrm{c}$ |
| 1153 |
> |
is larger than that exerted by the reciprocal-space portion of the |
| 1154 |
> |
Ewald summation. Bringing the image charges closer to the central |
| 1155 |
> |
particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the |
| 1156 |
> |
preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image |
| 1157 |
> |
charge interactions on the central particle and results in a further |
| 1158 |
> |
reduction of the densities. |
| 1159 |
|
|
| 1160 |
|
Because the strength of the image charge interactions has a noticeable |
| 1161 |
|
effect on the density, we would expect the use of electrostatic |
| 1169 |
|
important role in the resulting densities. |
| 1170 |
|
|
| 1171 |
|
As a final note, all of the above density calculations were performed |
| 1172 |
< |
with systems of 512 water molecules. Rick observed a system sized |
| 1172 |
> |
with systems of 512 water molecules. Rick observed a system size |
| 1173 |
|
dependence of the computed densities when using the Ewald summation, |
| 1174 |
|
most likely due to his tying of the convergence parameter to the box |
| 1175 |
|
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
| 1211 |
|
non-polarizable models.\cite{Sorenson00} This excellent agreement with |
| 1212 |
|
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To |
| 1213 |
|
check whether the choice of using the Ewald summation or the {\sc sf} |
| 1214 |
< |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K |
| 1215 |
< |
and 1atm were determined for the systems compared in the previous |
| 1216 |
< |
section. |
| 1214 |
> |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at |
| 1215 |
> |
298~K and 1~atm were determined for the systems compared in the |
| 1216 |
> |
previous section. |
| 1217 |
|
|
| 1218 |
|
\begin{figure} |
| 1219 |
|
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 1224 |
|
identical.} |
| 1225 |
|
\label{fig:t5peGofRs} |
| 1226 |
|
\end{figure} |
| 2079 |
– |
|
| 1227 |
|
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc |
| 1228 |
|
sf} technique with a various parameters are overlaid on the |
| 1229 |
< |
$g_\textrm{OO}(r)$ while using the Ewald summation. The differences in |
| 1230 |
< |
density do not appear to have any effect on the liquid structure as |
| 1231 |
< |
the $g_\textrm{OO}(r)$s are indistinguishable. These results indicate |
| 1232 |
< |
that the $g_\textrm{OO}(r)$ is insensitive to the choice of |
| 1233 |
< |
electrostatic correction. |
| 1229 |
> |
$g_\textrm{OO}(r)$ while using the Ewald summation in |
| 1230 |
> |
figure~\ref{fig:t5peGofRs}. The differences in density do not appear |
| 1231 |
> |
to have any effect on the liquid structure as the $g_\textrm{OO}(r)$s |
| 1232 |
> |
are indistinguishable. These results indicate that the |
| 1233 |
> |
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic |
| 1234 |
> |
correction. |
| 1235 |
|
|
| 1236 |
|
\subsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 1237 |
|
|
| 1246 |
|
good set for comparisons involving the {\sc sf} technique. |
| 1247 |
|
|
| 1248 |
|
The $\Delta H_\textrm{vap}$ is the enthalpy change required to |
| 1249 |
< |
transform one mol of substance from the liquid phase to the gas |
| 1249 |
> |
transform one mole of substance from the liquid phase to the gas |
| 1250 |
|
phase.\cite{Berry00} In molecular simulations, this quantity can be |
| 1251 |
|
determined via |
| 1252 |
|
\begin{equation} |
| 1361 |
|
|
| 1362 |
|
As observed for the density in section \ref{sec:t5peDensity}, the |
| 1363 |
|
property trends with temperature seen when using the Ewald summation |
| 1364 |
< |
are reproduced with the {\sc sf} technique. Differences include the |
| 1365 |
< |
calculated values of $\Delta H_\textrm{vap}$ under-predicting the Ewald |
| 1366 |
< |
values. This is to be expected due to the direct weakening of the |
| 1367 |
< |
electrostatic interaction through forced neutralization in {\sc |
| 1368 |
< |
sf}. This results in an increase of the intermolecular potential |
| 1369 |
< |
producing lower values from equation (\ref{eq:DeltaHVap}). The slopes of |
| 1370 |
< |
these values with temperature are similar to that seen using the Ewald |
| 1371 |
< |
summation; however, they are both steeper than the experimental trend, |
| 1372 |
< |
indirectly resulting in the inflated $C_p$ values at all temperatures. |
| 1364 |
> |
are reproduced with the {\sc sf} technique. One noticeable difference |
| 1365 |
> |
between the properties calculated using the two methods are the lower |
| 1366 |
> |
$\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be |
| 1367 |
> |
expected due to the direct weakening of the electrostatic interaction |
| 1368 |
> |
through forced neutralization. This results in an increase of the |
| 1369 |
> |
intermolecular potential producing lower values from equation |
| 1370 |
> |
(\ref{eq:DeltaHVap}). The slopes of these values with temperature are |
| 1371 |
> |
similar to that seen using the Ewald summation; however, they are both |
| 1372 |
> |
steeper than the experimental trend, indirectly resulting in the |
| 1373 |
> |
inflated $C_p$ values at all temperatures. |
| 1374 |
|
|
| 1375 |
< |
Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ |
| 1376 |
< |
values all overlap within error. As indicated for the $\Delta |
| 1377 |
< |
H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph, |
| 1378 |
< |
the deviations between experiment and simulation in this region are |
| 1379 |
< |
not the fault of the electrostatic summation methods but are due to |
| 1380 |
< |
the TIP5P class model itself. Like most rigid, non-polarizable, |
| 1381 |
< |
point-charge water models, the density decreases with temperature at a |
| 1382 |
< |
much faster rate than experiment (see figure |
| 1383 |
< |
\ref{fig:t5peDensities}). The reduced density leads to the inflated |
| 1375 |
> |
Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values |
| 1376 |
> |
all overlap within error. As indicated for the $\Delta H_\textrm{vap}$ |
| 1377 |
> |
and $C_p$ results discussed in the previous paragraph, the deviations |
| 1378 |
> |
between experiment and simulation in this region are not the fault of |
| 1379 |
> |
the electrostatic summation methods but are due to the geometry and |
| 1380 |
> |
parameters of the TIP5P class of water models. Like most rigid, |
| 1381 |
> |
non-polarizable, point-charge water models, the density decreases with |
| 1382 |
> |
temperature at a much faster rate than experiment (see figure |
| 1383 |
> |
\ref{fig:t5peDensities}). This reduced density leads to the inflated |
| 1384 |
|
compressibility and expansivity values at higher temperatures seen |
| 1385 |
|
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability |
| 1386 |
< |
and many-body effects are required in order for simulation to overcome |
| 1387 |
< |
these differences with experiment.\cite{Laasonen93,Donchev06} |
| 1386 |
> |
and many-body effects are required in order for water models to |
| 1387 |
> |
overcome differences between simulation-based and experimentally |
| 1388 |
> |
determined densities at these higher |
| 1389 |
> |
temperatures.\cite{Laasonen93,Donchev06} |
| 1390 |
|
|
| 1391 |
|
At temperatures below the freezing point for experimental water, the |
| 1392 |
|
differences between {\sc sf} and the Ewald summation results are more |
| 1394 |
|
indicate a more pronounced transition in the supercooled regime, |
| 1395 |
|
particularly in the case of {\sc sf} without damping. This points to |
| 1396 |
|
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 1397 |
< |
temperatures below 250~K in these simulations. Because the systems are |
| 1398 |
< |
locked in different regions of phase-space, comparisons between |
| 1399 |
< |
properties at these temperatures are not exactly fair. This |
| 1400 |
< |
observation is explored in more detail in section |
| 1401 |
< |
\ref{sec:t5peDynamics}. |
| 1397 |
> |
temperatures below 250~K in the {\sc sf} simulations, indicating that |
| 1398 |
> |
disorder in the reciprocal-space term of the Ewald summation might act |
| 1399 |
> |
to loosen up the local structure more than the image-charges in {\sc |
| 1400 |
> |
sf}. The damped {\sc sf} actually makes a better comparison with |
| 1401 |
> |
experiment in this region, particularly for the $\alpha_p$ values. The |
| 1402 |
> |
local interactions in the undamped {\sc sf} technique appear to be too |
| 1403 |
> |
strong since the property change is much more dramatic than the damped |
| 1404 |
> |
forms, while the Ewald summation appears to weight the |
| 1405 |
> |
reciprocal-space interactions at the expense the local interactions, |
| 1406 |
> |
disagreeing with the experimental results. This observation is |
| 1407 |
> |
explored in more detail in section \ref{sec:t5peDynamics}. |
| 1408 |
|
|
| 1409 |
|
The final thermodynamic property displayed in figure |
| 1410 |
|
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
| 1414 |
|
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the |
| 1415 |
|
converged $\epsilon$ values accumulated for the {\sc sf} |
| 1416 |
|
simulations. Lack of a damping function results in dielectric |
| 1417 |
< |
constants significantly smaller than that obtained using the Ewald |
| 1417 |
> |
constants significantly smaller than those obtained using the Ewald |
| 1418 |
|
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the |
| 1419 |
|
agreement considerably. It should be noted that the choice of the |
| 1420 |
|
``Ewald coefficient'' value also has a significant effect on the |
| 1426 |
|
sf}; however, the choice of cutoff radius also plays an important |
| 1427 |
|
role. In section \ref{sec:dampingDielectric}, this connection is |
| 1428 |
|
further explored as optimal damping coefficients for different choices |
| 1429 |
< |
of $R_\textrm{c}$ are determined for {\sc sf} for capturing the |
| 1430 |
< |
dielectric behavior. |
| 1429 |
> |
of $R_\textrm{c}$ are determined for {\sc sf} in order to best capture |
| 1430 |
> |
the dielectric behavior. |
| 1431 |
|
|
| 1432 |
|
\subsection{Dynamic Properties}\label{sec:t5peDynamics} |
| 1433 |
|
|
| 1434 |
|
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 1435 |
< |
method, 200~ps $NVE$ simulations were performed for each temperature at |
| 1436 |
< |
the average density reported by the $NPT$ simulations. The |
| 1437 |
< |
self-diffusion constants ($D$) were calculated with the Einstein |
| 1438 |
< |
relation using the mean square displacement (MSD), |
| 1435 |
> |
method, 200~ps $NVE$ simulations were performed for each temperature |
| 1436 |
> |
at the average density reported by the $NPT$ simulations. The |
| 1437 |
> |
self-diffusion constants ($D$) were calculated using the mean square |
| 1438 |
> |
displacement (MSD) form of the Einstein relation, |
| 1439 |
|
\begin{equation} |
| 1440 |
< |
D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t}, |
| 1440 |
> |
D = \lim_{t\rightarrow\infty} |
| 1441 |
> |
\frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t}, |
| 1442 |
|
\label{eq:MSD} |
| 1443 |
|
\end{equation} |
| 1444 |
|
where $t$ is time, and $\mathbf{r}_i$ is the position of particle |
| 1449 |
|
\begin{enumerate}[itemsep=0pt] |
| 1450 |
|
\item parabolic short-time ballistic motion, |
| 1451 |
|
\item linear diffusive regime, and |
| 1452 |
< |
\item poor statistic region at long-time. |
| 1452 |
> |
\item a region with poor statistics. |
| 1453 |
|
\end{enumerate} |
| 1454 |
< |
The slope from the linear region (region 2) is used to calculate $D$. |
| 1454 |
> |
The slope from the linear regime (region 2) is used to calculate $D$. |
| 1455 |
|
\begin{figure} |
| 1456 |
|
\centering |
| 1457 |
|
\includegraphics[width=3.5in]{./figures/ExampleMSD.pdf} |
| 1470 |
|
labeled frame axes.} |
| 1471 |
|
\label{fig:waterFrame} |
| 1472 |
|
\end{figure} |
| 1473 |
< |
In addition to translational diffusion, reorientational time constants |
| 1473 |
> |
In addition to translational diffusion, orientational relaxation times |
| 1474 |
|
were calculated for comparisons with the Ewald simulations and with |
| 1475 |
< |
experiments. These values were determined from 25~ps $NVE$ trajectories |
| 1476 |
< |
through calculation of the orientational time correlation function, |
| 1475 |
> |
experiments. These values were determined from 25~ps $NVE$ |
| 1476 |
> |
trajectories through calculation of the orientational time correlation |
| 1477 |
> |
function, |
| 1478 |
|
\begin{equation} |
| 1479 |
|
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t) |
| 1480 |
|
\cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle, |
| 1519 |
|
easier comparisons in the more relevant temperature regime.} |
| 1520 |
|
\label{fig:t5peDynamics} |
| 1521 |
|
\end{figure} |
| 1522 |
< |
Results for the diffusion constants and reorientational time constants |
| 1522 |
> |
Results for the diffusion constants and orientational relaxation times |
| 1523 |
|
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is |
| 1524 |
|
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using |
| 1525 |
|
the Ewald sum are reproduced with the {\sc sf} technique. The enhanced |
| 1528 |
|
insight into differences between the electrostatic summation |
| 1529 |
|
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to |
| 1530 |
|
diffuse a little faster than with the Ewald sum; however, use of light |
| 1531 |
< |
to moderate damping results in indistinguishable $D$ values. Though not |
| 1532 |
< |
apparent in this figure, {\sc sf} values at the lowest temperature are |
| 1533 |
< |
approximately an order of magnitude lower than with Ewald. These |
| 1531 |
> |
to moderate damping results in indistinguishable $D$ values. Though |
| 1532 |
> |
not apparent in this figure, {\sc sf} values at the lowest temperature |
| 1533 |
> |
are approximately an order of magnitude lower than with Ewald. These |
| 1534 |
|
values support the observation from section \ref{sec:t5peThermo} that |
| 1535 |
|
there appeared to be a change to a more glassy-like phase with the |
| 1536 |
|
{\sc sf} technique at these lower temperatures. |
| 1541 |
|
relaxes faster than experiment with the Ewald sum while tracking |
| 1542 |
|
experiment fairly well when using the {\sc sf} technique, independent |
| 1543 |
|
of the choice of damping constant. Their are several possible reasons |
| 1544 |
< |
for this deviation between techniques. The Ewald results were taken |
| 1545 |
< |
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick |
| 1546 |
< |
calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of |
| 1547 |
< |
0.2~\AA$^-1$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, placing |
| 1548 |
< |
the result more in line with that obtained using the Ewald sum. These |
| 1549 |
< |
results support this explanation; however, recomputing the results to |
| 1550 |
< |
meet a poorer statistical standard is counter-productive. Assuming the |
| 1551 |
< |
Ewald results are not the product of poor statistics, differences in |
| 1552 |
< |
techniques to integrate the orientational motion could also play a |
| 1553 |
< |
role. {\sc shake} is the most commonly used technique for |
| 1554 |
< |
approximating rigid-body orientational motion,\cite{Ryckaert77} where |
| 1555 |
< |
as in {\sc oopse}, we maintain and integrate the entire rotation |
| 1556 |
< |
matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake} |
| 1557 |
< |
is an iterative constraint technique, if the convergence tolerances |
| 1558 |
< |
are raised for increased performance, error will accumulate in the |
| 1559 |
< |
orientational motion. Finally, the Ewald results were calculated using |
| 1560 |
< |
the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf} |
| 1544 |
> |
for this deviation between techniques. The Ewald results were |
| 1545 |
> |
calculated using shorter (10ps) trajectories than the {\sc sf} results |
| 1546 |
> |
(25ps). A quick calculation from a 10~ps trajectory with {\sc sf} with |
| 1547 |
> |
an $\alpha$ of 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in |
| 1548 |
> |
$\tau_2^y$, placing the result more in line with that obtained using |
| 1549 |
> |
the Ewald sum. This example supports this explanation; however, |
| 1550 |
> |
recomputing the results to meet a poorer statistical standard is |
| 1551 |
> |
counter-productive. Assuming the Ewald results are not the product of |
| 1552 |
> |
poor statistics, differences in techniques to integrate the |
| 1553 |
> |
orientational motion could also play a role. {\sc shake} is the most |
| 1554 |
> |
commonly used technique for approximating rigid-body orientational |
| 1555 |
> |
motion,\cite{Ryckaert77} whereas in {\sc oopse}, we maintain and |
| 1556 |
> |
integrate the entire rotation matrix using the {\sc dlm} |
| 1557 |
> |
method.\cite{Meineke05} Since {\sc shake} is an iterative constraint |
| 1558 |
> |
technique, if the convergence tolerances are raised for increased |
| 1559 |
> |
performance, error will accumulate in the orientational |
| 1560 |
> |
motion. Finally, the Ewald results were calculated using the $NVT$ |
| 1561 |
> |
ensemble, while the $NVE$ ensemble was used for {\sc sf} |
| 1562 |
|
calculations. The additional mode of motion due to the thermostat will |
| 1563 |
|
alter the dynamics, resulting in differences between $NVT$ and $NVE$ |
| 1564 |
|
results. These differences are increasingly noticeable as the |
| 1570 |
|
neutralizing the cutoff sphere with charge-charge interaction shifting |
| 1571 |
|
and by damping the electrostatic interactions. Now we would like to |
| 1572 |
|
consider an extension of these techniques to include point multipole |
| 1573 |
< |
interactions. How will the shifting and damping need to develop in |
| 1573 |
> |
interactions. How will the shifting and damping need to be modified in |
| 1574 |
|
order to accommodate point multipoles? |
| 1575 |
|
|
| 1576 |
< |
Of the two techniques, the least to vary is shifting. Shifting is |
| 1576 |
> |
Of the two techniques, the easiest to adapt is shifting. Shifting is |
| 1577 |
|
employed to neutralize the cutoff sphere; however, in a system |
| 1578 |
|
composed purely of point multipoles, the cutoff sphere is already |
| 1579 |
|
neutralized. This means that shifting is not necessary between point |
| 1580 |
|
multipoles. In a mixed system of monopoles and multipoles, the |
| 1581 |
|
undamped {\sc sf} potential needs only to shift the force terms of the |
| 1582 |
|
monopole (and use the monopole potential of equation (\ref{eq:SFPot})) |
| 1583 |
< |
and smoothly cutoff the multipole interactions with a switching |
| 1583 |
> |
and smoothly truncate the multipole interactions with a switching |
| 1584 |
|
function. The switching function is required in order to conserve |
| 1585 |
< |
energy, because a discontinuity will exist at $R_\textrm{c}$ in the |
| 1586 |
< |
absence of shifting terms. |
| 1585 |
> |
energy, because a discontinuity will exist in both the potential and |
| 1586 |
> |
forces at $R_\textrm{c}$ in the absence of shifting terms. |
| 1587 |
|
|
| 1588 |
|
If we consider damping the {\sc sf} potential (Eq. (\ref{eq:DSFPot})), |
| 1589 |
|
then we need to incorporate the complimentary error function term into |
| 1591 |
|
replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole |
| 1592 |
|
expansion.\cite{Hirschfelder67} In the multipole expansion, rather |
| 1593 |
|
than considering only the interactions between single point charges, |
| 1594 |
< |
the electrostatic interactions is reformulated such that it describes |
| 1594 |
> |
the electrostatic interaction is reformulated such that it describes |
| 1595 |
|
the interaction between charge distributions about central sites of |
| 1596 |
|
the respective sets of charges. This procedure is what leads to the |
| 1597 |
|
familiar charge-dipole, |
| 1674 |
|
\end{equation} |
| 1675 |
|
Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional |
| 1676 |
|
term. Continuing with higher rank tensors, we can obtain the damping |
| 1677 |
< |
functions for higher multipoles as well as the forces. Each subsequent |
| 1677 |
> |
functions for higher multipole potentials and forces. Each subsequent |
| 1678 |
|
damping function includes one additional term, and we can simplify the |
| 1679 |
|
procedure for obtaining these terms by writing out the following |
| 1680 |
|
generating function, |
| 1708 |
|
c_1(r_{ij}), |
| 1709 |
|
\label{eq:dampDipoleDipole} |
| 1710 |
|
\end{equation} |
| 1711 |
< |
$c_2(r_{ij})$ and $c_1(r_{ij})$ respectively dampen these two |
| 1712 |
< |
parts. The forces for the damped dipole-dipole interaction, |
| 1711 |
> |
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts |
| 1712 |
> |
respectively. The forces for the damped dipole-dipole interaction, |
| 1713 |
|
\begin{equation} |
| 1714 |
|
\begin{split} |
| 1715 |
|
F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
| 1739 |
|
constant is calculated from the long-time fluctuations of the system's |
| 1740 |
|
accumulated dipole moment (Eq. (\ref{eq:staticDielectric})), so it is |
| 1741 |
|
going to be quite sensitive to the choice of damping parameter. We |
| 1742 |
< |
would like to choose an optimal damping constant for any particular |
| 1743 |
< |
cutoff radius choice that would properly capture the dielectric |
| 1744 |
< |
behavior of the liquid. |
| 1742 |
> |
would like to choose optimal damping constants such that any arbitrary |
| 1743 |
> |
choice of cutoff radius will properly capture the dielectric behavior |
| 1744 |
> |
of the liquid. |
| 1745 |
|
|
| 1746 |
|
In order to find these optimal values, we mapped out the static |
| 1747 |
|
dielectric constant as a function of both the damping parameter and |
| 1752 |
|
four-point transferable intermolecular potential (TIP4P) for water |
| 1753 |
|
targeted for use with the Ewald summation.\cite{Horn04} SSD/RF is the |
| 1754 |
|
reaction field modified variant of the soft sticky dipole (SSD) model |
| 1755 |
< |
for water\cite{Fennell04} This model is discussed in more detail in |
| 1755 |
> |
for water.\cite{Fennell04} This model is discussed in more detail in |
| 1756 |
|
the next chapter. One thing to note about it, electrostatic |
| 1757 |
|
interactions are handled via dipole-dipole interactions rather than |
| 1758 |
|
charge-charge interactions like the other three models. Damping of the |
| 1760 |
|
\ref{sec:dampingMultipoles}. Each of these systems were studied with |
| 1761 |
|
cutoff radii of 9, 10, 11, and 12~\AA\ and with damping parameter values |
| 1762 |
|
ranging from 0 to 0.35~\AA$^{-1}$. |
| 1763 |
+ |
|
| 1764 |
|
\begin{figure} |
| 1765 |
|
\centering |
| 1766 |
|
\includegraphics[width=\linewidth]{./figures/dielectricMap.pdf} |
| 1769 |
|
radius ($R_\textrm{c}$) and damping coefficient ($\alpha$).} |
| 1770 |
|
\label{fig:dielectricMap} |
| 1771 |
|
\end{figure} |
| 2611 |
– |
|
| 1772 |
|
The results of these calculations are displayed in figure |
| 1773 |
|
\ref{fig:dielectricMap} in the form of shaded contour plots. An |
| 1774 |
|
interesting aspect of all four contour plots is that the dielectric |
| 1785 |
|
with {\sc sf} these parameters give a dielectric constant of |
| 1786 |
|
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where |
| 1787 |
|
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and |
| 1788 |
< |
0.35~\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal |
| 1789 |
< |
to 63$\pm$1.\cite{Horn04} We did not perform calculations with these |
| 1790 |
< |
exact parameters, but interpolating between surrounding values gives a |
| 1791 |
< |
$\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric |
| 1792 |
< |
constant on $\alpha$ and $R_\textrm{c}$ with the {\sc sf} technique, |
| 1793 |
< |
it might be interesting to investigate the dielectric dependence of |
| 1794 |
< |
the real-space Ewald parameters. |
| 1788 |
> |
0.35~\AA$^{-1}$, and these parameters resulted in a dielectric |
| 1789 |
> |
constant equal to 63$\pm$1.\cite{Horn04} We did not perform |
| 1790 |
> |
calculations with these exact parameters, but interpolating between |
| 1791 |
> |
surrounding values gives a dielectric constant of 61$\pm$1. Since the |
| 1792 |
> |
dielectric constant is dependent on $\alpha$ and $R_\textrm{c}$ with |
| 1793 |
> |
the {\sc sf} technique, it might be interesting to investigate the |
| 1794 |
> |
dielectric dependence of the real-space Ewald parameters. |
| 1795 |
|
|
| 1796 |
|
Although it is tempting to choose damping parameters equivalent to |
| 1797 |
|
these Ewald examples, the results discussed in sections |
| 1798 |
< |
\ref{sec:EnergyResults} through \ref{sec:IndividualResults} indicate |
| 1799 |
< |
that values this high are destructive to both the energetics and |
| 1800 |
< |
dynamics. Ideally, $\alpha$ should not exceed 0.3~\AA$^{-1}$ for any of |
| 1801 |
< |
the cutoff values in this range. If the optimal damping parameter is |
| 1802 |
< |
chosen to be midway between 0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) |
| 1803 |
< |
at the 9~\AA\ cutoff, then the linear trend with $R_\textrm{c}$ will |
| 1804 |
< |
always keep $\alpha$ below 0.3~\AA$^{-1}$. This linear progression |
| 1805 |
< |
would give values of 0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for |
| 1806 |
< |
cutoff radii of 9, 10, 11, and 12~\AA. Setting this to be the default |
| 1807 |
< |
behavior for the damped {\sc sf} technique will result in consistent |
| 1808 |
< |
dielectric behavior for these and other condensed molecular systems, |
| 1809 |
< |
regardless of the chosen cutoff radius. The static dielectric |
| 1810 |
< |
constants for TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be |
| 1811 |
< |
approximately fixed at 74, 52, 58, and 89 respectively. These values |
| 1812 |
< |
are generally lower than the values reported in the literature; |
| 1813 |
< |
however, the relative dielectric behavior scales as expected when |
| 1814 |
< |
comparing the models to one another. |
| 1798 |
> |
\ref{sec:EnergyResults} through \ref{sec:FTDirResults} and appendix |
| 1799 |
> |
\ref{app:IndividualResults} indicate that values this high are |
| 1800 |
> |
destructive to both the energetics and dynamics. Ideally, $\alpha$ |
| 1801 |
> |
should not exceed 0.3~\AA$^{-1}$ for any of the cutoff values in this |
| 1802 |
> |
range. If the optimal damping parameter is chosen to be midway between |
| 1803 |
> |
0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) at the 9~\AA\ cutoff, |
| 1804 |
> |
then the linear trend with $R_\textrm{c}$ will always keep $\alpha$ |
| 1805 |
> |
below 0.3~\AA$^{-1}$. This linear progression would give values of |
| 1806 |
> |
0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for cutoff radii of 9, |
| 1807 |
> |
10, 11, and 12~\AA. Setting this to be the default behavior for the |
| 1808 |
> |
damped {\sc sf} technique will result in consistent dielectric |
| 1809 |
> |
behavior for these and other condensed molecular systems, regardless |
| 1810 |
> |
of the chosen cutoff radius. The static dielectric constants for |
| 1811 |
> |
TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be fixed at approximately |
| 1812 |
> |
74, 52, 58, and 89 respectively. These values are generally lower than |
| 1813 |
> |
the values reported in the literature; however, the relative |
| 1814 |
> |
dielectric behavior scales as expected when comparing the models to |
| 1815 |
> |
one another. |
| 1816 |
|
|
| 1817 |
|
\section{Conclusions}\label{sec:PairwiseConclusions} |
| 1818 |
|
|
| 1825 |
|
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1826 |
|
energetic and dynamic characteristics exhibited by simulations |
| 1827 |
|
employing lattice summation techniques. The cumulative energy |
| 1828 |
< |
difference results showed the undamped {\sc sf} and moderately damped |
| 1829 |
< |
{\sc sp} methods produced results nearly identical to the Ewald |
| 1828 |
> |
difference results showed that the undamped {\sc sf} and moderately |
| 1829 |
> |
damped {\sc sp} methods produce results nearly identical to the Ewald |
| 1830 |
|
summation. Similarly for the dynamic features, the undamped or |
| 1831 |
|
moderately damped {\sc sf} and moderately damped {\sc sp} methods |
| 1832 |
|
produce force and torque vector magnitude and directions very similar |
| 1839 |
|
As in all purely-pairwise cutoff methods, these methods are expected |
| 1840 |
|
to scale approximately {\it linearly} with system size, and they are |
| 1841 |
|
easily parallelizable. This should result in substantial reductions |
| 1842 |
< |
in the computational cost of performing large simulations. |
| 1842 |
> |
in the computational cost associated with large-scale simulations. |
| 1843 |
|
|
| 1844 |
|
Aside from the computational cost benefit, these techniques have |
| 1845 |
|
applicability in situations where the use of the Ewald sum can prove |
| 1858 |
|
systems containing point charges, most structural features will be |
| 1859 |
|
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1860 |
|
method with an electrostatic damping of 0.2~\AA$^{-1}$. These methods |
| 1861 |
< |
would also be appropriate for molecular dynamics simulations where the |
| 1862 |
< |
data of interest is either structural or short-time dynamical |
| 1861 |
> |
would also be appropriate in molecular dynamics simulations where the |
| 1862 |
> |
data of interest are either structural or short-time dynamical |
| 1863 |
|
quantities. For long-time dynamics and collective motions, the safest |
| 1864 |
|
pairwise method we have evaluated is the {\sc sf} method with an |
| 1865 |
|
electrostatic damping between 0.2 and 0.25~\AA$^{-1}$. It is also |
| 1868 |
|
$R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc |
| 1869 |
|
sf} method should use an $\alpha$ of 0.2175~\AA$^{-1}$ for an |
| 1870 |
|
$R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease by |
| 1871 |
< |
0.025~\AA$^{-1}$ for every 1~\AA\ increase in cutoff radius. |
| 1871 |
> |
0.025~\AA$^{-1}$ for every 1~\AA\ increase in the cutoff radius. |
| 1872 |
|
|
| 1873 |
|
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1874 |
< |
fact, it is the standard by which these simple pairwise sums have been |
| 1875 |
< |
judged. However, these results do suggest that in the typical |
| 1874 |
> |
fact, it is the standard by which these simple pairwise methods have |
| 1875 |
> |
been judged. However, these results do suggest that in the typical |
| 1876 |
|
simulations performed today, the Ewald summation may no longer be |
| 1877 |
|
required to obtain the level of accuracy most researchers have come to |
| 1878 |
|
expect. |
