| 65 |
|
\section{Introduction} |
| 66 |
|
|
| 67 |
|
In molecular simulations, proper accumulation of the electrostatic |
| 68 |
< |
interactions is considered one of the most essential and |
| 69 |
< |
computationally demanding tasks. The common molecular mechanics force |
| 70 |
< |
fields are founded on representation of the atomic sites centered on |
| 71 |
< |
full or partial charges shielded by Lennard-Jones type interactions. |
| 72 |
< |
This means that nearly every pair interaction involves an |
| 73 |
< |
charge-charge calculation. Coupled with $r^{-1}$ decay, the monopole |
| 74 |
< |
interactions quickly become a burden for molecular systems of all |
| 75 |
< |
sizes. For example, in small systems, the electrostatic pair |
| 76 |
< |
interaction may not have decayed appreciably within the box length |
| 77 |
< |
leading to an effect excluded from the pair interactions within a unit |
| 78 |
< |
box. In large systems, excessively large cutoffs need to be used to |
| 79 |
< |
accurately incorporate their effect, and since the computational cost |
| 80 |
< |
increases proportionally with the cutoff sphere, it quickly becomes |
| 81 |
< |
very time-consuming to perform these calculations. |
| 68 |
> |
interactions is essential and is one of the most |
| 69 |
> |
computationally-demanding tasks. The common molecular mechanics force |
| 70 |
> |
fields represent atomic sites with full or partial charges protected |
| 71 |
> |
by Lennard-Jones (short range) interactions. This means that nearly |
| 72 |
> |
every pair interaction involves a calculation of charge-charge forces. |
| 73 |
> |
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
| 74 |
> |
interactions quickly become the most expensive part of molecular |
| 75 |
> |
simulations. Historically, the electrostatic pair interaction would |
| 76 |
> |
not have decayed appreciably within the typical box lengths that could |
| 77 |
> |
be feasibly simulated. In the larger systems that are more typical of |
| 78 |
> |
modern simulations, large cutoffs should be used to incorporate |
| 79 |
> |
electrostatics correctly. |
| 80 |
|
|
| 81 |
< |
There have been many efforts to address this issue of both proper and |
| 82 |
< |
practical handling of electrostatic interactions, and these have |
| 83 |
< |
resulted in the availability of a variety of |
| 84 |
< |
techniques.\cite{Roux99,Sagui99,Tobias01} These are typically |
| 85 |
< |
classified as implicit methods (i.e., continuum dielectrics, static |
| 86 |
< |
dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e., |
| 87 |
< |
Ewald summations, interaction shifting or |
| 90 |
< |
trucation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 81 |
> |
There have been many efforts to address the proper and practical |
| 82 |
> |
handling of electrostatic interactions, and these have resulted in a |
| 83 |
> |
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are |
| 84 |
> |
typically classified as implicit methods (i.e., continuum dielectrics, |
| 85 |
> |
static dipolar fields),\cite{Born20,Grossfield00} explicit methods |
| 86 |
> |
(i.e., Ewald summations, interaction shifting or |
| 87 |
> |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 88 |
|
reaction field type methods, fast multipole |
| 89 |
|
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 90 |
< |
often preferred because they incorporate dynamic solvent molecules in |
| 91 |
< |
the system of interest, but these methods are sometimes difficult to |
| 92 |
< |
utilize because of their high computational cost.\cite{Roux99} In |
| 93 |
< |
addition to this cost, there has been some question of the inherent |
| 94 |
< |
periodicity of the explicit Ewald summation artificially influencing |
| 95 |
< |
systems dynamics.\cite{Tobias01} |
| 90 |
> |
often preferred because they physically incorporate solvent molecules |
| 91 |
> |
in the system of interest, but these methods are sometimes difficult |
| 92 |
> |
to utilize because of their high computational cost.\cite{Roux99} In |
| 93 |
> |
addition to the computational cost, there have been some questions |
| 94 |
> |
regarding possible artifacts caused by the inherent periodicity of the |
| 95 |
> |
explicit Ewald summation.\cite{Tobias01} |
| 96 |
|
|
| 97 |
< |
In this paper, we focus on the common mixed and explicit methods of |
| 98 |
< |
reaction filed and smooth particle mesh |
| 99 |
< |
Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods |
| 100 |
< |
devised by Wolf {\it et al.} which we further extend.\cite{Wolf99} |
| 101 |
< |
These new methods for handling electrostatics are quite |
| 102 |
< |
computationally efficient, since they involve only a simple |
| 103 |
< |
modification to the direct pairwise sum, and they lack the added |
| 104 |
< |
periodicity of the Ewald sum. Below, these methods are evaluated using |
| 105 |
< |
a variety of model systems and comparison methodologies to establish |
| 106 |
< |
their useability in molecular simulations. |
| 97 |
> |
In this paper, we focus on a new set of shifted methods devised by |
| 98 |
> |
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 99 |
> |
methods along with a few other mixed methods (i.e. reaction field) are |
| 100 |
> |
compared with the smooth particle mesh Ewald |
| 101 |
> |
sum,\cite{Onsager36,Essmann99} which is our reference method for |
| 102 |
> |
handling long-range electrostatic interactions. The new methods for |
| 103 |
> |
handling electrostatics have the potential to scale linearly with |
| 104 |
> |
increasing system size since they involve only a simple modification |
| 105 |
> |
to the direct pairwise sum. They also lack the added periodicity of |
| 106 |
> |
the Ewald sum, so they can be used for systems which are non-periodic |
| 107 |
> |
or which have one- or two-dimensional periodicity. Below, these |
| 108 |
> |
methods are evaluated using a variety of model systems to establish |
| 109 |
> |
their usability in molecular simulations. |
| 110 |
|
|
| 111 |
|
\subsection{The Ewald Sum} |
| 112 |
|
The complete accumulation electrostatic interactions in a system with |
| 113 |
|
periodic boundary conditions (PBC) requires the consideration of the |
| 114 |
< |
effect of all charges within a simulation box, as well as those in the |
| 115 |
< |
periodic replicas, |
| 114 |
> |
effect of all charges within a (cubic) simulation box as well as those |
| 115 |
> |
in the periodic replicas, |
| 116 |
|
\begin{equation} |
| 117 |
|
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], |
| 118 |
|
\label{eq:PBCSum} |
| 123 |
|
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
| 124 |
|
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
| 125 |
|
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
| 126 |
< |
$j$, and $\phi$ is Poisson's equation ($\phi(\mathbf{r}_{ij}) = q_i |
| 127 |
< |
q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge interactions). In the |
| 128 |
< |
case of monopole electrostatics, eq. (\ref{eq:PBCSum}) is |
| 129 |
< |
conditionally convergent and is discontiuous for non-neutral systems. |
| 126 |
> |
$j$, and $\phi$ is the solution to Poisson's equation |
| 127 |
> |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for |
| 128 |
> |
charge-charge interactions). In the case of monopole electrostatics, |
| 129 |
> |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for |
| 130 |
> |
non-neutral systems. |
| 131 |
|
|
| 132 |
< |
This electrostatic summation problem was originally studied by Ewald |
| 132 |
> |
The electrostatic summation problem was originally studied by Ewald |
| 133 |
|
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 134 |
|
took was to convert this conditionally convergent sum into two |
| 135 |
|
absolutely convergent summations: a short-ranged real-space summation |
| 141 |
|
\label{eq:EwaldSum} |
| 142 |
|
\end{equation} |
| 143 |
|
where $\alpha$ is a damping parameter, or separation constant, with |
| 144 |
< |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and equal |
| 145 |
< |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 146 |
< |
constant of the encompassing medium. The final two terms of |
| 144 |
> |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are |
| 145 |
> |
equal to $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the |
| 146 |
> |
dielectric constant of the surrounding medium. The final two terms of |
| 147 |
|
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 148 |
|
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 149 |
|
dipolar term was neglected in early applications in molecular |
| 150 |
|
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 151 |
|
Leeuw {\it et al.} to address situations where the unit cell has a |
| 152 |
< |
dipole moment and this dipole moment gets magnified through |
| 153 |
< |
replication of the periodic images.\cite{deLeeuw80,Smith81} If this |
| 154 |
< |
term is taken to be zero, the system is using conducting boundary |
| 152 |
> |
dipole moment which is magnified through replication of the periodic |
| 153 |
> |
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
| 154 |
> |
system is said to be using conducting (or ``tin-foil'') boundary |
| 155 |
|
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 156 |
|
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 157 |
< |
time. Initially, due to the small size of systems, the entire |
| 158 |
< |
simulation box was replicated to convergence. Currently, we balance a |
| 159 |
< |
spherical real-space cutoff with the reciprocal sum and consider the |
| 160 |
< |
surrounding dielectric. |
| 157 |
> |
time. Initially, due to the small sizes of the systems that could be |
| 158 |
> |
feasibly simulated, the entire simulation box was replicated to |
| 159 |
> |
convergence. In more modern simulations, the simulation boxes have |
| 160 |
> |
grown large enough that a real-space cutoff could potentially give |
| 161 |
> |
convergent behavior. Indeed, it has often been observed that the |
| 162 |
> |
reciprocal-space portion of the Ewald sum can be vanishingly |
| 163 |
> |
small compared to the real-space portion.\cite{XXX} |
| 164 |
> |
|
| 165 |
|
\begin{figure} |
| 166 |
|
\centering |
| 167 |
|
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 175 |
|
\label{fig:ewaldTime} |
| 176 |
|
\end{figure} |
| 177 |
|
|
| 178 |
< |
The Ewald summation in the straight-forward form is an |
| 179 |
< |
$\mathscr{O}(N^2)$ algorithm. The separation constant $(\alpha)$ |
| 180 |
< |
plays an important role in the computational cost balance between the |
| 181 |
< |
direct and reciprocal-space portions of the summation. The choice of |
| 182 |
< |
the magnitude of this value allows one to select whether the |
| 183 |
< |
real-space or reciprocal space portion of the summation is an |
| 184 |
< |
$\mathscr{O}(N^2)$ calcualtion (with the other being |
| 185 |
< |
$\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$ |
| 186 |
< |
and thoughtful algorithm development, this cost can be brought down to |
| 187 |
< |
$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to |
| 188 |
< |
reduce the cost of the Ewald summation further is to set $\alpha$ such |
| 189 |
< |
that the real-space interactions decay rapidly, allowing for a short |
| 190 |
< |
spherical cutoff, and then optimize the reciprocal space summation. |
| 191 |
< |
These optimizations usually involve the utilization of the fast |
| 187 |
< |
Fourier transform (FFT),\cite{Hockney81} leading to the |
| 178 |
> |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 179 |
> |
separation constant $(\alpha)$ plays an important role in balancing |
| 180 |
> |
the computational cost between the direct and reciprocal-space |
| 181 |
> |
portions of the summation. The choice of this value allows one to |
| 182 |
> |
select whether the real-space or reciprocal space portion of the |
| 183 |
> |
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
| 184 |
> |
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
| 185 |
> |
$\alpha$ and thoughtful algorithm development, this cost can be |
| 186 |
> |
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
| 187 |
> |
taken to reduce the cost of the Ewald summation even further is to set |
| 188 |
> |
$\alpha$ such that the real-space interactions decay rapidly, allowing |
| 189 |
> |
for a short spherical cutoff. Then the reciprocal space summation is |
| 190 |
> |
optimized. These optimizations usually involve utilization of the |
| 191 |
> |
fast Fourier transform (FFT),\cite{Hockney81} leading to the |
| 192 |
|
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 193 |
|
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 194 |
|
methods, the cost of the reciprocal-space portion of the Ewald |
| 195 |
< |
summation is from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$. |
| 195 |
> |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
| 196 |
> |
\log N)$. |
| 197 |
|
|
| 198 |
< |
These developments and optimizations have led the use of the Ewald |
| 199 |
< |
summation to become routine in simulations with periodic boundary |
| 200 |
< |
conditions. However, in certain systems the intrinsic three |
| 201 |
< |
dimensional periodicity can prove to be problematic, such as two |
| 202 |
< |
dimensional surfaces and membranes. The Ewald sum has been |
| 203 |
< |
reformulated to handle 2D |
| 204 |
< |
systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new |
| 205 |
< |
methods have been found to be computationally |
| 206 |
< |
expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the |
| 207 |
< |
full Ewald summation is a possible direction for enabling the handling |
| 203 |
< |
of 2D systems and the inclusion of the optimizations described |
| 204 |
< |
previously.\cite{Yeh99} |
| 198 |
> |
These developments and optimizations have made the use of the Ewald |
| 199 |
> |
summation routine in simulations with periodic boundary |
| 200 |
> |
conditions. However, in certain systems, such as vapor-liquid |
| 201 |
> |
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 202 |
> |
can prove problematic. The Ewald sum has been reformulated to handle |
| 203 |
> |
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the |
| 204 |
> |
new methods are computationally expensive.\cite{Spohr97,Yeh99} |
| 205 |
> |
Inclusion of a correction term in the Ewald summation is a possible |
| 206 |
> |
direction for handling 2D systems while still enabling the use of the |
| 207 |
> |
modern optimizations.\cite{Yeh99} |
| 208 |
|
|
| 209 |
|
Several studies have recognized that the inherent periodicity in the |
| 210 |
< |
Ewald sum can also have an effect on systems that have the same |
| 211 |
< |
dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 212 |
< |
Good examples are solvated proteins kept at high relative |
| 213 |
< |
concentration due to the periodicity of the electrostatics. In these |
| 210 |
> |
Ewald sum can also have an effect on three-dimensional |
| 211 |
> |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 212 |
> |
Solvated proteins are essentially kept at high concentration due to |
| 213 |
> |
the periodicity of the electrostatic summation method. In these |
| 214 |
|
systems, the more compact folded states of a protein can be |
| 215 |
|
artificially stabilized by the periodic replicas introduced by the |
| 216 |
< |
Ewald summation.\cite{Weber00} Thus, care ought to be taken when |
| 217 |
< |
considering the use of the Ewald summation where the intrinsic |
| 218 |
< |
perodicity may negatively affect the system dynamics. |
| 216 |
> |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
| 217 |
> |
considering the use of the Ewald summation where the assumed |
| 218 |
> |
periodicity would introduce spurious effects in the system dynamics. |
| 219 |
|
|
| 217 |
– |
|
| 220 |
|
\subsection{The Wolf and Zahn Methods} |
| 221 |
|
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 222 |
|
for the accurate accumulation of electrostatic interactions in an |
| 223 |
< |
efficient pairwise fashion and lacks the inherent periodicity of the |
| 224 |
< |
Ewald summation.\cite{Wolf99} Wolf \textit{et al.} observed that the |
| 225 |
< |
electrostatic interaction is effectively short-ranged in condensed |
| 226 |
< |
phase systems and that neutralization of the charge contained within |
| 227 |
< |
the cutoff radius is crucial for potential stability. They devised a |
| 228 |
< |
pairwise summation method that ensures charge neutrality and gives |
| 229 |
< |
results similar to those obtained with the Ewald summation. The |
| 230 |
< |
resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes |
| 231 |
< |
image-charges subtracted out through placement on the cutoff sphere |
| 232 |
< |
and a distance-dependent damping function (identical to that seen in |
| 233 |
< |
the real-space portion of the Ewald sum) to aid convergence |
| 223 |
> |
efficient pairwise fashion. This procedure lacks the inherent |
| 224 |
> |
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
| 225 |
> |
observed that the electrostatic interaction is effectively |
| 226 |
> |
short-ranged in condensed phase systems and that neutralization of the |
| 227 |
> |
charge contained within the cutoff radius is crucial for potential |
| 228 |
> |
stability. They devised a pairwise summation method that ensures |
| 229 |
> |
charge neutrality and gives results similar to those obtained with the |
| 230 |
> |
Ewald summation. The resulting shifted Coulomb potential |
| 231 |
> |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 232 |
> |
placement on the cutoff sphere and a distance-dependent damping |
| 233 |
> |
function (identical to that seen in the real-space portion of the |
| 234 |
> |
Ewald sum) to aid convergence |
| 235 |
|
\begin{equation} |
| 236 |
< |
V_{\textrm{Wolf}}(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\}. |
| 236 |
> |
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\}. |
| 237 |
|
\label{eq:WolfPot} |
| 238 |
|
\end{equation} |
| 239 |
|
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 258 |
|
force expressions for use in simulations involving water.\cite{Zahn02} |
| 259 |
|
In their work, they pointed out that the forces and derivative of |
| 260 |
|
the potential are not commensurate. Attempts to use both |
| 261 |
< |
Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 261 |
> |
eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 262 |
|
to poor energy conservation. They correctly observed that taking the |
| 263 |
|
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 264 |
|
derivatives gives forces for a different potential energy function |
| 265 |
< |
than the one shown in Eq. (\ref{eq:WolfPot}). |
| 265 |
> |
than the one shown in eq. (\ref{eq:WolfPot}). |
| 266 |
|
|
| 267 |
< |
Zahn \textit{et al.} proposed a modified form of this ``Wolf summation |
| 268 |
< |
method'' as a way to use this technique in Molecular Dynamics |
| 269 |
< |
simulations. Taking the integral of the forces shown in equation |
| 267 |
< |
\ref{eq:WolfForces}, they proposed a new damped Coulomb |
| 268 |
< |
potential, |
| 267 |
> |
Zahn \textit{et al.} introduced a modified form of this summation |
| 268 |
> |
method as a way to use the technique in Molecular Dynamics |
| 269 |
> |
simulations. They proposed a new damped Coulomb potential, |
| 270 |
|
\begin{equation} |
| 271 |
< |
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\}. |
| 271 |
> |
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\}, |
| 272 |
|
\label{eq:ZahnPot} |
| 273 |
|
\end{equation} |
| 274 |
< |
They showed that this potential does fairly well at capturing the |
| 274 |
> |
and showed that this potential does fairly well at capturing the |
| 275 |
|
structural and dynamic properties of water compared the same |
| 276 |
|
properties obtained using the Ewald sum. |
| 277 |
|
|
| 302 |
|
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 303 |
|
shifted potential, |
| 304 |
|
\begin{equation} |
| 305 |
< |
v_\textrm{SP}(r) = \begin{cases} |
| 305 |
> |
V_\textrm{SP}(r) = \begin{cases} |
| 306 |
|
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 307 |
|
R_\textrm{c} |
| 308 |
|
\end{cases}, |
| 310 |
|
\end{equation} |
| 311 |
|
and shifted force, |
| 312 |
|
\begin{equation} |
| 313 |
< |
v_\textrm{SF}(r) = \begin{cases} |
| 313 |
> |
V_\textrm{SF}(r) = \begin{cases} |
| 314 |
|
v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 315 |
|
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 316 |
|
\end{cases}, |
| 326 |
|
The forces associated with the shifted potential are simply the forces |
| 327 |
|
of the unshifted potential itself (when inside the cutoff sphere), |
| 328 |
|
\begin{equation} |
| 329 |
< |
f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 329 |
> |
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 330 |
|
\end{equation} |
| 331 |
|
and are zero outside. Inside the cutoff sphere, the forces associated |
| 332 |
|
with the shifted force form can be written, |
| 333 |
|
\begin{equation} |
| 334 |
< |
f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 334 |
> |
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 335 |
|
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 336 |
|
\end{equation} |
| 337 |
|
|
| 338 |
< |
If the potential ($v(r)$) is taken to be the normal Coulomb potential, |
| 338 |
> |
If the potential, $v(r)$, is taken to be the normal Coulomb potential, |
| 339 |
|
\begin{equation} |
| 340 |
|
v(r) = \frac{q_i q_j}{r}, |
| 341 |
|
\label{eq:Coulomb} |
| 343 |
|
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 344 |
|
al.}'s undamped prescription: |
| 345 |
|
\begin{equation} |
| 346 |
< |
v_\textrm{SP}(r) = |
| 346 |
> |
V_\textrm{SP}(r) = |
| 347 |
|
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 348 |
|
r\leqslant R_\textrm{c}, |
| 349 |
|
\label{eq:SPPot} |
| 350 |
|
\end{equation} |
| 351 |
|
with associated forces, |
| 352 |
|
\begin{equation} |
| 353 |
< |
f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 353 |
> |
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 354 |
|
\label{eq:SPForces} |
| 355 |
|
\end{equation} |
| 356 |
|
These forces are identical to the forces of the standard Coulomb |
| 365 |
|
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 366 |
|
will give, |
| 367 |
|
\begin{equation} |
| 368 |
< |
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}. |
| 368 |
> |
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}. |
| 369 |
|
\label{eq:SFPot} |
| 370 |
|
\end{equation} |
| 371 |
|
with associated forces, |
| 372 |
|
\begin{equation} |
| 373 |
< |
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}. |
| 373 |
> |
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}. |
| 374 |
|
\label{eq:SFForces} |
| 375 |
|
\end{equation} |
| 376 |
|
This formulation has the benefits that there are no discontinuities at |
| 377 |
< |
the cutoff distance, while the neutralizing image charges are present |
| 378 |
< |
in both the energy and force expressions. It would be simple to add |
| 379 |
< |
the self-neutralizing term back when computing the total energy of the |
| 377 |
> |
the cutoff radius, while the neutralizing image charges are present in |
| 378 |
> |
both the energy and force expressions. It would be simple to add the |
| 379 |
> |
self-neutralizing term back when computing the total energy of the |
| 380 |
|
system, thereby maintaining the agreement with the Madelung energies. |
| 381 |
|
A side effect of this treatment is the alteration in the shape of the |
| 382 |
|
potential that comes from the derivative term. Thus, a degree of |
| 384 |
|
to gain functionality in dynamics simulations. |
| 385 |
|
|
| 386 |
|
Wolf \textit{et al.} originally discussed the energetics of the |
| 387 |
< |
shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that |
| 388 |
< |
it was still insufficient for accurate determination of the energy |
| 389 |
< |
with reasonable cutoff distances. The calculated Madelung energies |
| 390 |
< |
fluctuate around the expected value with increasing cutoff radius, but |
| 391 |
< |
the oscillations converge toward the correct value.\cite{Wolf99} A |
| 387 |
> |
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
| 388 |
> |
insufficient for accurate determination of the energy with reasonable |
| 389 |
> |
cutoff distances. The calculated Madelung energies fluctuated around |
| 390 |
> |
the expected value as the cutoff radius was increased, but the |
| 391 |
> |
oscillations converged toward the correct value.\cite{Wolf99} A |
| 392 |
|
damping function was incorporated to accelerate the convergence; and |
| 393 |
< |
though alternative functional forms could be |
| 393 |
> |
though alternative forms for the damping function could be |
| 394 |
|
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 395 |
|
chosen to mirror the effective screening used in the Ewald summation. |
| 396 |
|
Incorporating this error function damping into the simple Coulomb |
| 399 |
|
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 400 |
|
\label{eq:dampCoulomb} |
| 401 |
|
\end{equation} |
| 402 |
< |
the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using |
| 402 |
< |
eq. (\ref{eq:shiftingForm}), |
| 402 |
> |
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 403 |
|
\begin{equation} |
| 404 |
< |
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}, |
| 404 |
> |
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}, |
| 405 |
|
\label{eq:DSPPot} |
| 406 |
|
\end{equation} |
| 407 |
|
with associated forces, |
| 408 |
|
\begin{equation} |
| 409 |
< |
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}. |
| 409 |
> |
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}. |
| 410 |
|
\label{eq:DSPForces} |
| 411 |
|
\end{equation} |
| 412 |
< |
Again, this damped shifted potential suffers from a discontinuity and |
| 413 |
< |
a lack of the image charges in the forces. To remedy these concerns, |
| 414 |
< |
one may derive a {\sc sf} variant by including the derivative |
| 415 |
< |
term in eq. (\ref{eq:shiftingForm}), |
| 412 |
> |
Again, this damped shifted potential suffers from a |
| 413 |
> |
force-discontinuity at the cutoff radius, and the image charges play |
| 414 |
> |
no role in the forces. To remedy these concerns, one may derive a |
| 415 |
> |
{\sc sf} variant by including the derivative term in |
| 416 |
> |
eq. (\ref{eq:shiftingForm}), |
| 417 |
|
\begin{equation} |
| 418 |
|
\begin{split} |
| 419 |
< |
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}. |
| 419 |
> |
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}. |
| 420 |
|
\label{eq:DSFPot} |
| 421 |
|
\end{split} |
| 422 |
|
\end{equation} |
| 423 |
|
The derivative of the above potential will lead to the following forces, |
| 424 |
|
\begin{equation} |
| 425 |
|
\begin{split} |
| 426 |
< |
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}. |
| 426 |
> |
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}. |
| 427 |
|
\label{eq:DSFForces} |
| 428 |
|
\end{split} |
| 429 |
|
\end{equation} |
| 430 |
< |
If the damping parameter $(\alpha)$ is chosen to be zero, the undamped |
| 431 |
< |
case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered |
| 432 |
< |
from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}). |
| 430 |
> |
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
| 431 |
> |
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
| 432 |
> |
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
| 433 |
|
|
| 434 |
|
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 435 |
|
derived by Zahn \textit{et al.}; however, there are two important |
| 441 |
|
portion is different. The missing $v_\textrm{c}$ term would not |
| 442 |
|
affect molecular dynamics simulations (although the computed energy |
| 443 |
|
would be expected to have sudden jumps as particle distances crossed |
| 444 |
< |
$R_c$). The sign problem would be a potential source of errors, |
| 445 |
< |
however. In fact, it introduces a discontinuity in the forces at the |
| 446 |
< |
cutoff, because the force function is shifted in the wrong direction |
| 447 |
< |
and doesn't cross zero at $R_\textrm{c}$. |
| 444 |
> |
$R_c$). The sign problem is a potential source of errors, however. |
| 445 |
> |
In fact, it introduces a discontinuity in the forces at the cutoff, |
| 446 |
> |
because the force function is shifted in the wrong direction and |
| 447 |
> |
doesn't cross zero at $R_\textrm{c}$. |
| 448 |
|
|
| 449 |
|
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 450 |
< |
electrostatic summation method that is continuous in both the |
| 451 |
< |
potential and forces and which incorporates the damping function |
| 452 |
< |
proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of this |
| 453 |
< |
paper, we will evaluate exactly how good these methods ({\sc sp}, {\sc |
| 454 |
< |
sf}, damping) are at reproducing the correct electrostatic summation |
| 455 |
< |
performed by the Ewald sum. |
| 450 |
> |
electrostatic summation method in which the potential and forces are |
| 451 |
> |
continuous at the cutoff radius and which incorporates the damping |
| 452 |
> |
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
| 453 |
> |
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
| 454 |
> |
{\sc sf}, damping) are at reproducing the correct electrostatic |
| 455 |
> |
summation performed by the Ewald sum. |
| 456 |
|
|
| 457 |
|
\subsection{Other alternatives} |
| 458 |
< |
In addition to the methods described above, we will consider some |
| 459 |
< |
other techniques that commonly get used in molecular simulations. The |
| 458 |
> |
In addition to the methods described above, we considered some other |
| 459 |
> |
techniques that are commonly used in molecular simulations. The |
| 460 |
|
simplest of these is group-based cutoffs. Though of little use for |
| 461 |
< |
non-neutral molecules, collecting atoms into neutral groups takes |
| 461 |
> |
charged molecules, collecting atoms into neutral groups takes |
| 462 |
|
advantage of the observation that the electrostatic interactions decay |
| 463 |
|
faster than those for monopolar pairs.\cite{Steinbach94} When |
| 464 |
< |
considering these molecules as groups, an orientational aspect is |
| 465 |
< |
introduced to the interactions. Consequently, as these molecular |
| 466 |
< |
particles move through $R_\textrm{c}$, the energy will drift upward |
| 467 |
< |
due to the anisotropy of the net molecular dipole |
| 468 |
< |
interactions.\cite{Rahman71} To maintain good energy conservation, |
| 469 |
< |
both the potential and derivative need to be smoothly switched to zero |
| 470 |
< |
at $R_\textrm{c}$.\cite{Adams79} This is accomplished using a |
| 471 |
< |
switching function, |
| 472 |
< |
\begin{equation} |
| 473 |
< |
S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\ |
| 473 |
< |
\frac{(R_\textrm{c}+2r-3r_\textrm{sw})(R_\textrm{c}-r)^2}{(R_\textrm{c}-r_\textrm{sw})^3} &\quad r_\textrm{sw}<r\leqslant R_\textrm{c} \\ |
| 474 |
< |
0 &\quad r>R_\textrm{c} |
| 475 |
< |
\end{cases}, |
| 476 |
< |
\end{equation} |
| 477 |
< |
where the above form is for a cubic function. If a smooth second |
| 478 |
< |
derivative is desired, a fifth (or higher) order polynomial can be |
| 479 |
< |
used.\cite{Andrea83} |
| 464 |
> |
considering these molecules as neutral groups, the relative |
| 465 |
> |
orientations of the molecules control the strength of the interactions |
| 466 |
> |
at the cutoff radius. Consequently, as these molecular particles move |
| 467 |
> |
through $R_\textrm{c}$, the energy will drift upward due to the |
| 468 |
> |
anisotropy of the net molecular dipole interactions.\cite{Rahman71} To |
| 469 |
> |
maintain good energy conservation, both the potential and derivative |
| 470 |
> |
need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79} |
| 471 |
> |
This is accomplished using a standard switching function. If a smooth |
| 472 |
> |
second derivative is desired, a fifth (or higher) order polynomial can |
| 473 |
> |
be used.\cite{Andrea83} |
| 474 |
|
|
| 475 |
|
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
| 476 |
< |
and to incorporate their effect, a method like Reaction Field ({\sc |
| 477 |
< |
rf}) can be used. The original theory for {\sc rf} was originally |
| 478 |
< |
developed by Onsager,\cite{Onsager36} and it was applied in |
| 479 |
< |
simulations for the study of water by Barker and Watts.\cite{Barker73} |
| 480 |
< |
In application, it is simply an extension of the group-based cutoff |
| 481 |
< |
method where the net dipole within the cutoff sphere polarizes an |
| 482 |
< |
external dielectric, which reacts back on the central dipole. The |
| 483 |
< |
same switching function considerations for group-based cutoffs need to |
| 484 |
< |
made for {\sc rf}, with the additional pre-specification of a |
| 485 |
< |
dielectric constant. |
| 476 |
> |
and to incorporate the effects of the surroundings, a method like |
| 477 |
> |
Reaction Field ({\sc rf}) can be used. The original theory for {\sc |
| 478 |
> |
rf} was originally developed by Onsager,\cite{Onsager36} and it was |
| 479 |
> |
applied in simulations for the study of water by Barker and |
| 480 |
> |
Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply |
| 481 |
> |
an extension of the group-based cutoff method where the net dipole |
| 482 |
> |
within the cutoff sphere polarizes an external dielectric, which |
| 483 |
> |
reacts back on the central dipole. The same switching function |
| 484 |
> |
considerations for group-based cutoffs need to made for {\sc rf}, with |
| 485 |
> |
the additional pre-specification of a dielectric constant. |
| 486 |
|
|
| 487 |
|
\section{Methods} |
| 488 |
|
|
| 492 |
|
techniques utilize pairwise summations of interactions between |
| 493 |
|
particle sites, but they use these summations in different ways. |
| 494 |
|
|
| 495 |
< |
In MC, the potential energy difference between two subsequent |
| 496 |
< |
configurations dictates the progression of MC sampling. Going back to |
| 497 |
< |
the origins of this method, the acceptance criterion for the canonical |
| 498 |
< |
ensemble laid out by Metropolis \textit{et al.} states that a |
| 499 |
< |
subsequent configuration is accepted if $\Delta E < 0$ or if $\xi < |
| 500 |
< |
\exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and |
| 501 |
< |
1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an |
| 502 |
< |
alternate method for handling the long-range electrostatics will |
| 503 |
< |
ensure proper sampling from the ensemble. |
| 495 |
> |
In MC, the potential energy difference between configurations dictates |
| 496 |
> |
the progression of MC sampling. Going back to the origins of this |
| 497 |
> |
method, the acceptance criterion for the canonical ensemble laid out |
| 498 |
> |
by Metropolis \textit{et al.} states that a subsequent configuration |
| 499 |
> |
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
| 500 |
> |
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
| 501 |
> |
Maintaining the correct $\Delta E$ when using an alternate method for |
| 502 |
> |
handling the long-range electrostatics will ensure proper sampling |
| 503 |
> |
from the ensemble. |
| 504 |
|
|
| 505 |
|
In MD, the derivative of the potential governs how the system will |
| 506 |
|
progress in time. Consequently, the force and torque vectors on each |
| 513 |
|
vectors will diverge from each other more rapidly. |
| 514 |
|
|
| 515 |
|
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 516 |
+ |
|
| 517 |
|
The pairwise summation techniques (outlined in section |
| 518 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 519 |
|
studying the energy differences between conformations. We took the |
| 520 |
|
SPME-computed energy difference between two conformations to be the |
| 521 |
|
correct behavior. An ideal performance by an alternative method would |
| 522 |
< |
reproduce these energy differences exactly. Since none of the methods |
| 523 |
< |
provide exact energy differences, we used linear least squares |
| 524 |
< |
regressions of the $\Delta E$ values between configurations using SPME |
| 525 |
< |
against $\Delta E$ values using tested methods provides a quantitative |
| 526 |
< |
comparison of this agreement. Unitary results for both the |
| 527 |
< |
correlation and correlation coefficient for these regressions indicate |
| 528 |
< |
equivalent energetic results between the method under consideration |
| 529 |
< |
and electrostatics handled using SPME. Sample correlation plots for |
| 530 |
< |
two alternate methods are shown in Fig. \ref{fig:linearFit}. |
| 522 |
> |
reproduce these energy differences exactly (even if the absolute |
| 523 |
> |
energies calculated by the methods are different). Since none of the |
| 524 |
> |
methods provide exact energy differences, we used linear least squares |
| 525 |
> |
regressions of energy gap data to evaluate how closely the methods |
| 526 |
> |
mimicked the Ewald energy gaps. Unitary results for both the |
| 527 |
> |
correlation (slope) and correlation coefficient for these regressions |
| 528 |
> |
indicate perfect agreement between the alternative method and SPME. |
| 529 |
> |
Sample correlation plots for two alternate methods are shown in |
| 530 |
> |
Fig. \ref{fig:linearFit}. |
| 531 |
|
|
| 532 |
|
\begin{figure} |
| 533 |
|
\centering |
| 534 |
|
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
| 535 |
< |
\caption{Example least squares regressions of the configuration energy differences for SPC/E water systems. The upper plot shows a data set with a poor correlation coefficient ($R^2$), while the lower plot shows a data set with a good correlation coefficient.} |
| 536 |
< |
\label{fig:linearFit} |
| 535 |
> |
\caption{Example least squares regressions of the configuration energy |
| 536 |
> |
differences for SPC/E water systems. The upper plot shows a data set |
| 537 |
> |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 538 |
> |
shows a data set with a good correlation coefficient.} |
| 539 |
> |
\label{fig:linearFit} |
| 540 |
|
\end{figure} |
| 541 |
|
|
| 542 |
|
Each system type (detailed in section \ref{sec:RepSims}) was |
| 543 |
|
represented using 500 independent configurations. Additionally, we |
| 544 |
< |
used seven different system types, so each of the alternate |
| 544 |
> |
used seven different system types, so each of the alternative |
| 545 |
|
(non-Ewald) electrostatic summation methods was evaluated using |
| 546 |
|
873,250 configurational energy differences. |
| 547 |
|
|
| 571 |
|
investigated through measurement of the angle ($\theta$) formed |
| 572 |
|
between those computed from the particular method and those from SPME, |
| 573 |
|
\begin{equation} |
| 574 |
< |
\theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right), |
| 574 |
> |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 575 |
|
\end{equation} |
| 576 |
< |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 577 |
< |
force vector computed using method $M$. |
| 576 |
> |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 577 |
> |
vector computed using method M. |
| 578 |
|
|
| 579 |
|
Each of these $\theta$ values was accumulated in a distribution |
| 580 |
< |
function, weighted by the area on the unit sphere. Non-linear |
| 580 |
> |
function and weighted by the area on the unit sphere. Non-linear |
| 581 |
|
Gaussian fits were used to measure the width of the resulting |
| 582 |
|
distributions. |
| 583 |
|
|
| 584 |
|
\begin{figure} |
| 585 |
|
\centering |
| 586 |
|
\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 587 |
< |
\caption{Sample fit of the angular distribution of the force vectors over all of the studied systems. Gaussian fits were used to obtain values for the variance in force and torque vectors used in the following figure.} |
| 587 |
> |
\caption{Sample fit of the angular distribution of the force vectors |
| 588 |
> |
accumulated using all of the studied systems. Gaussian fits were used |
| 589 |
> |
to obtain values for the variance in force and torque vectors.} |
| 590 |
|
\label{fig:gaussian} |
| 591 |
|
\end{figure} |
| 592 |
|
|
| 603 |
|
when using the reference method (SPME). |
| 604 |
|
|
| 605 |
|
\subsection{Short-time Dynamics} |
| 606 |
< |
Evaluation of the short-time dynamics of charged systems was performed |
| 607 |
< |
by considering the 1000 K NaCl crystal system while using a subset of the |
| 608 |
< |
best performing pairwise methods. The NaCl crystal was chosen to |
| 609 |
< |
avoid possible complications involving the propagation techniques of |
| 610 |
< |
orientational motion in molecular systems. All systems were started |
| 611 |
< |
with the same initial positions and velocities. Simulations were |
| 612 |
< |
performed under the microcanonical ensemble, and velocity |
| 606 |
> |
|
| 607 |
> |
The effects of the alternative electrostatic summation methods on the |
| 608 |
> |
short-time dynamics of charged systems were evaluated by considering a |
| 609 |
> |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 610 |
> |
performing pairwise methods was used in this comparison. The NaCl |
| 611 |
> |
crystal was chosen to avoid possible complications from the treatment |
| 612 |
> |
of orientational motion in molecular systems. All systems were |
| 613 |
> |
started with the same initial positions and velocities. Simulations |
| 614 |
> |
were performed under the microcanonical ensemble, and velocity |
| 615 |
|
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 616 |
|
of the trajectories, |
| 617 |
|
\begin{equation} |
| 625 |
|
functions was used for comparisons. |
| 626 |
|
|
| 627 |
|
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 628 |
< |
Evaluation of the long-time dynamics of charged systems was performed |
| 629 |
< |
by considering the NaCl crystal system, again while using a subset of |
| 630 |
< |
the best performing pairwise methods. To enhance the atomic motion, |
| 631 |
< |
these crystals were equilibrated at 1000 K, near the experimental |
| 632 |
< |
$T_m$ for NaCl. Simulations were performed under the microcanonical |
| 633 |
< |
ensemble, and velocity information was saved every 5 fs over 100 ps |
| 634 |
< |
trajectories. The power spectrum ($I(\omega)$) was obtained via |
| 633 |
< |
Fourier transform of the velocity autocorrelation function |
| 634 |
< |
\begin{equation} |
| 635 |
< |
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 628 |
> |
|
| 629 |
> |
The effects of the same subset of alternative electrostatic methods on |
| 630 |
> |
the {\it long-time} dynamics of charged systems were evaluated using |
| 631 |
> |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 632 |
> |
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 633 |
> |
autocorrelation function, \begin{equation} I(\omega) = |
| 634 |
> |
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 635 |
|
\label{eq:powerSpec} |
| 636 |
|
\end{equation} |
| 637 |
|
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 638 |
|
NaCl crystal is composed of two different atom types, the average of |
| 639 |
< |
the two resulting power spectra was used for comparisons. |
| 639 |
> |
the two resulting power spectra was used for comparisons. Simulations |
| 640 |
> |
were performed under the microcanonical ensemble, and velocity |
| 641 |
> |
information was saved every 5 fs over 100 ps trajectories. |
| 642 |
|
|
| 643 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 644 |
< |
A variety of common and representative simulations were analyzed to |
| 645 |
< |
determine the relative effectiveness of the pairwise summation |
| 646 |
< |
techniques in reproducing the energetics and dynamics exhibited by |
| 647 |
< |
SPME. The studied systems were as follows: |
| 644 |
> |
A variety of representative simulations were analyzed to determine the |
| 645 |
> |
relative effectiveness of the pairwise summation techniques in |
| 646 |
> |
reproducing the energetics and dynamics exhibited by SPME. We wanted |
| 647 |
> |
to span the space of modern simulations (i.e. from liquids of neutral |
| 648 |
> |
molecules to ionic crystals), so the systems studied were: |
| 649 |
|
\begin{enumerate} |
| 650 |
< |
\item Liquid Water |
| 651 |
< |
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 652 |
< |
\item NaCl Crystal |
| 653 |
< |
\item NaCl Melt |
| 654 |
< |
\item Low Ionic Strength Solution of NaCl in Water |
| 655 |
< |
\item High Ionic Strength Solution of NaCl in Water |
| 656 |
< |
\item 6 \AA\ Radius Sphere of Argon in Water |
| 650 |
> |
\item liquid water (SPC/E),\cite{Berendsen87} |
| 651 |
> |
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 652 |
> |
\item NaCl crystals, |
| 653 |
> |
\item NaCl melts, |
| 654 |
> |
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 655 |
> |
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 656 |
> |
\item a 6 \AA\ radius sphere of Argon in water. |
| 657 |
|
\end{enumerate} |
| 658 |
|
By utilizing the pairwise techniques (outlined in section |
| 659 |
|
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 660 |
< |
charged particles, and mixtures of the two, we can comment on possible |
| 661 |
< |
system dependence and/or universal applicability of the techniques. |
| 660 |
> |
charged particles, and mixtures of the two, we hope to discern under |
| 661 |
> |
which conditions it will be possible to use one of the alternative |
| 662 |
> |
summation methodologies instead of the Ewald sum. |
| 663 |
|
|
| 664 |
< |
Generation of the system configurations was dependent on the system |
| 665 |
< |
type. For the solid and liquid water configurations, configuration |
| 666 |
< |
snapshots were taken at regular intervals from higher temperature 1000 |
| 667 |
< |
SPC/E water molecule trajectories and each equilibrated individually. |
| 668 |
< |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 669 |
< |
ions and were selected and equilibrated in the same fashion as the |
| 670 |
< |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 671 |
< |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 672 |
< |
respectively. Ion and water positions were then randomly swapped, and |
| 673 |
< |
the resulting configurations were again equilibrated individually. |
| 674 |
< |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 675 |
< |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 676 |
< |
water configurations were converted to argon |
| 664 |
> |
For the solid and liquid water configurations, configurations were |
| 665 |
> |
taken at regular intervals from high temperature trajectories of 1000 |
| 666 |
> |
SPC/E water molecules. Each configuration was equilibrated |
| 667 |
> |
independently at a lower temperature (300~K for the liquid, 200~K for |
| 668 |
> |
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 669 |
> |
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 670 |
> |
these systems were selected and equilibrated in the same manner as the |
| 671 |
> |
water systems. The equilibrated temperatures were 1000~K for the NaCl |
| 672 |
> |
crystal and 7000~K for the liquid. The ionic solutions were made by |
| 673 |
> |
solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water |
| 674 |
> |
molecules. Ion and water positions were then randomly swapped, and |
| 675 |
> |
the resulting configurations were again equilibrated individually. |
| 676 |
> |
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 677 |
> |
of all the SPC/E waters within 6 \AA\ of the center of the |
| 678 |
> |
equilibrated water configurations were converted to argon |
| 679 |
|
(Fig. \ref{fig:argonSlice}). |
| 680 |
|
|
| 681 |
+ |
These procedures guaranteed us a set of representative configurations |
| 682 |
+ |
from chemically-relevant systems sampled from an appropriate |
| 683 |
+ |
ensemble. Force field parameters for the ions and Argon were taken |
| 684 |
+ |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 685 |
+ |
|
| 686 |
|
\begin{figure} |
| 687 |
|
\centering |
| 688 |
|
\includegraphics[width = \linewidth]{./slice.pdf} |
| 689 |
< |
\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.} |
| 689 |
> |
\caption{A slice from the center of a water box used in a charge void |
| 690 |
> |
simulation. The darkened region represents the boundary sphere within |
| 691 |
> |
which the water molecules were converted to argon atoms.} |
| 692 |
|
\label{fig:argonSlice} |
| 693 |
|
\end{figure} |
| 694 |
|
|
| 695 |
< |
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 696 |
< |
Electrostatic summation method comparisons were performed using SPME, |
| 697 |
< |
the {\sc sp} and {\sc sf} methods - both with damping |
| 698 |
< |
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, |
| 699 |
< |
moderate, and strong damping respectively), reaction field with an |
| 700 |
< |
infinite dielectric constant, and an unmodified cutoff. Group-based |
| 701 |
< |
cutoffs with a fifth-order polynomial switching function were |
| 702 |
< |
necessary for the reaction field simulations and were utilized in the |
| 703 |
< |
SP, SF, and pure cutoff methods for comparison to the standard lack of |
| 704 |
< |
group-based cutoffs with a hard truncation. The SPME calculations |
| 705 |
< |
were performed using the TINKER implementation of SPME,\cite{Ponder87} |
| 706 |
< |
while all other method calculations were performed using the OOPSE |
| 707 |
< |
molecular mechanics package.\cite{Meineke05} |
| 695 |
> |
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 696 |
> |
We compared the following alternative summation methods with results |
| 697 |
> |
from the reference method (SPME): |
| 698 |
> |
\begin{itemize} |
| 699 |
> |
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 700 |
> |
and 0.3 \AA$^{-1}$, |
| 701 |
> |
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 702 |
> |
and 0.3 \AA$^{-1}$, |
| 703 |
> |
\item reaction field with an infinite dielectric constant, and |
| 704 |
> |
\item an unmodified cutoff. |
| 705 |
> |
\end{itemize} |
| 706 |
> |
Group-based cutoffs with a fifth-order polynomial switching function |
| 707 |
> |
were utilized for the reaction field simulations. Additionally, we |
| 708 |
> |
investigated the use of these cutoffs with the SP, SF, and pure |
| 709 |
> |
cutoff. The SPME electrostatics were performed using the TINKER |
| 710 |
> |
implementation of SPME,\cite{Ponder87} while all other method |
| 711 |
> |
calculations were performed using the OOPSE molecular mechanics |
| 712 |
> |
package.\cite{Meineke05} All other portions of the energy calculation |
| 713 |
> |
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 714 |
> |
manner across all systems and configurations. |
| 715 |
|
|
| 716 |
< |
These methods were additionally evaluated with three different cutoff |
| 717 |
< |
radii (9, 12, and 15 \AA) to investigate possible cutoff radius |
| 718 |
< |
dependence. It should be noted that the damping parameter chosen in |
| 719 |
< |
SPME, or so called ``Ewald Coefficient", has a significant effect on |
| 720 |
< |
the energies and forces calculated. Typical molecular mechanics |
| 721 |
< |
packages default this to a value dependent on the cutoff radius and a |
| 722 |
< |
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 723 |
< |
tolerances are typically associated with increased accuracy, but this |
| 724 |
< |
usually means more time spent calculating the reciprocal-space portion |
| 725 |
< |
of the summation.\cite{Perram88,Essmann95} The default TINKER |
| 726 |
< |
tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 727 |
< |
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 728 |
< |
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 716 |
> |
The althernative methods were also evaluated with three different |
| 717 |
> |
cutoff radii (9, 12, and 15 \AA). It should be noted that the damping |
| 718 |
> |
parameter chosen in SPME, or so called ``Ewald Coefficient'', has a |
| 719 |
> |
significant effect on the energies and forces calculated. Typical |
| 720 |
> |
molecular mechanics packages set this to a value dependent on the |
| 721 |
> |
cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$ |
| 722 |
> |
kcal/mol). Smaller tolerances are typically associated with increased |
| 723 |
> |
accuracy at the expense of increased time spent calculating the |
| 724 |
> |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 725 |
> |
The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 726 |
> |
in all SPME calculations, resulting in Ewald Coefficients of 0.4200, |
| 727 |
> |
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 728 |
> |
respectively. |
| 729 |
|
|
| 730 |
|
\section{Results and Discussion} |
| 731 |
|
|
| 739 |
|
\begin{figure} |
| 740 |
|
\centering |
| 741 |
|
\includegraphics[width=5.5in]{./delEplot.pdf} |
| 742 |
< |
\caption{Statistical analysis of the quality of configurational energy differences for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate $\Delta E$ values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 742 |
> |
\caption{Statistical analysis of the quality of configurational energy |
| 743 |
> |
differences for a given electrostatic method compared with the |
| 744 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 745 |
> |
indicate $\Delta E$ values indistinguishable from those obtained using |
| 746 |
> |
SPME. Different values of the cutoff radius are indicated with |
| 747 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 748 |
> |
inverted triangles).} |
| 749 |
|
\label{fig:delE} |
| 750 |
|
\end{figure} |
| 751 |
|
|
| 752 |
< |
In this figure, it is apparent that it is unreasonable to expect |
| 753 |
< |
realistic results using an unmodified cutoff. This is not all that |
| 754 |
< |
surprising since this results in large energy fluctuations as atoms |
| 755 |
< |
move in and out of the cutoff radius. These fluctuations can be |
| 756 |
< |
alleviated to some degree by using group based cutoffs with a |
| 757 |
< |
switching function.\cite{Steinbach94} The Group Switch Cutoff row |
| 733 |
< |
doesn't show a significant improvement in this plot because the salt |
| 734 |
< |
and salt solution systems contain non-neutral groups, see the |
| 735 |
< |
accompanying supporting information for a comparison where all groups |
| 736 |
< |
are neutral. |
| 752 |
> |
The most striking feature of this plot is how well the Shifted Force |
| 753 |
> |
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 754 |
> |
differences. For the undamped {\sc sf} method, and the |
| 755 |
> |
moderately-damped {\sc sp} methods, the results are nearly |
| 756 |
> |
indistinguishable from the Ewald results. The other common methods do |
| 757 |
> |
significantly less well. |
| 758 |
|
|
| 759 |
< |
Correcting the resulting charged cutoff sphere is one of the purposes |
| 760 |
< |
of the damped Coulomb summation proposed by Wolf \textit{et |
| 761 |
< |
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 762 |
< |
seen in the Shifted-Potental rows. While the undamped case of this |
| 763 |
< |
method is a significant improvement over the pure cutoff, it still |
| 764 |
< |
doesn't correlate that well with SPME. Inclusion of potential damping |
| 765 |
< |
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 759 |
> |
The unmodified cutoff method is essentially unusable. This is not |
| 760 |
> |
surprising since hard cutoffs give large energy fluctuations as atoms |
| 761 |
> |
or molecules move in and out of the cutoff |
| 762 |
> |
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
| 763 |
> |
some degree by using group based cutoffs with a switching |
| 764 |
> |
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 765 |
> |
significant improvement using the group-switched cutoff because the |
| 766 |
> |
salt and salt solution systems contain non-neutral groups. Interested |
| 767 |
> |
readers can consult the accompanying supporting information for a |
| 768 |
> |
comparison where all groups are neutral. |
| 769 |
> |
|
| 770 |
> |
For the {\sc sp} method, inclusion of potential damping improves the |
| 771 |
> |
agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 772 |
|
an excellent correlation and quality of fit with the SPME results, |
| 773 |
< |
particularly with a cutoff radius greater than 12 \AA . Use of a |
| 774 |
< |
larger damping parameter is more helpful for the shortest cutoff |
| 775 |
< |
shown, but it has a detrimental effect on simulations with larger |
| 776 |
< |
cutoffs. In the {\sc sf} sets, increasing damping results in |
| 777 |
< |
progressively poorer correlation. Overall, the undamped case is the |
| 778 |
< |
best performing set, as the correlation and quality of fits are |
| 779 |
< |
consistently superior regardless of the cutoff distance. This result |
| 780 |
< |
is beneficial in that the undamped case is less computationally |
| 781 |
< |
prohibitive do to the lack of complimentary error function calculation |
| 782 |
< |
when performing the electrostatic pair interaction. The reaction |
| 783 |
< |
field results illustrates some of that method's limitations, primarily |
| 784 |
< |
that it was developed for use in homogenous systems; although it does |
| 785 |
< |
provide results that are an improvement over those from an unmodified |
| 786 |
< |
cutoff. |
| 773 |
> |
particularly with a cutoff radius greater than 12 |
| 774 |
> |
\AA . Use of a larger damping parameter is more helpful for the |
| 775 |
> |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 776 |
> |
with larger cutoffs. |
| 777 |
> |
|
| 778 |
> |
In the {\sc sf} sets, increasing damping results in progressively |
| 779 |
> |
worse correlation with Ewald. Overall, the undamped case is the best |
| 780 |
> |
performing set, as the correlation and quality of fits are |
| 781 |
> |
consistently superior regardless of the cutoff distance. The undamped |
| 782 |
> |
case is also less computationally demanding (because no evaluation of |
| 783 |
> |
the complementary error function is required). |
| 784 |
> |
|
| 785 |
> |
The reaction field results illustrates some of that method's |
| 786 |
> |
limitations, primarily that it was developed for use in homogenous |
| 787 |
> |
systems; although it does provide results that are an improvement over |
| 788 |
> |
those from an unmodified cutoff. |
| 789 |
|
|
| 790 |
|
\subsection{Magnitudes of the Force and Torque Vectors} |
| 791 |
|
|
| 980 |
|
increased, these peaks are smoothed out, and approach the SPME |
| 981 |
|
curve. The damping acts as a distance dependent Gaussian screening of |
| 982 |
|
the point charges for the pairwise summation methods; thus, the |
| 983 |
< |
collisions are more elastic in the undamped {\sc sf} potental, and the |
| 983 |
> |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 984 |
|
stiffness of the potential is diminished as the electrostatic |
| 985 |
|
interactions are softened by the damping function. With $\alpha$ |
| 986 |
|
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
