| 1 |
|
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES} |
| 2 |
|
|
| 3 |
|
In molecular simulations, proper accumulation of electrostatic |
| 4 |
< |
interactions is essential and is one of the most |
| 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 every |
| 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$. |
| 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} |
| 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. |
| 168 |
|
|
| 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 |
| 338 |
< |
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 |
|
|
| 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 |
| 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 |
|
|
| 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} |
| 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 |
|
|
| 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 |
|
|
| 514 |
– |
Results and discussion for the individual analysis of each of the |
| 515 |
– |
system types appear in appendix \ref{app:IndividualResults}, while the |
| 516 |
– |
cumulative results over all the investigated systems appear below in |
| 517 |
– |
sections \ref{sec:EnergyResults}. |
| 518 |
– |
|
| 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} |
| 612 |
– |
|
| 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 |
| 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 |
< |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 758 |
< |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 759 |
< |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 760 |
< |
This is perhaps expected because this method explicitly incorporates a |
| 761 |
< |
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 |
| 794 |
|
\end{figure} |
| 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 reproduce 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 |
| 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 (Fig. \ref{fig:frcTrqAng}) are compared |
| 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 |
|
|
| 835 |
|
\end{figure} |
| 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) |
| 845 |
< |
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 |
| 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. |
| 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 |
|
|
| 937 |
|
|
| 938 |
|
\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics} |
| 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 |
< |
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}$, |
| 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 |
< |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 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. |
| 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 |
| 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. |
| 1039 |
< |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1040 |
< |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1041 |
< |
However, even without any electrostatic damping, the {\sc sf} method |
| 1042 |
< |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1043 |
< |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1044 |
< |
would predict the lowest frequency peak near 325~cm$^{-1}$. {\it |
| 1045 |
< |
Most} of the collective behavior in the crystal is accurately captured |
| 1046 |
< |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1047 |
< |
obtained using moderate damping in addition to the shifting at the |
| 1048 |
< |
cutoff distance. |
| 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 |
|
|
