| 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 an |
| 81 |
< |
impractical task to perform these calculations. |
| 80 |
> |
increases proportionally with the cutoff sphere, it quickly becomes |
| 81 |
> |
very time-consuming to perform these calculations. |
| 82 |
|
|
| 83 |
+ |
There have been many efforts to address this issue of both proper and |
| 84 |
+ |
practical handling of electrostatic interactions, and these have |
| 85 |
+ |
resulted in the availability of a variety of |
| 86 |
+ |
techniques.\cite{Roux99,Sagui99,Tobias01} These are typically |
| 87 |
+ |
classified as implicit methods (i.e., continuum dielectrics, static |
| 88 |
+ |
dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e., |
| 89 |
+ |
Ewald summations, interaction shifting or |
| 90 |
+ |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 91 |
+ |
reaction field type methods, fast multipole |
| 92 |
+ |
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 93 |
+ |
often preferred because they incorporate dynamic solvent molecules in |
| 94 |
+ |
the system of interest, but these methods are sometimes difficult to |
| 95 |
+ |
utilize because of their high computational cost.\cite{Roux99} In |
| 96 |
+ |
addition to this cost, there has been some question of the inherent |
| 97 |
+ |
periodicity of the explicit Ewald summation artificially influencing |
| 98 |
+ |
systems dynamics.\cite{Tobias01} |
| 99 |
+ |
|
| 100 |
+ |
In this paper, we focus on the common mixed and explicit methods of |
| 101 |
+ |
reaction filed and smooth particle mesh |
| 102 |
+ |
Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods |
| 103 |
+ |
devised by Wolf {\it et al.} which we further extend.\cite{Wolf99} |
| 104 |
+ |
These new methods for handling electrostatics are quite |
| 105 |
+ |
computationally efficient, since they involve only a simple |
| 106 |
+ |
modification to the direct pairwise sum, and they lack the added |
| 107 |
+ |
periodicity of the Ewald sum. Below, these methods are evaluated using |
| 108 |
+ |
a variety of model systems and comparison methodologies to establish |
| 109 |
+ |
their usability in molecular simulations. |
| 110 |
+ |
|
| 111 |
|
\subsection{The Ewald Sum} |
| 112 |
< |
The complete accumulation electrostatic interactions in a system with periodic boundary conditions (PBC) requires the consideration of the effect of all charges within a simulation box, as well as those in the periodic replicas, |
| 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, |
| 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} |
| 119 |
|
\end{equation} |
| 120 |
< |
where the sum over $\mathbf{n}$ is a sum over all periodic box replicas |
| 121 |
< |
with integer coordinates $\mathbf{n} = (l,m,n)$, and the prime indicates |
| 122 |
< |
$i = j$ are neglected for $\mathbf{n} = 0$.\cite{deLeeuw80} Within the |
| 123 |
< |
sum, $N$ is the number of electrostatic particles, $\mathbf{r}_{ij}$ is |
| 124 |
< |
$\mathbf{r}_j - \mathbf{r}_i$, $L$ is the cell length, $\bm{\Omega}_{i,j}$ are |
| 125 |
< |
the Euler angles for $i$ and $j$, and $\phi$ is Poisson's equation |
| 126 |
< |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge |
| 127 |
< |
interactions). In the case of monopole electrostatics, |
| 128 |
< |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is discontiuous |
| 129 |
< |
for non-neutral systems. |
| 120 |
> |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
| 121 |
> |
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
| 122 |
> |
prime indicates $i = j$ are neglected for $\mathbf{n} = |
| 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 discontinuous for non-neutral systems. |
| 130 |
|
|
| 131 |
|
This electrostatic summation problem was originally studied by Ewald |
| 132 |
|
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 176 |
|
direct and reciprocal-space portions of the summation. The choice of |
| 177 |
|
the magnitude of this value allows one to select whether the |
| 178 |
|
real-space or reciprocal space portion of the summation is an |
| 179 |
< |
$\mathscr{O}(N^2)$ calcualtion (with the other being |
| 179 |
> |
$\mathscr{O}(N^2)$ calculation (with the other being |
| 180 |
|
$\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$ |
| 181 |
|
and thoughtful algorithm development, this cost can be brought down to |
| 182 |
|
$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to |
| 212 |
|
artificially stabilized by the periodic replicas introduced by the |
| 213 |
|
Ewald summation.\cite{Weber00} Thus, care ought to be taken when |
| 214 |
|
considering the use of the Ewald summation where the intrinsic |
| 215 |
< |
perodicity may negatively affect the system dynamics. |
| 215 |
> |
periodicity may negatively affect the system dynamics. |
| 216 |
|
|
| 217 |
|
|
| 218 |
|
\subsection{The Wolf and Zahn Methods} |
| 230 |
|
and a distance-dependent damping function (identical to that seen in |
| 231 |
|
the real-space portion of the Ewald sum) to aid convergence |
| 232 |
|
\begin{equation} |
| 233 |
< |
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\}. |
| 233 |
> |
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\}. |
| 234 |
|
\label{eq:WolfPot} |
| 235 |
|
\end{equation} |
| 236 |
|
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 573 |
|
investigated through measurement of the angle ($\theta$) formed |
| 574 |
|
between those computed from the particular method and those from SPME, |
| 575 |
|
\begin{equation} |
| 576 |
< |
\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, |
| 576 |
> |
\theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right), |
| 577 |
|
\end{equation} |
| 578 |
|
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 579 |
|
force vector computed using method $M$. |
| 661 |
|
Generation of the system configurations was dependent on the system |
| 662 |
|
type. For the solid and liquid water configurations, configuration |
| 663 |
|
snapshots were taken at regular intervals from higher temperature 1000 |
| 664 |
< |
SPC/E water molecule trajectories and each equilibrated individually. |
| 665 |
< |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 666 |
< |
ions and were selected and equilibrated in the same fashion as the |
| 667 |
< |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 668 |
< |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 669 |
< |
respectively. Ion and water positions were then randomly swapped, and |
| 670 |
< |
the resulting configurations were again equilibrated individually. |
| 671 |
< |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 672 |
< |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 673 |
< |
water configurations were converted to argon |
| 674 |
< |
(Fig. \ref{fig:argonSlice}). |
| 664 |
> |
SPC/E water molecule trajectories and each equilibrated |
| 665 |
> |
individually.\cite{Berendsen87} The solid and liquid NaCl systems |
| 666 |
> |
consisted of 500 Na+ and 500 Cl- ions and were selected and |
| 667 |
> |
equilibrated in the same fashion as the water systems. For the low |
| 668 |
> |
and high ionic strength NaCl solutions, 4 and 40 ions were first |
| 669 |
> |
solvated in a 1000 water molecule boxes respectively. Ion and water |
| 670 |
> |
positions were then randomly swapped, and the resulting configurations |
| 671 |
> |
were again equilibrated individually. Finally, for the Argon/Water |
| 672 |
> |
"charge void" systems, the identities of all the SPC/E waters within 6 |
| 673 |
> |
\AA\ of the center of the equilibrated water configurations were |
| 674 |
> |
converted to argon (Fig. \ref{fig:argonSlice}). |
| 675 |
|
|
| 676 |
|
\begin{figure} |
| 677 |
|
\centering |
| 726 |
|
|
| 727 |
|
In this figure, it is apparent that it is unreasonable to expect |
| 728 |
|
realistic results using an unmodified cutoff. This is not all that |
| 729 |
< |
surprising since this results in large energy fluctuations as atoms |
| 730 |
< |
move in and out of the cutoff radius. These fluctuations can be |
| 731 |
< |
alleviated to some degree by using group based cutoffs with a |
| 732 |
< |
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 |
| 729 |
> |
surprising since this results in large energy fluctuations as atoms or |
| 730 |
> |
molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79} |
| 731 |
> |
These fluctuations can be alleviated to some degree by using group |
| 732 |
> |
based cutoffs with a switching |
| 733 |
> |
function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff |
| 734 |
> |
row doesn't show a significant improvement in this plot because the |
| 735 |
> |
salt and salt solution systems contain non-neutral groups, see the |
| 736 |
|
accompanying supporting information for a comparison where all groups |
| 737 |
|
are neutral. |
| 738 |
|
|
| 739 |
|
Correcting the resulting charged cutoff sphere is one of the purposes |
| 740 |
|
of the damped Coulomb summation proposed by Wolf \textit{et |
| 741 |
|
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 742 |
< |
seen in the Shifted-Potental rows. While the undamped case of this |
| 742 |
> |
seen in the {\sc sp} rows. While the undamped case of this |
| 743 |
|
method is a significant improvement over the pure cutoff, it still |
| 744 |
|
doesn't correlate that well with SPME. Inclusion of potential damping |
| 745 |
|
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 952 |
|
increased, these peaks are smoothed out, and approach the SPME |
| 953 |
|
curve. The damping acts as a distance dependent Gaussian screening of |
| 954 |
|
the point charges for the pairwise summation methods; thus, the |
| 955 |
< |
collisions are more elastic in the undamped {\sc sf} potental, and the |
| 955 |
> |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 956 |
|
stiffness of the potential is diminished as the electrostatic |
| 957 |
|
interactions are softened by the damping function. With $\alpha$ |
| 958 |
|
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
