| 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 |
|
|
| 655 |
|
Generation of the system configurations was dependent on the system |
| 656 |
|
type. For the solid and liquid water configurations, configuration |
| 657 |
|
snapshots were taken at regular intervals from higher temperature 1000 |
| 658 |
< |
SPC/E water molecule trajectories and each equilibrated individually. |
| 659 |
< |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 660 |
< |
ions and were selected and equilibrated in the same fashion as the |
| 661 |
< |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 662 |
< |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 663 |
< |
respectively. Ion and water positions were then randomly swapped, and |
| 664 |
< |
the resulting configurations were again equilibrated individually. |
| 665 |
< |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 666 |
< |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 667 |
< |
water configurations were converted to argon |
| 668 |
< |
(Fig. \ref{fig:argonSlice}). |
| 658 |
> |
SPC/E water molecule trajectories and each equilibrated |
| 659 |
> |
individually.\cite{Berendsen87} The solid and liquid NaCl systems |
| 660 |
> |
consisted of 500 Na+ and 500 Cl- ions and were selected and |
| 661 |
> |
equilibrated in the same fashion as the water systems. For the low |
| 662 |
> |
and high ionic strength NaCl solutions, 4 and 40 ions were first |
| 663 |
> |
solvated in a 1000 water molecule boxes respectively. Ion and water |
| 664 |
> |
positions were then randomly swapped, and the resulting configurations |
| 665 |
> |
were again equilibrated individually. Finally, for the Argon/Water |
| 666 |
> |
"charge void" systems, the identities of all the SPC/E waters within 6 |
| 667 |
> |
\AA\ of the center of the equilibrated water configurations were |
| 668 |
> |
converted to argon (Fig. \ref{fig:argonSlice}). |
| 669 |
|
|
| 670 |
|
\begin{figure} |
| 671 |
|
\centering |
| 720 |
|
|
| 721 |
|
In this figure, it is apparent that it is unreasonable to expect |
| 722 |
|
realistic results using an unmodified cutoff. This is not all that |
| 723 |
< |
surprising since this results in large energy fluctuations as atoms |
| 724 |
< |
move in and out of the cutoff radius. These fluctuations can be |
| 725 |
< |
alleviated to some degree by using group based cutoffs with a |
| 726 |
< |
switching function.\cite{Steinbach94} The Group Switch Cutoff row |
| 727 |
< |
doesn't show a significant improvement in this plot because the salt |
| 728 |
< |
and salt solution systems contain non-neutral groups, see the |
| 723 |
> |
surprising since this results in large energy fluctuations as atoms or |
| 724 |
> |
molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79} |
| 725 |
> |
These fluctuations can be alleviated to some degree by using group |
| 726 |
> |
based cutoffs with a switching |
| 727 |
> |
function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff |
| 728 |
> |
row doesn't show a significant improvement in this plot because the |
| 729 |
> |
salt and salt solution systems contain non-neutral groups, see the |
| 730 |
|
accompanying supporting information for a comparison where all groups |
| 731 |
|
are neutral. |
| 732 |
|
|
| 733 |
|
Correcting the resulting charged cutoff sphere is one of the purposes |
| 734 |
|
of the damped Coulomb summation proposed by Wolf \textit{et |
| 735 |
|
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 736 |
< |
seen in the Shifted-Potental rows. While the undamped case of this |
| 736 |
> |
seen in the {\sc sp} rows. While the undamped case of this |
| 737 |
|
method is a significant improvement over the pure cutoff, it still |
| 738 |
|
doesn't correlate that well with SPME. Inclusion of potential damping |
| 739 |
|
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 946 |
|
increased, these peaks are smoothed out, and approach the SPME |
| 947 |
|
curve. The damping acts as a distance dependent Gaussian screening of |
| 948 |
|
the point charges for the pairwise summation methods; thus, the |
| 949 |
< |
collisions are more elastic in the undamped {\sc sf} potental, and the |
| 949 |
> |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 950 |
|
stiffness of the potential is diminished as the electrostatic |
| 951 |
|
interactions are softened by the damping function. With $\alpha$ |
| 952 |
|
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
