| 1 |
|
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
| 2 |
< |
\documentclass[12pt]{article} |
| 2 |
> |
%\documentclass[aps,prb,preprint]{revtex4} |
| 3 |
> |
\documentclass[11pt]{article} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
< |
\usepackage{amsmath} |
| 5 |
> |
\usepackage{amsmath,bm} |
| 6 |
|
\usepackage{amssymb} |
| 6 |
– |
%\usepackage{ifsym} |
| 7 |
|
\usepackage{epsf} |
| 8 |
|
\usepackage{times} |
| 9 |
< |
\usepackage{mathptm} |
| 9 |
> |
\usepackage{mathptmx} |
| 10 |
|
\usepackage{setspace} |
| 11 |
|
\usepackage{tabularx} |
| 12 |
|
\usepackage{graphicx} |
| 13 |
|
\usepackage{booktabs} |
| 14 |
|
\usepackage{bibentry} |
| 15 |
|
\usepackage{mathrsfs} |
| 16 |
– |
%\usepackage{berkeley} |
| 16 |
|
\usepackage[ref]{overcite} |
| 17 |
|
\pagestyle{plain} |
| 18 |
|
\pagenumbering{arabic} |
| 20 |
|
\topmargin -21pt \headsep 10pt |
| 21 |
|
\textheight 9.0in \textwidth 6.5in |
| 22 |
|
\brokenpenalty=10000 |
| 23 |
< |
\renewcommand{\baselinestretch}{1.2} |
| 23 |
> |
%\renewcommand{\baselinestretch}{1.2} |
| 24 |
> |
\renewcommand{\baselinestretch}{2} |
| 25 |
|
\renewcommand\citemid{\ } % no comma in optional reference note |
| 26 |
+ |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions |
| 27 |
+ |
\let\Caption\caption |
| 28 |
+ |
\renewcommand\caption[1]{% |
| 29 |
+ |
\Caption[#1]{}% |
| 30 |
+ |
} |
| 31 |
|
|
| 32 |
+ |
|
| 33 |
|
\begin{document} |
| 34 |
|
|
| 35 |
< |
\title{Is the Ewald Summation necessary in typical molecular simulations: Alternatives to the accepted standard of cutoff policies} |
| 35 |
> |
\title{Is the Ewald summation still necessary? \\ |
| 36 |
> |
Pairwise alternatives to the accepted standard for |
| 37 |
> |
long-range electrostatics in molecular simulations} |
| 38 |
|
|
| 39 |
< |
\author{Christopher J. Fennell and J. Daniel Gezelter \\ |
| 39 |
> |
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 40 |
> |
gezelter@nd.edu} \\ |
| 41 |
|
Department of Chemistry and Biochemistry\\ |
| 42 |
|
University of Notre Dame\\ |
| 43 |
|
Notre Dame, Indiana 46556} |
| 46 |
|
|
| 47 |
|
\maketitle |
| 48 |
|
%\doublespacing |
| 49 |
< |
\nobibliography{} |
| 49 |
> |
|
| 50 |
|
\begin{abstract} |
| 51 |
< |
A new method for accumulating electrostatic interactions was derived from the previous efforts described in \bibentry{Wolf99} and \bibentry{Zahn02} as a possible replacement for lattice sum methods in molecular simulations. Comparisons were performed with this and other pairwise electrostatic summation techniques against the smooth particle mesh Ewald (SPME) summation to see how well they reproduce the energetics and dynamics of a variety of simulation types. The newly derived Shifted-Force technique shows a remarkable ability to reproduce the behavior exhibited in simulations using SPME with an $\mathscr{O}(N)$ computational cost, equivalent to merely the real-space portion of the lattice summation. |
| 51 |
> |
We investigate pairwise electrostatic interaction methods and show |
| 52 |
> |
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
| 53 |
> |
alternatives to the Ewald summation for typical modern molecular |
| 54 |
> |
simulations. These methods are extended from the damped and |
| 55 |
> |
cutoff-neutralized Coulombic sum originally proposed by |
| 56 |
> |
[D. Wolf, P. Keblinski, S.~R. Phillpot, and J. Eggebrecht, {\it J. Chem. Phys.} {\bf 110}, 8255 (1999)] One of these, the damped shifted force method, shows |
| 57 |
> |
a remarkable ability to reproduce the energetic and dynamic |
| 58 |
> |
characteristics exhibited by simulations employing lattice summation |
| 59 |
> |
techniques. Comparisons were performed with this and other pairwise |
| 60 |
> |
methods against the smooth particle mesh Ewald ({\sc spme}) summation |
| 61 |
> |
to see how well they reproduce the energetics and dynamics of a |
| 62 |
> |
variety of simulation types. |
| 63 |
|
\end{abstract} |
| 64 |
|
|
| 65 |
+ |
\newpage |
| 66 |
+ |
|
| 67 |
|
%\narrowtext |
| 68 |
|
|
| 69 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 69 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 70 |
|
% BODY OF TEXT |
| 71 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 71 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 72 |
|
|
| 73 |
|
\section{Introduction} |
| 74 |
|
|
| 75 |
< |
In molecular simulations, proper accumulation of the electrostatic interactions is considered one of the most essential and computationally demanding tasks. |
| 75 |
> |
In molecular simulations, proper accumulation of the electrostatic |
| 76 |
> |
interactions is essential and is one of the most |
| 77 |
> |
computationally-demanding tasks. The common molecular mechanics force |
| 78 |
> |
fields represent atomic sites with full or partial charges protected |
| 79 |
> |
by Lennard-Jones (short range) interactions. This means that nearly |
| 80 |
> |
every pair interaction involves a calculation of charge-charge forces. |
| 81 |
> |
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
| 82 |
> |
interactions quickly become the most expensive part of molecular |
| 83 |
> |
simulations. Historically, the electrostatic pair interaction would |
| 84 |
> |
not have decayed appreciably within the typical box lengths that could |
| 85 |
> |
be feasibly simulated. In the larger systems that are more typical of |
| 86 |
> |
modern simulations, large cutoffs should be used to incorporate |
| 87 |
> |
electrostatics correctly. |
| 88 |
|
|
| 89 |
< |
blah blah blah Ewald Sum Important blah blah blah |
| 89 |
> |
There have been many efforts to address the proper and practical |
| 90 |
> |
handling of electrostatic interactions, and these have resulted in a |
| 91 |
> |
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are |
| 92 |
> |
typically classified as implicit methods (i.e., continuum dielectrics, |
| 93 |
> |
static dipolar fields),\cite{Born20,Grossfield00} explicit methods |
| 94 |
> |
(i.e., Ewald summations, interaction shifting or |
| 95 |
> |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 96 |
> |
reaction field type methods, fast multipole |
| 97 |
> |
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 98 |
> |
often preferred because they physically incorporate solvent molecules |
| 99 |
> |
in the system of interest, but these methods are sometimes difficult |
| 100 |
> |
to utilize because of their high computational cost.\cite{Roux99} In |
| 101 |
> |
addition to the computational cost, there have been some questions |
| 102 |
> |
regarding possible artifacts caused by the inherent periodicity of the |
| 103 |
> |
explicit Ewald summation.\cite{Tobias01} |
| 104 |
|
|
| 105 |
< |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined for accumulation of electrostatic interactions in a simple pairwise fashion.\cite{Wolf99} They took the observation that the electrostatic interaction is short-ranged in systems of charges and that charge neutralization within the cutoff spheres is crucial for potential stability, and they devised a pairwise summation method that ensures charge neutrality and gives results similar to those obtained using the Ewald summation. The resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through placement on the cutoff sphere and a distance-dependent damping function (identical to that seen in the real-space portion of the Ewald sum) to hasten energetic convergence |
| 105 |
> |
In this paper, we focus on a new set of pairwise methods devised by |
| 106 |
> |
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 107 |
> |
methods along with a few other mixed methods (i.e. reaction field) are |
| 108 |
> |
compared with the smooth particle mesh Ewald |
| 109 |
> |
sum,\cite{Onsager36,Essmann99} which is our reference method for |
| 110 |
> |
handling long-range electrostatic interactions. The new methods for |
| 111 |
> |
handling electrostatics have the potential to scale linearly with |
| 112 |
> |
increasing system size since they involve only a simple modification |
| 113 |
> |
to the direct pairwise sum. They also lack the added periodicity of |
| 114 |
> |
the Ewald sum, so they can be used for systems which are non-periodic |
| 115 |
> |
or which have one- or two-dimensional periodicity. Below, these |
| 116 |
> |
methods are evaluated using a variety of model systems to |
| 117 |
> |
establish their usability in molecular simulations. |
| 118 |
> |
|
| 119 |
> |
\subsection{The Ewald Sum} |
| 120 |
> |
The complete accumulation of the electrostatic interactions in a system with |
| 121 |
> |
periodic boundary conditions (PBC) requires the consideration of the |
| 122 |
> |
effect of all charges within a (cubic) simulation box as well as those |
| 123 |
> |
in the periodic replicas, |
| 124 |
|
\begin{equation} |
| 125 |
< |
V(r_{ij})= \frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 125 |
> |
V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime \left[ \sum_{i=1}^N\sum_{j=1}^N \phi\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) \right], |
| 126 |
> |
\label{eq:PBCSum} |
| 127 |
> |
\end{equation} |
| 128 |
> |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
| 129 |
> |
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
| 130 |
> |
prime indicates $i = j$ are neglected for $\mathbf{n} = |
| 131 |
> |
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
| 132 |
> |
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
| 133 |
> |
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
| 134 |
> |
$j$, and $\phi$ is the solution to Poisson's equation |
| 135 |
> |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for |
| 136 |
> |
charge-charge interactions). In the case of monopole electrostatics, |
| 137 |
> |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for |
| 138 |
> |
non-neutral systems. |
| 139 |
> |
|
| 140 |
> |
The electrostatic summation problem was originally studied by Ewald |
| 141 |
> |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 142 |
> |
took was to convert this conditionally convergent sum into two |
| 143 |
> |
absolutely convergent summations: a short-ranged real-space summation |
| 144 |
> |
and a long-ranged reciprocal-space summation, |
| 145 |
> |
\begin{equation} |
| 146 |
> |
\begin{split} |
| 147 |
> |
V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
| 148 |
> |
\end{split} |
| 149 |
> |
\label{eq:EwaldSum} |
| 150 |
> |
\end{equation} |
| 151 |
> |
where $\alpha$ is the damping or convergence parameter with units of |
| 152 |
> |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 153 |
> |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 154 |
> |
constant of the surrounding medium. The final two terms of |
| 155 |
> |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 156 |
> |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 157 |
> |
dipolar term was neglected in early applications in molecular |
| 158 |
> |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 159 |
> |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 160 |
> |
dipole moment which is magnified through replication of the periodic |
| 161 |
> |
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
| 162 |
> |
system is said to be using conducting (or ``tin-foil'') boundary |
| 163 |
> |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 164 |
> |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 165 |
> |
time. Initially, due to the small system sizes that could be |
| 166 |
> |
simulated feasibly, the entire simulation box was replicated to |
| 167 |
> |
convergence. In more modern simulations, the systems have grown large |
| 168 |
> |
enough that a real-space cutoff could potentially give convergent |
| 169 |
> |
behavior. Indeed, it has been observed that with the choice of a |
| 170 |
> |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 171 |
> |
rapidly convergent and small relative to the real-space |
| 172 |
> |
portion.\cite{Karasawa89,Kolafa92} |
| 173 |
> |
|
| 174 |
> |
\begin{figure} |
| 175 |
> |
\centering |
| 176 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 177 |
> |
\caption{The change in the need for the Ewald sum with |
| 178 |
> |
increasing computational power. A:~Initially, only small systems |
| 179 |
> |
could be studied, and the Ewald sum replicated the simulation box to |
| 180 |
> |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 181 |
> |
convergence for the larger systems of charges that are common today.} |
| 182 |
> |
\label{fig:ewaldTime} |
| 183 |
> |
\end{figure} |
| 184 |
> |
|
| 185 |
> |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 186 |
> |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 187 |
> |
the computational cost between the direct and reciprocal-space |
| 188 |
> |
portions of the summation. The choice of this value allows one to |
| 189 |
> |
select whether the real-space or reciprocal space portion of the |
| 190 |
> |
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
| 191 |
> |
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 192 |
> |
$\alpha$ and thoughtful algorithm development, this cost can be |
| 193 |
> |
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 194 |
> |
taken to reduce the cost of the Ewald summation even further is to set |
| 195 |
> |
$\alpha$ such that the real-space interactions decay rapidly, allowing |
| 196 |
> |
for a short spherical cutoff. Then the reciprocal space summation is |
| 197 |
> |
optimized. These optimizations usually involve utilization of the |
| 198 |
> |
fast Fourier transform (FFT),\cite{Hockney81} leading to the |
| 199 |
> |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 200 |
> |
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 201 |
> |
methods, the cost of the reciprocal-space portion of the Ewald |
| 202 |
> |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
| 203 |
> |
\log N)$. |
| 204 |
> |
|
| 205 |
> |
These developments and optimizations have made the use of the Ewald |
| 206 |
> |
summation routine in simulations with periodic boundary |
| 207 |
> |
conditions. However, in certain systems, such as vapor-liquid |
| 208 |
> |
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 209 |
> |
can prove problematic. The Ewald sum has been reformulated to handle |
| 210 |
> |
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the |
| 211 |
> |
new methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 212 |
> |
recently, there have been several successful efforts toward reducing |
| 213 |
> |
the computational cost of 2D lattice summations, often enabling the |
| 214 |
> |
use of the mentioned |
| 215 |
> |
optimizations.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 216 |
> |
|
| 217 |
> |
Several studies have recognized that the inherent periodicity in the |
| 218 |
> |
Ewald sum can also have an effect on three-dimensional |
| 219 |
> |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 220 |
> |
Solvated proteins are essentially kept at high concentration due to |
| 221 |
> |
the periodicity of the electrostatic summation method. In these |
| 222 |
> |
systems, the more compact folded states of a protein can be |
| 223 |
> |
artificially stabilized by the periodic replicas introduced by the |
| 224 |
> |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 225 |
> |
considering the use of the Ewald summation where the assumed |
| 226 |
> |
periodicity would introduce spurious effects in the system dynamics. |
| 227 |
> |
|
| 228 |
> |
\subsection{The Wolf and Zahn Methods} |
| 229 |
> |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 230 |
> |
for the accurate accumulation of electrostatic interactions in an |
| 231 |
> |
efficient pairwise fashion. This procedure lacks the inherent |
| 232 |
> |
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
| 233 |
> |
observed that the electrostatic interaction is effectively |
| 234 |
> |
short-ranged in condensed phase systems and that neutralization of the |
| 235 |
> |
charge contained within the cutoff radius is crucial for potential |
| 236 |
> |
stability. They devised a pairwise summation method that ensures |
| 237 |
> |
charge neutrality and gives results similar to those obtained with the |
| 238 |
> |
Ewald summation. The resulting shifted Coulomb potential includes |
| 239 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 240 |
> |
and a distance-dependent damping function (identical to that seen in |
| 241 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 242 |
> |
\begin{equation} |
| 243 |
> |
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 244 |
|
\label{eq:WolfPot} |
| 245 |
|
\end{equation} |
| 246 |
< |
In order to use this potential in molecular dynamics simulations, the derivative of this potential was taken, followed by evaluation of the limit to give the following forces, |
| 246 |
> |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 247 |
> |
potential. However, neutralizing the charge contained within each |
| 248 |
> |
cutoff sphere requires the placement of a self-image charge on the |
| 249 |
> |
surface of the cutoff sphere. This additional self-term in the total |
| 250 |
> |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
| 251 |
> |
Madelung energies for many crystals. |
| 252 |
> |
|
| 253 |
> |
In order to use their charge-neutralized potential in molecular |
| 254 |
> |
dynamics simulations, Wolf \textit{et al.} suggested taking the |
| 255 |
> |
derivative of this potential prior to evaluation of the limit. This |
| 256 |
> |
procedure gives an expression for the forces, |
| 257 |
|
\begin{equation} |
| 258 |
< |
F(r_{ij}) = q_iq_j\left[\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right)-\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right]. |
| 258 |
> |
F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{\left[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}\right]-\left[\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right]\right\}, |
| 259 |
|
\label{eq:WolfForces} |
| 260 |
|
\end{equation} |
| 261 |
< |
More recently, Zahn \textit{et al.} investigated this electrostatic summation method for use in simulations involving water.\cite{Zahn02} In their work, they point out that the method as proposed is problematic for use in Molecular Dynamics simulations, because the integral of the forces and derivative of the potential are not equivalent. This comes about from the procedure of taking the limit shown in equation \ref{eq:WolfPot} after calculating the derivatives.\cite{Wolf99} Zahn \textit{et al.} proposed a shifted force adaptation of this "Wolf summation method" as a way to use this technique in Molecular Dynamics simulations. Taking the integral of the forces shown in equation \ref{eq:WolfForces}, they obtained a new shifted damped Coulomb potential |
| 261 |
> |
that incorporates both image charges and damping of the electrostatic |
| 262 |
> |
interaction. |
| 263 |
> |
|
| 264 |
> |
More recently, Zahn \textit{et al.} investigated these potential and |
| 265 |
> |
force expressions for use in simulations involving water.\cite{Zahn02} |
| 266 |
> |
In their work, they pointed out that the forces and derivative of |
| 267 |
> |
the potential are not commensurate. Attempts to use both |
| 268 |
> |
eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 269 |
> |
to poor energy conservation. They correctly observed that taking the |
| 270 |
> |
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 271 |
> |
derivatives gives forces for a different potential energy function |
| 272 |
> |
than the one shown in eq. (\ref{eq:WolfPot}). |
| 273 |
> |
|
| 274 |
> |
Zahn \textit{et al.} introduced a modified form of this summation |
| 275 |
> |
method as a way to use the technique in Molecular Dynamics |
| 276 |
> |
simulations. They proposed a new damped Coulomb potential, |
| 277 |
|
\begin{equation} |
| 278 |
< |
V\left(r_{ij}\right)=q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. |
| 278 |
> |
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}, |
| 279 |
|
\label{eq:ZahnPot} |
| 280 |
|
\end{equation} |
| 281 |
< |
They showed that this new potential does well in capturing the structural and dynamic properties of water in their simulations. |
| 281 |
> |
and showed that this potential does fairly well at capturing the |
| 282 |
> |
structural and dynamic properties of water compared the same |
| 283 |
> |
properties obtained using the Ewald sum. |
| 284 |
|
|
| 285 |
< |
While implementing these methods for use in our own work, we discovered the potential presented in equation \ref{eq:ZahnPot} is still not entirely correct. The derivative of this equation leads to a sign error in the forces, resulting in erroneous dynamics. We can apply the standard shifted force potential form for Lennard-Jones potentials, |
| 285 |
> |
\subsection{Simple Forms for Pairwise Electrostatics} |
| 286 |
> |
|
| 287 |
> |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 288 |
> |
al.} are constructed using two different (and separable) computational |
| 289 |
> |
tricks: \begin{enumerate} |
| 290 |
> |
\item shifting through the use of image charges, and |
| 291 |
> |
\item damping the electrostatic interaction. |
| 292 |
> |
\end{enumerate} Wolf \textit{et al.} treated the |
| 293 |
> |
development of their summation method as a progressive application of |
| 294 |
> |
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 295 |
> |
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 296 |
> |
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
| 297 |
> |
both techniques. It is possible, however, to separate these |
| 298 |
> |
tricks and study their effects independently. |
| 299 |
> |
|
| 300 |
> |
Starting with the original observation that the effective range of the |
| 301 |
> |
electrostatic interaction in condensed phases is considerably less |
| 302 |
> |
than $r^{-1}$, either the cutoff sphere neutralization or the |
| 303 |
> |
distance-dependent damping technique could be used as a foundation for |
| 304 |
> |
a new pairwise summation method. Wolf \textit{et al.} made the |
| 305 |
> |
observation that charge neutralization within the cutoff sphere plays |
| 306 |
> |
a significant role in energy convergence; therefore we will begin our |
| 307 |
> |
analysis with the various shifted forms that maintain this charge |
| 308 |
> |
neutralization. We can evaluate the methods of Wolf |
| 309 |
> |
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 310 |
> |
shifted potential, |
| 311 |
|
\begin{equation} |
| 312 |
< |
V^\textrm{SF}(r_{ij}) = \begin{cases} v(r_{ij})-v_\textrm{c}-\left(\frac{\textrm{d}v(r_{ij})}{\textrm{d}r_{ij}}\right)_{r_{ij}=R_\textrm{c}}(r_{ij}-R_\textrm{c}) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 313 |
< |
\end{cases}, |
| 312 |
> |
V_\textrm{SP}(r) = \begin{cases} |
| 313 |
> |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 314 |
> |
R_\textrm{c} |
| 315 |
> |
\end{cases}, |
| 316 |
> |
\label{eq:shiftingPotForm} |
| 317 |
|
\end{equation} |
| 318 |
< |
where $v(r_{ij})$ unshifted form of the potential, and $v_c$ is a constant term that insures the potential goes to zero at the cutoff radius.\cite{Allen87} Using the simple damped Coulomb potential as the starting point, |
| 318 |
> |
and shifted force, |
| 319 |
|
\begin{equation} |
| 320 |
< |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 321 |
< |
\label{eq:dampCoulomb} |
| 322 |
< |
\end{equation} |
| 323 |
< |
the resulting shifted force potential is |
| 320 |
> |
V_\textrm{SF}(r) = \begin{cases} |
| 321 |
> |
v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 322 |
> |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 323 |
> |
\end{cases}, |
| 324 |
> |
\label{eq:shiftingForm} |
| 325 |
> |
\end{equation} |
| 326 |
> |
functions where $v(r)$ is the unshifted form of the potential, and |
| 327 |
> |
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
| 328 |
> |
that both the potential and the forces goes to zero at the cutoff |
| 329 |
> |
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
| 330 |
> |
potential is smooth at the cutoff radius |
| 331 |
> |
($R_\textrm{c}$).\cite{Allen87} |
| 332 |
> |
|
| 333 |
> |
The forces associated with the shifted potential are simply the forces |
| 334 |
> |
of the unshifted potential itself (when inside the cutoff sphere), |
| 335 |
|
\begin{equation} |
| 336 |
< |
V^\mathrm{SF}\left(r_{ij}\right)=q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}}+\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. |
| 336 |
> |
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 337 |
> |
\end{equation} |
| 338 |
> |
and are zero outside. Inside the cutoff sphere, the forces associated |
| 339 |
> |
with the shifted force form can be written, |
| 340 |
> |
\begin{equation} |
| 341 |
> |
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 342 |
> |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 343 |
> |
\end{equation} |
| 344 |
> |
|
| 345 |
> |
If the potential, $v(r)$, is taken to be the normal Coulomb potential, |
| 346 |
> |
\begin{equation} |
| 347 |
> |
v(r) = \frac{q_i q_j}{r}, |
| 348 |
> |
\label{eq:Coulomb} |
| 349 |
> |
\end{equation} |
| 350 |
> |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 351 |
> |
al.}'s undamped prescription: |
| 352 |
> |
\begin{equation} |
| 353 |
> |
V_\textrm{SP}(r) = |
| 354 |
> |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 355 |
> |
r\leqslant R_\textrm{c}, |
| 356 |
> |
\label{eq:SPPot} |
| 357 |
> |
\end{equation} |
| 358 |
> |
with associated forces, |
| 359 |
> |
\begin{equation} |
| 360 |
> |
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 361 |
> |
\label{eq:SPForces} |
| 362 |
> |
\end{equation} |
| 363 |
> |
These forces are identical to the forces of the standard Coulomb |
| 364 |
> |
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 365 |
> |
\textit{et al.} as undesirable. They pointed out that the effect of |
| 366 |
> |
the image charges is neglected in the forces when this form is |
| 367 |
> |
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
| 368 |
> |
molecular dynamics. Additionally, there is a discontinuity in the |
| 369 |
> |
forces at the cutoff radius which results in energy drift during MD |
| 370 |
> |
simulations. |
| 371 |
> |
|
| 372 |
> |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 373 |
> |
will give, |
| 374 |
> |
\begin{equation} |
| 375 |
> |
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}. |
| 376 |
|
\label{eq:SFPot} |
| 377 |
|
\end{equation} |
| 378 |
< |
Equation \ref{eq:SFPot} is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term is simply equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in equation \ref{eq:ZahnPot}, resulting in a discontinuity in the potential as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 378 |
> |
with associated forces, |
| 379 |
> |
\begin{equation} |
| 380 |
> |
F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}. |
| 381 |
> |
\label{eq:SFForces} |
| 382 |
> |
\end{equation} |
| 383 |
> |
This formulation has the benefits that there are no discontinuities at |
| 384 |
> |
the cutoff radius, while the neutralizing image charges are present in |
| 385 |
> |
both the energy and force expressions. It would be simple to add the |
| 386 |
> |
self-neutralizing term back when computing the total energy of the |
| 387 |
> |
system, thereby maintaining the agreement with the Madelung energies. |
| 388 |
> |
A side effect of this treatment is the alteration in the shape of the |
| 389 |
> |
potential that comes from the derivative term. Thus, a degree of |
| 390 |
> |
clarity about agreement with the empirical potential is lost in order |
| 391 |
> |
to gain functionality in dynamics simulations. |
| 392 |
|
|
| 393 |
< |
It is important to note that shifted force techniques have a drawback in that they alter the shape of the original potential. We thereby lose a degree of clarity about the original formulation of the potential in order to gain functionality in dynamics simulations. An alternative direction would be use the derivatives of the original potential for the forces. This was addressed by Wolf \textit{et al.} as undesirable, because the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. In Monte Carlo simulations, this argument is mute, because forces are not evaluated. We decided to consider both the Shifted-Force technique described above and this Shifted-Potential technique to determine their usability in the evaluation of both energetic and dynamic results in simulations with electrostatics. |
| 393 |
> |
Wolf \textit{et al.} originally discussed the energetics of the |
| 394 |
> |
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
| 395 |
> |
insufficient for accurate determination of the energy with reasonable |
| 396 |
> |
cutoff distances. The calculated Madelung energies fluctuated around |
| 397 |
> |
the expected value as the cutoff radius was increased, but the |
| 398 |
> |
oscillations converged toward the correct value.\cite{Wolf99} A |
| 399 |
> |
damping function was incorporated to accelerate the convergence; and |
| 400 |
> |
though alternative forms for the damping function could be |
| 401 |
> |
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 402 |
> |
chosen to mirror the effective screening used in the Ewald summation. |
| 403 |
> |
Incorporating this error function damping into the simple Coulomb |
| 404 |
> |
potential, |
| 405 |
> |
\begin{equation} |
| 406 |
> |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 407 |
> |
\label{eq:dampCoulomb} |
| 408 |
> |
\end{equation} |
| 409 |
> |
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 410 |
> |
\begin{equation} |
| 411 |
> |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, |
| 412 |
> |
\label{eq:DSPPot} |
| 413 |
> |
\end{equation} |
| 414 |
> |
with associated forces, |
| 415 |
> |
\begin{equation} |
| 416 |
> |
F_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}. |
| 417 |
> |
\label{eq:DSPForces} |
| 418 |
> |
\end{equation} |
| 419 |
> |
Again, this damped shifted potential suffers from a |
| 420 |
> |
force-discontinuity at the cutoff radius, and the image charges play |
| 421 |
> |
no role in the forces. To remedy these concerns, one may derive a |
| 422 |
> |
{\sc sf} variant by including the derivative term in |
| 423 |
> |
eq. (\ref{eq:shiftingForm}), |
| 424 |
> |
\begin{equation} |
| 425 |
> |
\begin{split} |
| 426 |
> |
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}. |
| 427 |
> |
\label{eq:DSFPot} |
| 428 |
> |
\end{split} |
| 429 |
> |
\end{equation} |
| 430 |
> |
The derivative of the above potential will lead to the following forces, |
| 431 |
> |
\begin{equation} |
| 432 |
> |
\begin{split} |
| 433 |
> |
F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}. |
| 434 |
> |
\label{eq:DSFForces} |
| 435 |
> |
\end{split} |
| 436 |
> |
\end{equation} |
| 437 |
> |
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
| 438 |
> |
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
| 439 |
> |
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
| 440 |
|
|
| 441 |
< |
A variety of simulation situations were assembled and analyzed to determine the relative effectiveness of the adapted Wolf spherical truncation schemes at reproducing the results obtained using a smooth particle mesh Ewald (SPME) summation technique.\cite{Essmann95} In addition to the Shifted-Potential and Shifted-Force adapted Wolf methods, both reaction field and uncorrected cutoff methods were included for comparison purposes. The general usability of these methods in both Monte Carlo and Molecular Dynamics calculations was assessed through statistical analysis over the combined results from all of the following studied systems: |
| 442 |
< |
\begin{enumerate} |
| 443 |
< |
\item Liquid Water |
| 444 |
< |
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 445 |
< |
\item NaCl Crystal |
| 446 |
< |
\item NaCl Melt |
| 447 |
< |
\item Low Ionic Strength Solution of NaCl in Water |
| 448 |
< |
\item High Ionic Strength Solution of NaCl in Water |
| 449 |
< |
\item 6 \AA\ Radius Sphere of Argon in Water |
| 450 |
< |
\end{enumerate} |
| 451 |
< |
By studying these methods in systems composed entirely of neutral groups, charged particles, and mixtures of the two, we can either comment on possible system dependence or universal applicability of the techniques. |
| 441 |
> |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 442 |
> |
derived by Zahn \textit{et al.}; however, there are two important |
| 443 |
> |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 444 |
> |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 445 |
> |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 446 |
> |
in the Zahn potential, resulting in a potential discontinuity as |
| 447 |
> |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 448 |
> |
portion is different. The missing $v_\textrm{c}$ term would not |
| 449 |
> |
affect molecular dynamics simulations (although the computed energy |
| 450 |
> |
would be expected to have sudden jumps as particle distances crossed |
| 451 |
> |
$R_c$). The sign problem is a potential source of errors, however. |
| 452 |
> |
In fact, it introduces a discontinuity in the forces at the cutoff, |
| 453 |
> |
because the force function is shifted in the wrong direction and |
| 454 |
> |
doesn't cross zero at $R_\textrm{c}$. |
| 455 |
|
|
| 456 |
< |
\section{Methods} |
| 456 |
> |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 457 |
> |
electrostatic summation method in which the potential and forces are |
| 458 |
> |
continuous at the cutoff radius and which incorporates the damping |
| 459 |
> |
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
| 460 |
> |
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
| 461 |
> |
{\sc sf}, damping) are at reproducing the correct electrostatic |
| 462 |
> |
summation performed by the Ewald sum. |
| 463 |
|
|
| 464 |
< |
In each of the simulated systems, 500 distinct configurations were generated, and the electrostatic summation methods were compared via sequential application on each of these fixed configurations. The methods compared include SPME, the aforementioned Shifted Potential and Shifted Force methods - both with damping parameters ($\alpha$) of 0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, light, moderate, and heavy damping respectively), reaction field with an infinite dielectric constant, and an unmodified cutoff. Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation. The SPME calculations were performed using the TINKER implementation of SPME, while all all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Ponder87,Meineke05} |
| 464 |
> |
\subsection{Other alternatives} |
| 465 |
> |
In addition to the methods described above, we considered some other |
| 466 |
> |
techniques that are commonly used in molecular simulations. The |
| 467 |
> |
simplest of these is group-based cutoffs. Though of little use for |
| 468 |
> |
charged molecules, collecting atoms into neutral groups takes |
| 469 |
> |
advantage of the observation that the electrostatic interactions decay |
| 470 |
> |
faster than those for monopolar pairs.\cite{Steinbach94} When |
| 471 |
> |
considering these molecules as neutral groups, the relative |
| 472 |
> |
orientations of the molecules control the strength of the interactions |
| 473 |
> |
at the cutoff radius. Consequently, as these molecular particles move |
| 474 |
> |
through $R_\textrm{c}$, the energy will drift upward due to the |
| 475 |
> |
anisotropy of the net molecular dipole interactions.\cite{Rahman71} To |
| 476 |
> |
maintain good energy conservation, both the potential and derivative |
| 477 |
> |
need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79} |
| 478 |
> |
This is accomplished using a standard switching function. If a smooth |
| 479 |
> |
second derivative is desired, a fifth (or higher) order polynomial can |
| 480 |
> |
be used.\cite{Andrea83} |
| 481 |
|
|
| 482 |
< |
Generation of the system configurations was dependent on the system type. For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 SPC/E water molecule trajectories and each equilibrated individually. The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- ions and were selected and equilibrated in the same fashion as the water systems. For the low and high ionic strength NaCl solutions, 4 and 40 ions were first solvated in a 1000 water molecule boxes respectively. Ion and water positions were then randomly swapped, and the resulting configurations were again equilibrated individually. Finally, for the Argon/Water "charge void" systems, the identities of all the SPC/E waters within 6 \AA\ of the center of the equilibrated water configurations were converted to argon (Fig. \ref{fig:argonSlice}). |
| 482 |
> |
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
| 483 |
> |
and to incorporate the effects of the surroundings, a method like |
| 484 |
> |
Reaction Field ({\sc rf}) can be used. The original theory for {\sc |
| 485 |
> |
rf} was originally developed by Onsager,\cite{Onsager36} and it was |
| 486 |
> |
applied in simulations for the study of water by Barker and |
| 487 |
> |
Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply |
| 488 |
> |
an extension of the group-based cutoff method where the net dipole |
| 489 |
> |
within the cutoff sphere polarizes an external dielectric, which |
| 490 |
> |
reacts back on the central dipole. The same switching function |
| 491 |
> |
considerations for group-based cutoffs need to made for {\sc rf}, with |
| 492 |
> |
the additional pre-specification of a dielectric constant. |
| 493 |
|
|
| 494 |
< |
\begin{figure} |
| 112 |
< |
\centering |
| 113 |
< |
\includegraphics[width=3.25in]{./slice.pdf} |
| 114 |
< |
\caption{A slice from the center of a water box used in a charge void simulation. The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.} |
| 115 |
< |
\label{fig:argonSlice} |
| 116 |
< |
\end{figure} |
| 494 |
> |
\section{Methods} |
| 495 |
|
|
| 496 |
< |
All of these comparisons were performed with three different cutoff radii (9, 12, and 15 \AA) to investigate the cutoff radius dependence of the various techniques. It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated. Typical molecular mechanics packages default this to a value dependent on the cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller tolerances are typically associated with increased accuracy.\cite{Essmann95} The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 496 |
> |
In classical molecular mechanics simulations, there are two primary |
| 497 |
> |
techniques utilized to obtain information about the system of |
| 498 |
> |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 499 |
> |
techniques utilize pairwise summations of interactions between |
| 500 |
> |
particle sites, but they use these summations in different ways. |
| 501 |
|
|
| 502 |
< |
\section{Results and Discussion} |
| 502 |
> |
In MC, the potential energy difference between configurations dictates |
| 503 |
> |
the progression of MC sampling. Going back to the origins of this |
| 504 |
> |
method, the acceptance criterion for the canonical ensemble laid out |
| 505 |
> |
by Metropolis \textit{et al.} states that a subsequent configuration |
| 506 |
> |
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
| 507 |
> |
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
| 508 |
> |
Maintaining the correct $\Delta E$ when using an alternate method for |
| 509 |
> |
handling the long-range electrostatics will ensure proper sampling |
| 510 |
> |
from the ensemble. |
| 511 |
|
|
| 512 |
< |
\subsection{$\Delta E$ Comparison} |
| 513 |
< |
In order to evaluate the performance of the adapted Wolf Shifted Potential and Shifted Force electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations need to be compared to the results using SPME. Considering the SPME results to be the correct or desired behavior, ideal performance of a tested method is taken to be agreement between the energy differences calculated. Linear least squares regression of the $\Delta E$ values between configurations using SPME against $\Delta E$ values using tested methods provides a quantitative comparison of this agreement. Unitary results for both the correlation and correlation coefficient for these regressions indicate equivalent energetic results between the methods. The correlation is the slope of the plotted data while the correlation coefficient ($R^2$) is a measure of the of the data scatter around the fitted line and gives an idea of the quality of the fit (Fig. \ref{fig:linearFit}). |
| 512 |
> |
In MD, the derivative of the potential governs how the system will |
| 513 |
> |
progress in time. Consequently, the force and torque vectors on each |
| 514 |
> |
body in the system dictate how the system evolves. If the magnitude |
| 515 |
> |
and direction of these vectors are similar when using alternate |
| 516 |
> |
electrostatic summation techniques, the dynamics in the short term |
| 517 |
> |
will be indistinguishable. Because error in MD calculations is |
| 518 |
> |
cumulative, one should expect greater deviation at longer times, |
| 519 |
> |
although methods which have large differences in the force and torque |
| 520 |
> |
vectors will diverge from each other more rapidly. |
| 521 |
|
|
| 522 |
< |
\begin{figure} |
| 126 |
< |
\centering |
| 127 |
< |
\includegraphics[width=3.25in]{./linearFit.pdf} |
| 128 |
< |
\caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system. } |
| 129 |
< |
\label{fig:linearFit} |
| 130 |
< |
\end{figure} |
| 522 |
> |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 523 |
|
|
| 524 |
< |
With 500 independent configurations, 124,750 $\Delta E$ data points are used in a regression of a single system. Results and discussion for the individual analysis of each of the system types appear in the supporting information. To probe the applicability of each method in the general case, all the different system types were included in a single regression. The results for this regression are shown in figure \ref{fig:delE}. |
| 524 |
> |
The pairwise summation techniques (outlined in section |
| 525 |
> |
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 526 |
> |
studying the energy differences between conformations. We took the |
| 527 |
> |
{\sc spme}-computed energy difference between two conformations to be the |
| 528 |
> |
correct behavior. An ideal performance by an alternative method would |
| 529 |
> |
reproduce these energy differences exactly (even if the absolute |
| 530 |
> |
energies calculated by the methods are different). Since none of the |
| 531 |
> |
methods provide exact energy differences, we used linear least squares |
| 532 |
> |
regressions of energy gap data to evaluate how closely the methods |
| 533 |
> |
mimicked the Ewald energy gaps. Unitary results for both the |
| 534 |
> |
correlation (slope) and correlation coefficient for these regressions |
| 535 |
> |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 536 |
> |
Sample correlation plots for two alternate methods are shown in |
| 537 |
> |
Fig. \ref{fig:linearFit}. |
| 538 |
|
|
| 539 |
|
\begin{figure} |
| 540 |
|
\centering |
| 541 |
< |
\includegraphics[width=3.25in]{./delEplot.pdf} |
| 542 |
< |
\caption{The results from the statistical analysis of the $\Delta$E results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii. Results close to a value of 1 (dashed line) indicate $\Delta E$ values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME. Reaction Field results do not include NaCl crystal or melt configurations.} |
| 543 |
< |
\label{fig:delE} |
| 541 |
> |
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
| 542 |
> |
\caption{Example least squares regressions of the configuration energy |
| 543 |
> |
differences for SPC/E water systems. The upper plot shows a data set |
| 544 |
> |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 545 |
> |
shows a data set with a good correlation coefficient.} |
| 546 |
> |
\label{fig:linearFit} |
| 547 |
|
\end{figure} |
| 548 |
|
|
| 549 |
< |
In figure \ref{fig:delE}, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff. This is not all that surprising since this results in large energy fluctuations as atoms move in and out of the cutoff radius. These fluctuations can be alleviated to some degree by using group based cutoffs with a switching function. The Group Switch Cutoff row doesn't show a significant improvement in this plot because the salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral. Correcting the resulting charged cutoff sphere is one of the purposes of the shifted potential proposed by Wolf \textit{et al.}, and this correction indeed improves the results as seen in the Shifted Potental rows. While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME. Inclusion of potential damping improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows an excellent correlation and quality of fit with the SPME results, particularly with a cutoff radius greater than 12 \AA . Use of a larger damping parameter is more helpful for the shortest cutoff shown, but it has a detrimental effect on simulations with larger cutoffs. In the Shifted Force sets, increasing damping results in progressively poorer correlation. Overall, the undamped case is the best performing set, as the correlation and quality of fits are consistently superior regardless of the cutoff distance. This result is beneficial in that the undamped case is less computationally prohibitive do to the lack of complimentary error function calculation when performing the electrostatic pair interaction. The reaction field results illustrates some of that method's limitations, primarily that it was developed for use in homogenous systems; although it does provide results that are an improvement over those from an unmodified cutoff. |
| 549 |
> |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 550 |
> |
were represented using 500 independent configurations. Thus, each of |
| 551 |
> |
the alternative (non-Ewald) electrostatic summation methods was |
| 552 |
> |
evaluated using an accumulated 873,250 configurational energy |
| 553 |
> |
differences. |
| 554 |
|
|
| 555 |
< |
\subsection{Force Magnitude Comparison} |
| 555 |
> |
Results and discussion for the individual analysis of each of the |
| 556 |
> |
system types appear in the supporting information,\cite{EPAPSdeposit} |
| 557 |
> |
while the cumulative results over all the investigated systems appears |
| 558 |
> |
below in section \ref{sec:EnergyResults}. |
| 559 |
|
|
| 560 |
< |
While studying the energy differences provides insight into how comparable these methods are energetically, if we want to use these methods in Molecular Dynamics simulations, we also need to consider their effect on forces and torques. Both the magnitude and the direction of the force and torque vectors of each of the bodies in the system can be compared to those observed while using SPME. Analysis of the magnitude of these vectors can be performed in the manner described previously for comparing $\Delta E$ values, only instead of a single value between two system configurations, there is a value for each particle in each configuration. For a system of 1000 water molecules and 40 ions, there are 1040 force vectors and 1000 torque vectors. With 500 configurations, this results in excess of 500,000 data samples for each system type. Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the force and torque vector magnitude results for the accumulated analysis over all the system types. |
| 560 |
> |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 561 |
> |
We evaluated the pairwise methods (outlined in section |
| 562 |
> |
\ref{sec:ESMethods}) for use in MD simulations by |
| 563 |
> |
comparing the force and torque vectors with those obtained using the |
| 564 |
> |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 565 |
> |
direction of these vectors on each of the bodies in the system were |
| 566 |
> |
analyzed. For the magnitude of these vectors, linear least squares |
| 567 |
> |
regression analyses were performed as described previously for |
| 568 |
> |
comparing $\Delta E$ values. Instead of a single energy difference |
| 569 |
> |
between two system configurations, we compared the magnitudes of the |
| 570 |
> |
forces (and torques) on each molecule in each configuration. For a |
| 571 |
> |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 572 |
> |
vectors and 1000 torque vectors. With 500 configurations, this |
| 573 |
> |
results in 520,000 force and 500,000 torque vector comparisons. |
| 574 |
> |
Additionally, data from seven different system types was aggregated |
| 575 |
> |
before the comparison was made. |
| 576 |
> |
|
| 577 |
> |
The {\it directionality} of the force and torque vectors was |
| 578 |
> |
investigated through measurement of the angle ($\theta$) formed |
| 579 |
> |
between those computed from the particular method and those from {\sc spme}, |
| 580 |
> |
\begin{equation} |
| 581 |
> |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 582 |
> |
\end{equation} |
| 583 |
> |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 584 |
> |
vector computed using method M. Each of these $\theta$ values was |
| 585 |
> |
accumulated in a distribution function and weighted by the area on the |
| 586 |
> |
unit sphere. Since this distribution is a measure of angular error |
| 587 |
> |
between two different electrostatic summation methods, there is no |
| 588 |
> |
{\it a priori} reason for the profile to adhere to any specific |
| 589 |
> |
shape. Thus, gaussian fits were used to measure the width of the |
| 590 |
> |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 591 |
> |
each of these fits and was used to compare distribution widths. |
| 592 |
> |
Values of $\sigma^2$ near zero indicate vector directions |
| 593 |
> |
indistinguishable from those calculated when using the reference |
| 594 |
> |
method ({\sc spme}). |
| 595 |
|
|
| 596 |
+ |
\subsection{Short-time Dynamics} |
| 597 |
+ |
|
| 598 |
+ |
The effects of the alternative electrostatic summation methods on the |
| 599 |
+ |
short-time dynamics of charged systems were evaluated by considering a |
| 600 |
+ |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 601 |
+ |
performing pairwise methods was used in this comparison. The NaCl |
| 602 |
+ |
crystal was chosen to avoid possible complications from the treatment |
| 603 |
+ |
of orientational motion in molecular systems. All systems were |
| 604 |
+ |
started with the same initial positions and velocities. Simulations |
| 605 |
+ |
were performed under the microcanonical ensemble, and velocity |
| 606 |
+ |
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 607 |
+ |
of the trajectories, |
| 608 |
+ |
\begin{equation} |
| 609 |
+ |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 610 |
+ |
\label{eq:vCorr} |
| 611 |
+ |
\end{equation} |
| 612 |
+ |
Velocity autocorrelation functions require detailed short time data, |
| 613 |
+ |
thus velocity information was saved every 2 fs over 10 ps |
| 614 |
+ |
trajectories. Because the NaCl crystal is composed of two different |
| 615 |
+ |
atom types, the average of the two resulting velocity autocorrelation |
| 616 |
+ |
functions was used for comparisons. |
| 617 |
+ |
|
| 618 |
+ |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 619 |
+ |
|
| 620 |
+ |
The effects of the same subset of alternative electrostatic methods on |
| 621 |
+ |
the {\it long-time} dynamics of charged systems were evaluated using |
| 622 |
+ |
the same model system (NaCl crystals at 1000~K). The power spectrum |
| 623 |
+ |
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 624 |
+ |
autocorrelation function, \begin{equation} I(\omega) = |
| 625 |
+ |
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 626 |
+ |
\label{eq:powerSpec} |
| 627 |
+ |
\end{equation} |
| 628 |
+ |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 629 |
+ |
NaCl crystal is composed of two different atom types, the average of |
| 630 |
+ |
the two resulting power spectra was used for comparisons. Simulations |
| 631 |
+ |
were performed under the microcanonical ensemble, and velocity |
| 632 |
+ |
information was saved every 5~fs over 100~ps trajectories. |
| 633 |
+ |
|
| 634 |
+ |
\subsection{Representative Simulations}\label{sec:RepSims} |
| 635 |
+ |
A variety of representative molecular simulations were analyzed to |
| 636 |
+ |
determine the relative effectiveness of the pairwise summation |
| 637 |
+ |
techniques in reproducing the energetics and dynamics exhibited by |
| 638 |
+ |
{\sc spme}. We wanted to span the space of typical molecular |
| 639 |
+ |
simulations (i.e. from liquids of neutral molecules to ionic |
| 640 |
+ |
crystals), so the systems studied were: |
| 641 |
+ |
\begin{enumerate} |
| 642 |
+ |
\item liquid water (SPC/E),\cite{Berendsen87} |
| 643 |
+ |
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 644 |
+ |
\item NaCl crystals, |
| 645 |
+ |
\item NaCl melts, |
| 646 |
+ |
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 647 |
+ |
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 648 |
+ |
\item a 6 \AA\ radius sphere of Argon in water. |
| 649 |
+ |
\end{enumerate} |
| 650 |
+ |
By utilizing the pairwise techniques (outlined in section |
| 651 |
+ |
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 652 |
+ |
charged particles, and mixtures of the two, we hope to discern under |
| 653 |
+ |
which conditions it will be possible to use one of the alternative |
| 654 |
+ |
summation methodologies instead of the Ewald sum. |
| 655 |
+ |
|
| 656 |
+ |
For the solid and liquid water configurations, configurations were |
| 657 |
+ |
taken at regular intervals from high temperature trajectories of 1000 |
| 658 |
+ |
SPC/E water molecules. Each configuration was equilibrated |
| 659 |
+ |
independently at a lower temperature (300~K for the liquid, 200~K for |
| 660 |
+ |
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 661 |
+ |
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 662 |
+ |
these systems were selected and equilibrated in the same manner as the |
| 663 |
+ |
water systems. In order to introduce measurable fluctuations in the |
| 664 |
+ |
configuration energy differences, the crystalline simulations were |
| 665 |
+ |
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 666 |
+ |
NaCl configurations needed to represent a fully disordered array of |
| 667 |
+ |
point charges, so the high temperature of 7000~K was selected for |
| 668 |
+ |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 669 |
+ |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 670 |
+ |
water positions were then randomly swapped, and the resulting |
| 671 |
+ |
configurations were again equilibrated individually. Finally, for the |
| 672 |
+ |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 673 |
+ |
waters within 6 \AA\ of the center of the equilibrated water |
| 674 |
+ |
configurations were converted to argon. |
| 675 |
+ |
|
| 676 |
+ |
These procedures guaranteed us a set of representative configurations |
| 677 |
+ |
from chemically-relevant systems sampled from appropriate |
| 678 |
+ |
ensembles. Force field parameters for the ions and Argon were taken |
| 679 |
+ |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 680 |
+ |
|
| 681 |
+ |
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 682 |
+ |
We compared the following alternative summation methods with results |
| 683 |
+ |
from the reference method ({\sc spme}): |
| 684 |
+ |
\begin{itemize} |
| 685 |
+ |
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 686 |
+ |
and 0.3 \AA$^{-1}$, |
| 687 |
+ |
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 688 |
+ |
and 0.3 \AA$^{-1}$, |
| 689 |
+ |
\item reaction field with an infinite dielectric constant, and |
| 690 |
+ |
\item an unmodified cutoff. |
| 691 |
+ |
\end{itemize} |
| 692 |
+ |
Group-based cutoffs with a fifth-order polynomial switching function |
| 693 |
+ |
were utilized for the reaction field simulations. Additionally, we |
| 694 |
+ |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 695 |
+ |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 696 |
+ |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 697 |
+ |
were performed using the {\sc oopse} molecular mechanics |
| 698 |
+ |
package.\cite{Meineke05} All other portions of the energy calculation |
| 699 |
+ |
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 700 |
+ |
manner across all systems and configurations. |
| 701 |
+ |
|
| 702 |
+ |
The alternative methods were also evaluated with three different |
| 703 |
+ |
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 704 |
+ |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 705 |
+ |
real-space and reciprocal-space portions of the Ewald calculation. |
| 706 |
+ |
Typical molecular mechanics packages set this to a value dependent on |
| 707 |
+ |
the cutoff radius and a tolerance (typically less than $1 \times |
| 708 |
+ |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 709 |
+ |
increasing accuracy at the expense of computational time spent on the |
| 710 |
+ |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 711 |
+ |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 712 |
+ |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 713 |
+ |
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 714 |
+ |
respectively. |
| 715 |
+ |
|
| 716 |
+ |
\section{Results and Discussion} |
| 717 |
+ |
|
| 718 |
+ |
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 719 |
+ |
In order to evaluate the performance of the pairwise electrostatic |
| 720 |
+ |
summation methods for Monte Carlo simulations, the energy differences |
| 721 |
+ |
between configurations were compared to the values obtained when using |
| 722 |
+ |
{\sc spme}. The results for the subsequent regression analysis are shown in |
| 723 |
+ |
figure \ref{fig:delE}. |
| 724 |
+ |
|
| 725 |
|
\begin{figure} |
| 726 |
|
\centering |
| 727 |
< |
\includegraphics[width=3.25in]{./frcMagplot.pdf} |
| 728 |
< |
\caption{The results from the statistical analysis of the force vector magnitude results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii. Results close to a value of 1 (dashed line) indicate force vector magnitude values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME.} |
| 729 |
< |
\label{fig:frcMag} |
| 727 |
> |
\includegraphics[width=5.5in]{./delEplot.pdf} |
| 728 |
> |
\caption{Statistical analysis of the quality of configurational energy |
| 729 |
> |
differences for a given electrostatic method compared with the |
| 730 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 731 |
> |
indicate $\Delta E$ values indistinguishable from those obtained using |
| 732 |
> |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 733 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 734 |
> |
inverted triangles).} |
| 735 |
> |
\label{fig:delE} |
| 736 |
|
\end{figure} |
| 737 |
|
|
| 738 |
< |
The results in figure \ref{fig:frcMag} for the most part parallel those seen in the previous look at the $\Delta E$ results. The unmodified cutoff results are poor, but using group based cutoffs and a switching function provides a improvement much more significant than what was seen with $\Delta E$. Looking at the Shifted Potential sets, the slope and $R^2$ improve with the use of damping to an optimal result of 0.2 \AA $^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, while beneficial for simulations with a cutoff radius of 9 \AA\ , is detrimental to simulations with larger cutoff radii. The undamped Shifted Force method gives forces in line with those obtained using SPME, and use of a damping function results in minor improvement. The reaction field results are surprisingly good, considering the poor quality of the fits for the $\Delta E$ results. There is still a considerable degree of scatter in the data, but it correlates well in general. |
| 738 |
> |
The most striking feature of this plot is how well the Shifted Force |
| 739 |
> |
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 740 |
> |
differences. For the undamped {\sc sf} method, and the |
| 741 |
> |
moderately-damped {\sc sp} methods, the results are nearly |
| 742 |
> |
indistinguishable from the Ewald results. The other common methods do |
| 743 |
> |
significantly less well. |
| 744 |
|
|
| 745 |
< |
\subsection{Torque Magnitude Comparison} |
| 745 |
> |
The unmodified cutoff method is essentially unusable. This is not |
| 746 |
> |
surprising since hard cutoffs give large energy fluctuations as atoms |
| 747 |
> |
or molecules move in and out of the cutoff |
| 748 |
> |
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
| 749 |
> |
some degree by using group based cutoffs with a switching |
| 750 |
> |
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 751 |
> |
significant improvement using the group-switched cutoff because the |
| 752 |
> |
salt and salt solution systems contain non-neutral groups. Interested |
| 753 |
> |
readers can consult the accompanying supporting information for a |
| 754 |
> |
comparison where all groups are neutral.\cite{EPAPSdeposit} |
| 755 |
|
|
| 756 |
+ |
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 757 |
+ |
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 758 |
+ |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 759 |
+ |
results, particularly with a cutoff radius greater than 12 |
| 760 |
+ |
\AA . Use of a larger damping parameter is more helpful for the |
| 761 |
+ |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 762 |
+ |
with larger cutoffs. |
| 763 |
+ |
|
| 764 |
+ |
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 765 |
+ |
worse} correlation with Ewald. Overall, the undamped case is the best |
| 766 |
+ |
performing set, as the correlation and quality of fits are |
| 767 |
+ |
consistently superior regardless of the cutoff distance. The undamped |
| 768 |
+ |
case is also less computationally demanding (because no evaluation of |
| 769 |
+ |
the complementary error function is required). |
| 770 |
+ |
|
| 771 |
+ |
The reaction field results illustrates some of that method's |
| 772 |
+ |
limitations, primarily that it was developed for use in homogenous |
| 773 |
+ |
systems; although it does provide results that are an improvement over |
| 774 |
+ |
those from an unmodified cutoff. |
| 775 |
+ |
|
| 776 |
+ |
\subsection{Magnitudes of the Force and Torque Vectors} |
| 777 |
+ |
|
| 778 |
+ |
Evaluation of pairwise methods for use in Molecular Dynamics |
| 779 |
+ |
simulations requires consideration of effects on the forces and |
| 780 |
+ |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
| 781 |
+ |
regression results for the force and torque vector magnitudes, |
| 782 |
+ |
respectively. The data in these figures was generated from an |
| 783 |
+ |
accumulation of the statistics from all of the system types. |
| 784 |
+ |
|
| 785 |
|
\begin{figure} |
| 786 |
|
\centering |
| 787 |
< |
\includegraphics[width=3.25in]{./trqMagplot.pdf} |
| 788 |
< |
\caption{The results from the statistical analysis of the torque vector magnitude results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii. Results close to a value of 1 (dashed line) indicate torque vector magnitude values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME. Torques are only accumulated on the rigid water molecules, so these results exclude NaCl the systems.} |
| 789 |
< |
\label{fig:trqMag} |
| 787 |
> |
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 788 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 789 |
> |
magnitudes for a given electrostatic method compared with the |
| 790 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 791 |
> |
indicate force magnitude values indistinguishable from those obtained |
| 792 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 793 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 794 |
> |
inverted triangles).} |
| 795 |
> |
\label{fig:frcMag} |
| 796 |
|
\end{figure} |
| 797 |
|
|
| 798 |
< |
The torque vector magnitude results in figure \ref{fig:trqMag} are similar to those seen for the forces, but more clearly show the improved behavior with increasing cutoff radius. Moderate damping is beneficial to the Shifted Potential and unnecessary with the Shifted Force method, and they also show that over-damping adversely effects all cutoff radii rather than showing an improvement for systems with short cutoffs. The reaction field method performs well when calculating the torques, better than the Shifted Force method over this limited data set. |
| 798 |
> |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 799 |
> |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 800 |
> |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 801 |
> |
This is perhaps expected because this method explicitly incorporates a |
| 802 |
> |
smooth transition in the forces at the cutoff radius as well as the |
| 803 |
> |
neutralizing image charges. |
| 804 |
|
|
| 805 |
< |
\subsection{Force and Torque Direction Comparison} |
| 805 |
> |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 806 |
> |
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 807 |
> |
poor, but using group based cutoffs and a switching function provides |
| 808 |
> |
an improvement much more significant than what was seen with $\Delta |
| 809 |
> |
E$. |
| 810 |
|
|
| 811 |
< |
Having force and torque vectors with magnitudes that are well correlated to SPME is good, but if they are not pointing in the proper direction the results will be incorrect. These vector directions were investigated through measurement of the angle formed between them and those from SPME. The dot product of these unit vectors provides a theta value that is accumulated in a distribution function, weighted by the area on the unit sphere. Narrow distributions of theta values indicates similar to identical results between the tested method and SPME. To measure the narrowness of the resulting distributions, non-linear Gaussian fits were performed. |
| 811 |
> |
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 812 |
> |
method is generating usable forces. Further increases in damping, |
| 813 |
> |
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 814 |
> |
detrimental to simulations with larger cutoff radii. |
| 815 |
|
|
| 816 |
+ |
The reaction field results are surprisingly good, considering the poor |
| 817 |
+ |
quality of the fits for the $\Delta E$ results. There is still a |
| 818 |
+ |
considerable degree of scatter in the data, but the forces correlate |
| 819 |
+ |
well with the Ewald forces in general. We note that the reaction |
| 820 |
+ |
field calculations do not include the pure NaCl systems, so these |
| 821 |
+ |
results are partly biased towards conditions in which the method |
| 822 |
+ |
performs more favorably. |
| 823 |
+ |
|
| 824 |
|
\begin{figure} |
| 825 |
|
\centering |
| 826 |
< |
\includegraphics[width=3.25in]{./gaussFit.pdf} |
| 827 |
< |
\caption{Example fitting of the angular distribution of the force vectors over all of the studied systems. The solid and dotted lines show Gaussian and Voigt fits of the distribution data respectively. Even though the Voigt profile make for a more accurate fit, the Gaussian was used due to more versatile statistical results.} |
| 828 |
< |
\label{fig:gaussian} |
| 826 |
> |
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 827 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 828 |
> |
magnitudes for a given electrostatic method compared with the |
| 829 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 830 |
> |
indicate torque magnitude values indistinguishable from those obtained |
| 831 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 832 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 833 |
> |
inverted triangles).} |
| 834 |
> |
\label{fig:trqMag} |
| 835 |
|
\end{figure} |
| 836 |
|
|
| 837 |
< |
Figure \ref{fig:gaussian} shows an example distribution and the non-linear fit applied. The solid line is a Gaussian profile, while the dotted line is a Voigt profile, a convolution of a Gaussian and a Lorentzian profile. Since this distribution is a measure of angular error between two different electrostatic summation methods, there is particular reason for it to adhere to a particular shape. Because of this and the Gaussian profile's more statistically meaningful properties, Gaussian fitting was used to compare all the methods considered in this study. The results (Fig. \ref{fig:frcTrqAng}) are compared through the variance ($\sigma^2$) of these non-linear fits. |
| 837 |
> |
Molecular torques were only available from the systems which contained |
| 838 |
> |
rigid molecules (i.e. the systems containing water). The data in |
| 839 |
> |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 840 |
> |
|
| 841 |
> |
Torques appear to be much more sensitive to charges at a longer |
| 842 |
> |
distance. The striking feature in comparing the new electrostatic |
| 843 |
> |
methods with {\sc spme} is how much the agreement improves with increasing |
| 844 |
> |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 845 |
> |
appears to be reproducing the {\sc spme} torques most accurately. |
| 846 |
> |
|
| 847 |
> |
Water molecules are dipolar, and the reaction field method reproduces |
| 848 |
> |
the effect of the surrounding polarized medium on each of the |
| 849 |
> |
molecular bodies. Therefore it is not surprising that reaction field |
| 850 |
> |
performs best of all of the methods on molecular torques. |
| 851 |
> |
|
| 852 |
> |
\subsection{Directionality of the Force and Torque Vectors} |
| 853 |
> |
|
| 854 |
> |
It is clearly important that a new electrostatic method can reproduce |
| 855 |
> |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 856 |
> |
sum. However, the {\it directionality} of these vectors will also be |
| 857 |
> |
vital in calculating dynamical quantities accurately. Force and |
| 858 |
> |
torque directionalities were investigated by measuring the angles |
| 859 |
> |
formed between these vectors and the same vectors calculated using |
| 860 |
> |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 861 |
> |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 862 |
> |
distributions of the combined set over all system types. |
| 863 |
|
|
| 864 |
|
\begin{figure} |
| 865 |
|
\centering |
| 866 |
< |
\includegraphics[width=3.25in]{./frcTrqAngplot.pdf} |
| 867 |
< |
\caption{The results from the statistical analysis of the force and torque vector angular distributions for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii. Plotted values are the variance ($\sigma^2$) of the Gaussian non-linear fits. Results close to a value of 0 (dashed line) indicate force or torque vector directions from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME. Torques are only accumulated on the rigid water molecules, so the torque vector angle results exclude NaCl the systems.} |
| 866 |
> |
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 867 |
> |
\caption{Statistical analysis of the width of the angular distribution |
| 868 |
> |
that the force and torque vectors from a given electrostatic method |
| 869 |
> |
make with their counterparts obtained using the reference Ewald sum. |
| 870 |
> |
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 871 |
> |
indicate force and torque directions indistinguishable from those |
| 872 |
> |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 873 |
> |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 874 |
> |
and 15\AA\ = inverted triangles).} |
| 875 |
|
\label{fig:frcTrqAng} |
| 876 |
|
\end{figure} |
| 877 |
|
|
| 878 |
< |
Both the force and torque $\sigma^2$ results from the analysis of the total accumulated system data are tabulated in figure \ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case show the improvement afforded by choosing a longer simulation cutoff. Increasing the cutoff from 9 to 12 \AA\ typically results in a halving of $\sigma^2$, with a similar improvement going from 12 to 15 \AA . The undamped Shifted Force, Group Based Cutoff, and Reaction Field methods all do equivalently well at capturing the direction of both the force and torque vectors. Using damping improves the angular behavior significantly for the Shifted Potential and moderately for the Shifted Force methods. Increasing the damping too far is destructive for both methods, particularly to the torque vectors. Again it is important to recognize that the force vectors cover all particles in the systems, while torque vectors are only available for neutral molecular groups. Damping appears to have a more beneficial effect on non-neutral bodies, and this observation is investigated further in the accompanying supporting information. |
| 878 |
> |
Both the force and torque $\sigma^2$ results from the analysis of the |
| 879 |
> |
total accumulated system data are tabulated in figure |
| 880 |
> |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 881 |
> |
sp}) method would be essentially unusable for molecular dynamics |
| 882 |
> |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 883 |
> |
method, however, is generating force and torque vectors which are |
| 884 |
> |
within a few degrees of the Ewald results even with weak (or no) |
| 885 |
> |
damping. |
| 886 |
|
|
| 887 |
+ |
All of the sets (aside from the over-damped case) show the improvement |
| 888 |
+ |
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 889 |
+ |
from 9 to 12 \AA\ typically results in a halving of the width of the |
| 890 |
+ |
distribution, with a similar improvement when going from 12 to 15 |
| 891 |
+ |
\AA . |
| 892 |
+ |
|
| 893 |
+ |
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 894 |
+ |
all do equivalently well at capturing the direction of both the force |
| 895 |
+ |
and torque vectors. Using the electrostatic damping improves the |
| 896 |
+ |
angular behavior significantly for the {\sc sp} and moderately for the |
| 897 |
+ |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 898 |
+ |
it is important to recognize that the force vectors cover all |
| 899 |
+ |
particles in all seven systems, while torque vectors are only |
| 900 |
+ |
available for neutral molecular groups. Damping is more beneficial to |
| 901 |
+ |
charged bodies, and this observation is investigated further in the |
| 902 |
+ |
accompanying supporting information.\cite{EPAPSdeposit} |
| 903 |
+ |
|
| 904 |
+ |
Although not discussed previously, group based cutoffs can be applied |
| 905 |
+ |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 906 |
+ |
will reintroduce small discontinuities at the cutoff radius, but the |
| 907 |
+ |
effects of these can be minimized by utilizing a switching function. |
| 908 |
+ |
Though there are no significant benefits or drawbacks observed in |
| 909 |
+ |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 910 |
+ |
is a measurable improvement in the directionality of the forces and |
| 911 |
+ |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 912 |
+ |
obtained using group based cutoffs along with the results seen in |
| 913 |
+ |
figure \ref{fig:frcTrqAng}. The {\sc sp} (with an $\alpha$ of 0.2 |
| 914 |
+ |
\AA$^{-1}$ or smaller) shows much narrower angular distributions when |
| 915 |
+ |
using group-based cutoffs. The {\sc sf} method likewise shows |
| 916 |
+ |
improvement in the undamped and lightly damped cases. |
| 917 |
+ |
|
| 918 |
|
\begin{table}[htbp] |
| 919 |
< |
\centering |
| 920 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 919 |
> |
\centering |
| 920 |
> |
\caption{Statistical analysis of the angular |
| 921 |
> |
distributions that the force (upper) and torque (lower) vectors |
| 922 |
> |
from a given electrostatic method make with their counterparts |
| 923 |
> |
obtained using the reference Ewald sum. Calculations were |
| 924 |
> |
performed both with (Y) and without (N) group based cutoffs and a |
| 925 |
> |
switching function. The $\alpha$ values have units of \AA$^{-1}$ |
| 926 |
> |
and the variance values have units of degrees$^2$.} |
| 927 |
> |
|
| 928 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 929 |
|
\\ |
| 930 |
|
\toprule |
| 955 |
|
\label{tab:groupAngle} |
| 956 |
|
\end{table} |
| 957 |
|
|
| 958 |
< |
Although not discussed previously, group based cutoffs can be applied to both the Shifted Potential and Force methods. Use off a switching function corrects for the discontinuities that arise when atoms of a group exit the cutoff before the group's center of mass. Though there are no significant benefit or drawbacks observed in $\Delta E$ and vector magnitude results when doing this, there is a measurable improvement in the vector angle results. Table \ref{tab:groupAngle} shows the angular variance values obtained using group based cutoffs and a switching function alongside the standard results seen in figure \ref{fig:frcTrqAng} for comparison purposes. The Shifted Potential shows much narrower angular distributions for both the force and torque vectors when using an $\alpha$ of 0.2 \AA$^{-1}$ or less, while Shifted Force shows improvements in the undamped and lightly damped cases. Thus, by calculating the electrostatic interactions in terms of molecular pairs rather than atomic pairs, the direction of the force and torque vectors are determined more accurately. |
| 958 |
> |
One additional trend in table \ref{tab:groupAngle} is that the |
| 959 |
> |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 960 |
> |
increases, something that is more obvious with group-based cutoffs. |
| 961 |
> |
The complimentary error function inserted into the potential weakens |
| 962 |
> |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 963 |
> |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 964 |
> |
electrostatic interaction and to remove it completely. Kast |
| 965 |
> |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 966 |
> |
values for these types of electrostatic summation methods by fitting |
| 967 |
> |
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 968 |
> |
0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ |
| 969 |
> |
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 970 |
> |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 971 |
> |
these findings, choices this high would introduce error in the |
| 972 |
> |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 973 |
> |
observations, empirical damping up to 0.2 \AA$^{-1}$ is beneficial, |
| 974 |
> |
but damping may be unnecessary when using the {\sc sf} method. |
| 975 |
|
|
| 976 |
< |
One additional trend to recognize in table \ref{tab:groupAngle} is that the $\sigma^2$ values for both Shifted Potential and Shifted Force converge as $\alpha$ increases, something that is easier to see when using group based cutoffs. Looking back on figures \ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this behavior clearly at large $\alpha$ and cutoff values. The reason for this is that the complimentary error function inserted into the potential weakens the electrostatic interaction as $\alpha$ increases. Thus, at larger values of $\alpha$, both the summation method types progress toward non-interacting functions, so care is required in choosing large damping functions lest one generate an undesirable loss in the pair interaction. Kast \textit{et al.} developed a method for choosing appropriate $\alpha$ values for these types of electrostatic summation methods by fitting to $g(r)$ data, and their methods indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear to be reasonable choices to obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on these findings, choices this high would introduce error in the molecular torques, particularly for the shorter cutoffs. Based on the above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be beneficial, but is arguably unnecessary when using the Shifted-Force method. |
| 976 |
> |
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 977 |
|
|
| 978 |
< |
\subsection{Crystal Power Spectrum} |
| 978 |
> |
Zahn {\it et al.} investigated the structure and dynamics of water |
| 979 |
> |
using eqs. (\ref{eq:ZahnPot}) and |
| 980 |
> |
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 981 |
> |
that a method similar (but not identical with) the damped {\sc sf} |
| 982 |
> |
method resulted in properties very similar to those obtained when |
| 983 |
> |
using the Ewald summation. The properties they studied (pair |
| 984 |
> |
distribution functions, diffusion constants, and velocity and |
| 985 |
> |
orientational correlation functions) may not be particularly sensitive |
| 986 |
> |
to the long-range and collective behavior that governs the |
| 987 |
> |
low-frequency behavior in crystalline systems. Additionally, the |
| 988 |
> |
ionic crystals are the worst case scenario for the pairwise methods |
| 989 |
> |
because they lack the reciprocal space contribution contained in the |
| 990 |
> |
Ewald summation. |
| 991 |
|
|
| 992 |
< |
In the previous studies using a Shifted-Force variant of the damped Wolf coulomb potential, the structure and dynamics of water were investigated rather extensively.\cite{Zahn02,Kast03} Their results indicated that the damped Shifted-Force method results in properties very similar to those obtained when using the Ewald summation. Considering the statistical results shown above, the good performance of this method is not that surprising. Rather than consider the same systems and simply recapitulate their results, we decided to look at the solid state dynamical behavior obtained using the best performing summation methods from the above results. |
| 992 |
> |
We are using two separate measures to probe the effects of these |
| 993 |
> |
alternative electrostatic methods on the dynamics in crystalline |
| 994 |
> |
materials. For short- and intermediate-time dynamics, we are |
| 995 |
> |
computing the velocity autocorrelation function, and for long-time |
| 996 |
> |
and large length-scale collective motions, we are looking at the |
| 997 |
> |
low-frequency portion of the power spectrum. |
| 998 |
|
|
| 999 |
< |
Using the NaCl crystal as the model system, trajectories were obtained using SPME; Shifted-Force with $\alpha$ values of 0, 0.1 and 0.2 \AA$^{-1}$; and Shifted-Potential with an $\alpha$ value of 0.2 \AA$^{-1}$. To enhance the atomic motion, these simulations were run at 1000 K, near the experimental $T_m$ for NaCl. The velocity autocorrelation function (Eq. \ref{eq:vCorr})was computed on each of the trajectories. |
| 1000 |
< |
\begin{equation} |
| 1001 |
< |
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 1002 |
< |
\label{eq:vCorr} |
| 1003 |
< |
\end{equation} |
| 1004 |
< |
Velocity autocorrelation functions require detailed short time data and long trajectories for good statistics, thus velocity information was saved every 5 fs over 100 ps trajectories. The power spectrum ($I(\omega)$) is obtained via discrete Fourier transform of the autocorrelation function |
| 1005 |
< |
\begin{equation} |
| 1006 |
< |
I(\omega) = \sum^{N-1}_{\omega=0}C_v(t)e^{-i\omega t/N}, |
| 1007 |
< |
\label{eq:powerSpec} |
| 1008 |
< |
\end{equation} |
| 1009 |
< |
where $N$ is the number of time samples in $C_v(t)$ and the frequency, $\omega=0,\ 1,\ ...,\ N-1$. The resulting spectra (Fig. \ref{fig:normalModes}) show the normal mode frequencies for the crystal under the simulated conditions. |
| 999 |
> |
\begin{figure} |
| 1000 |
> |
\centering |
| 1001 |
> |
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 1002 |
> |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1003 |
> |
1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1004 |
> |
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 1005 |
> |
the first minimum. The times to first collision are nearly identical, |
| 1006 |
> |
but differences can be seen in the peaks and troughs, where the |
| 1007 |
> |
undamped and weakly damped methods are stiffer than the moderately |
| 1008 |
> |
damped and {\sc spme} methods.} |
| 1009 |
> |
\label{fig:vCorrPlot} |
| 1010 |
> |
\end{figure} |
| 1011 |
|
|
| 1012 |
+ |
The short-time decay of the velocity autocorrelation function through |
| 1013 |
+ |
the first collision are nearly identical in figure |
| 1014 |
+ |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1015 |
+ |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1016 |
+ |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1017 |
+ |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1018 |
+ |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1019 |
+ |
approaches the {\sc spme} results. With $\alpha$ values of 0.2 \AA$^{-1}$, |
| 1020 |
+ |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1021 |
+ |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1022 |
+ |
and {\sc sp} potentials become nearly identical with increased |
| 1023 |
+ |
damping. However, this appears to indicate that once damping is |
| 1024 |
+ |
utilized, the details of the form of the potential (and forces) |
| 1025 |
+ |
constructed out of the damped electrostatic interaction are less |
| 1026 |
+ |
important. |
| 1027 |
+ |
|
| 1028 |
+ |
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 1029 |
+ |
|
| 1030 |
+ |
To evaluate how the differences between the methods affect the |
| 1031 |
+ |
collective long-time motion, we computed power spectra from long-time |
| 1032 |
+ |
traces of the velocity autocorrelation function. The power spectra for |
| 1033 |
+ |
the best-performing alternative methods are shown in |
| 1034 |
+ |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1035 |
+ |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1036 |
+ |
ringing resulting from data truncation. This procedure had no |
| 1037 |
+ |
noticeable effect on peak location or magnitude. |
| 1038 |
+ |
|
| 1039 |
|
\begin{figure} |
| 1040 |
|
\centering |
| 1041 |
< |
\includegraphics[width = 3.25in]{./nModeFTPlotDot.pdf} |
| 1042 |
< |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, Shifted-Force ($\alpha$ = 0, 0.1, \& 0.2), and Shifted-Potential ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differentiate.} |
| 1043 |
< |
\label{fig:normalModes} |
| 1041 |
> |
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1042 |
> |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1043 |
> |
functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf} |
| 1044 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
| 1045 |
> |
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
| 1046 |
> |
spectra differ.} |
| 1047 |
> |
\label{fig:methodPS} |
| 1048 |
|
\end{figure} |
| 1049 |
|
|
| 1050 |
< |
Figure \ref{fig:normalModes} shows the power spectra for the NaCl crystals (from averaged Na and Cl ion velocity autocorrelation functions) using the stated electrostatic summation methods. While high frequency peaks of all the spectra overlap, showing the same general features, the low frequency region shows how the summation methods differ. The normal modes at frequencies below 100 cm$^{-1}$ are shifted up when using undamped or weakly damped Shifted-Force. When using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential methods give near identical normal mode behavior as the Ewald method (which has a damping value of 0.3119). The damping acts as a distance dependent Gaussian screening of the point charges in the system. This weakening of the electrostatic interaction with distance explains why the low level normal modes are at lower frequencies for the moderately damped methods than for undamped or weakly damped methods. Consider damping on a simple real-space electrostatic potential in the form |
| 1051 |
< |
\begin{equation} |
| 1052 |
< |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r_{ij}})}{r_{ij}}\right]S(r), |
| 1053 |
< |
\end{equation} |
| 1054 |
< |
where $S(r)$ is a switching function that smoothly zeroes the potential at the cutoff radius. Figure \ref{fig:dampInc} shows how the low frequency normal modes are dependent on the damping used in the direct electrostatic sum. As the damping increases, the normal modes drop to lower frequencies. Incidentally, use of an $\alpha$ of 0.25 \AA$^{-1}$ on a simple electrostatic summation results in low frequency normal mode dynamics equivalent to a simulation using SPME. When the coefficient lowers to 0.15 \AA$^{-1}$ and below, the normal modes shift to higher frequency in exponential fashion. Though not shown, the spectrum for the simple undamped electrostatic potential is blue-shifted such that the lowest normal mode resides near 325 cm$^{-1}$. In light of these results, the undamped Shifted-Force method producing low-lying normal modes within 10 cm$^{-1}$ of SPME is quite respectable; however, it appears as though moderate damping is required for accurate reproduction of crystal dynamics. |
| 1050 |
> |
While the high frequency regions of the power spectra for the |
| 1051 |
> |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1052 |
> |
the low frequency region shows how the summation methods differ. |
| 1053 |
> |
Considering the low-frequency inset (expanded in the upper frame of |
| 1054 |
> |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1055 |
> |
correlated motions are blue-shifted when using undamped or weakly |
| 1056 |
> |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1057 |
> |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1058 |
> |
correlated motion to the Ewald method (which has a convergence |
| 1059 |
> |
parameter of 0.3119 \AA$^{-1}$). This weakening of the electrostatic |
| 1060 |
> |
interaction with increased damping explains why the long-ranged |
| 1061 |
> |
correlated motions are at lower frequencies for the moderately damped |
| 1062 |
> |
methods than for undamped or weakly damped methods. |
| 1063 |
> |
|
| 1064 |
> |
To isolate the role of the damping constant, we have computed the |
| 1065 |
> |
spectra for a single method ({\sc sf}) with a range of damping |
| 1066 |
> |
constants and compared this with the {\sc spme} spectrum. |
| 1067 |
> |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1068 |
> |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1069 |
> |
However, even without any electrostatic damping, the {\sc sf} method |
| 1070 |
> |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1071 |
> |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1072 |
> |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1073 |
> |
Most} of the collective behavior in the crystal is accurately captured |
| 1074 |
> |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1075 |
> |
obtained using moderate damping in addition to the shifting at the |
| 1076 |
> |
cutoff distance. |
| 1077 |
> |
|
| 1078 |
|
\begin{figure} |
| 1079 |
|
\centering |
| 1080 |
< |
\includegraphics[width = 3.25in]{./alphaCompare.pdf} |
| 1081 |
< |
\caption{Normal modes for an NaCl crystal at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$)ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the normal modes are red-shifted towards and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 1080 |
> |
\includegraphics[width = \linewidth]{./increasedDamping.pdf} |
| 1081 |
> |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1082 |
> |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
| 1083 |
> |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1084 |
> |
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
| 1085 |
> |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1086 |
> |
result in underestimates of frequencies of the long-wavelength |
| 1087 |
> |
motions.} |
| 1088 |
|
\label{fig:dampInc} |
| 1089 |
|
\end{figure} |
| 1090 |
|
|
| 1091 |
|
\section{Conclusions} |
| 1092 |
|
|
| 1093 |
< |
This investigation of pairwise electrostatic summation techniques shows that there are viable and more computationally efficient electrostatic summation techniques than the Ewald summation, chiefly methods derived from the damped Coulombic sum originally proposed by Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the Shifted-Force method, reformulated above, shows a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. The cumulative energy difference results showed the undamped Shifted-Force and moderately damped Shifted-Potential methods produced results nearly identical to SPME. Similarly for the dynamic features, the un- to moderately damped Shifted-Force and moderately damped Shifted-Potential methods produce force and torque vector magnitude and directions very similar to the expected values. These results translate into long-time dynamic behavior equivalent to that produced in simulations using SPME. |
| 1093 |
> |
This investigation of pairwise electrostatic summation techniques |
| 1094 |
> |
shows that there are viable and computationally efficient alternatives |
| 1095 |
> |
to the Ewald summation. These methods are derived from the damped and |
| 1096 |
> |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 1097 |
> |
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
| 1098 |
> |
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
| 1099 |
> |
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1100 |
> |
energetic and dynamic characteristics exhibited by simulations |
| 1101 |
> |
employing lattice summation techniques. The cumulative energy |
| 1102 |
> |
difference results showed the undamped {\sc sf} and moderately damped |
| 1103 |
> |
{\sc sp} methods produced results nearly identical to {\sc spme}. Similarly |
| 1104 |
> |
for the dynamic features, the undamped or moderately damped {\sc sf} |
| 1105 |
> |
and moderately damped {\sc sp} methods produce force and torque vector |
| 1106 |
> |
magnitude and directions very similar to the expected values. These |
| 1107 |
> |
results translate into long-time dynamic behavior equivalent to that |
| 1108 |
> |
produced in simulations using {\sc spme}. |
| 1109 |
|
|
| 1110 |
< |
Aside from the computational cost benefit, these techniques have applicability in situations where the use of the Ewald sum can prove problematic. Primary among them is their use in interfacial systems, where the unmodified lattice sum techniques artificially accentuate the periodicity of the system in an undesirable manner. There have been alterations to the standard Ewald techniques, via corrections and reformulations, to compensate for these systems; but these pairwise techniques require no modifications, making them natural tools to tackle these problems. Additionally, this transferability gives it benefits over other pairwise methods, like reaction field, because estimations of physical properties, like the dielectric constant, are unnecessary. |
| 1110 |
> |
As in all purely-pairwise cutoff methods, these methods are expected |
| 1111 |
> |
to scale approximately {\it linearly} with system size, and they are |
| 1112 |
> |
easily parallelizable. This should result in substantial reductions |
| 1113 |
> |
in the computational cost of performing large simulations. |
| 1114 |
|
|
| 1115 |
< |
These results don't deprecate the use of the Ewald summation; in fact, it is the standard to which these simple pairwise sums are judged. However, these results do speak to the necessity of the Ewald summation in all molecular simulations. That a simple pairwise technique can be substituted to gain nearly all the physical effects provided by the full lattice sum makes us question whether the minimal perturbations bestowed through added complexity and increased cost are worth it. |
| 1115 |
> |
Aside from the computational cost benefit, these techniques have |
| 1116 |
> |
applicability in situations where the use of the Ewald sum can prove |
| 1117 |
> |
problematic. Of greatest interest is their potential use in |
| 1118 |
> |
interfacial systems, where the unmodified lattice sum techniques |
| 1119 |
> |
artificially accentuate the periodicity of the system in an |
| 1120 |
> |
undesirable manner. There have been alterations to the standard Ewald |
| 1121 |
> |
techniques, via corrections and reformulations, to compensate for |
| 1122 |
> |
these systems; but the pairwise techniques discussed here require no |
| 1123 |
> |
modifications, making them natural tools to tackle these problems. |
| 1124 |
> |
Additionally, this transferability gives them benefits over other |
| 1125 |
> |
pairwise methods, like reaction field, because estimations of physical |
| 1126 |
> |
properties (e.g. the dielectric constant) are unnecessary. |
| 1127 |
|
|
| 1128 |
+ |
If a researcher is using Monte Carlo simulations of large chemical |
| 1129 |
+ |
systems containing point charges, most structural features will be |
| 1130 |
+ |
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1131 |
+ |
method with an electrostatic damping of 0.2 \AA$^{-1}$. These methods |
| 1132 |
+ |
would also be appropriate for molecular dynamics simulations where the |
| 1133 |
+ |
data of interest is either structural or short-time dynamical |
| 1134 |
+ |
quantities. For long-time dynamics and collective motions, the safest |
| 1135 |
+ |
pairwise method we have evaluated is the {\sc sf} method with an |
| 1136 |
+ |
electrostatic damping between 0.2 and 0.25 |
| 1137 |
+ |
\AA$^{-1}$. |
| 1138 |
+ |
|
| 1139 |
+ |
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1140 |
+ |
fact, it is the standard by which these simple pairwise sums have been |
| 1141 |
+ |
judged. However, these results do suggest that in the typical |
| 1142 |
+ |
simulations performed today, the Ewald summation may no longer be |
| 1143 |
+ |
required to obtain the level of accuracy most researchers have come to |
| 1144 |
+ |
expect. |
| 1145 |
+ |
|
| 1146 |
|
\section{Acknowledgments} |
| 1147 |
+ |
Support for this project was provided by the National Science |
| 1148 |
+ |
Foundation under grant CHE-0134881. The authors would like to thank |
| 1149 |
+ |
Steve Corcelli and Ed Maginn for helpful discussions and comments. |
| 1150 |
|
|
| 1151 |
|
\newpage |
| 1152 |
|
|
| 1153 |
< |
\bibliographystyle{achemso} |
| 1153 |
> |
\bibliographystyle{jcp2} |
| 1154 |
|
\bibliography{electrostaticMethods} |
| 1155 |
|
|
| 1156 |
|
|
