| 57 |
|
\begin{figure} |
| 58 |
|
\centering |
| 59 |
|
\includegraphics[width = 3.25in]{./ewaldProgression.pdf} |
| 60 |
< |
\caption{How the application of the Ewald summation has changed with the increase in computer power. Initially, only small numbers of particles could be studied, and the Ewald sum acted to replicate the unit cell charge distribution out to infinity. Now, much larger systems of charges are investigated with fixed distance cutoffs. The calculated structure factor is used to sum out to great distance, and a surrounding dielectric term is included.} |
| 60 |
> |
\caption{How the application of the Ewald summation has changed with the increase in computer power. Initially, only small numbers of particles could be studied, and the Ewald sum acted to replicate the unit cell charge distribution out to convergence. Now, much larger systems of charges are investigated with fixed distance cutoffs. The calculated structure factor is used to sum out to great distance, and a surrounding dielectric term is included.} |
| 61 |
|
\label{fig:ewaldTime} |
| 62 |
|
\end{figure} |
| 63 |
|
|
| 69 |
|
\end{equation} |
| 70 |
|
In order to use this potential in molecular dynamics simulations, Wolf \textit{et al.} suggested taking the derivative of this potential, followed by evaluation of the limit to give the following forces, |
| 71 |
|
\begin{equation} |
| 72 |
< |
F^{\textrm{Wolf}}(r_{ij}) = q_iq_j\left[\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right)-\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right]. |
| 72 |
> |
F^{\textrm{Wolf}}(r_{ij}) = q_iq_j\left\{\left[-\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}\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\}. |
| 73 |
|
\label{eq:WolfForces} |
| 74 |
|
\end{equation} |
| 75 |
|
More recently, Zahn \textit{et al.} investigated this electrostatic summation method for use in simulations involving water.\cite{Zahn02} In their work, they point out that the method as proposed is problematic for use in Molecular Dynamics simulations, because the forces and derivative of the potential are not equivalent. This comes about from the procedure of taking the limit shown in equation \ref{eq:WolfPot} after calculating the derivatives.\cite{Wolf99} Zahn \textit{et al.} proposed a shifted force adaptation of this ``Wolf summation method" as a way to use this technique in Molecular Dynamics simulations. Taking the integral of the forces shown in equation \ref{eq:WolfForces}, they obtained a new shifted damped Coulomb potential |
| 80 |
|
They showed that this new potential does well in capturing the structural and dynamic properties present when using the Ewald sum with the models of water used in their simulations. |
| 81 |
|
|
| 82 |
|
\subsection{Simple Forms for Pairwise Electrostatics} |
| 83 |
< |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et al.} are constructed using two different (and separable) computational tricks: shifting through use of image charges and damping of the electrostatic interaction. |
| 83 |
> |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et al.} are constructed using two different (and separable) computational tricks: shifting through use of image charges and damping of the electrostatic interaction. Wolf \textit{et al.} treated the development of their summation method as a progressive application of these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded their shifted force adaptation \ref{eq:ZahnPot} on what they called "the formally incorrect prescription for the derivation of damped Coulomb pair forces".\cite{Zahn02} Below, we consider the ideas encompassing these electrostatic summation method formulations and clarify their development. |
| 84 |
|
|
| 85 |
< |
While implementing these methods for use in our own work, we discovered the potential presented in equation \ref{eq:ZahnPot} is still not entirely correct. The derivative of this equation leads to a sign error in the forces, resulting in erroneous dynamics. We can apply the standard shifted force potential, |
| 85 |
> |
Starting with the original observation that the effective range of the electrostatic interaction in condensed phases is considerably less than the $r^{-1}$ in vacuum, either the shifting or the distance-dependent damping technique could be used as a foundation for the summation method. Wolf \textit{et al.} made the additional observation that charge neutralization within the cutoff sphere plays a significant role in energy convergence; thus, shifting through the use of image charges was taken as the initial step. Using these image charges, the electrostatic summation is forced to converge at the cutoff radius. We can incorporate the methods of Wolf \textit{et al.} and Zahn \textit{et al.} by considering the standard shifted force potential |
| 86 |
|
\begin{equation} |
| 87 |
< |
V^\textrm{SF}(r_{ij}) = \begin{cases} v(r_{ij})-v_\textrm{c}-\left(\frac{\textrm{d}v(r_{ij})}{\textrm{d}r_{ij}}\right)_{r_{ij}=R_\textrm{c}}(r_{ij}-R_\textrm{c}) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 87 |
> |
V^\textrm{SF}(r_{ij}) = \begin{cases} v(r_{ij})-v_\textrm{c}-\left[\frac{\textrm{d}v(r_{ij})}{\textrm{d}r_{ij}}\right]_{r_{ij}=R_\textrm{c}}(r_{ij}-R_\textrm{c}) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 88 |
|
\end{cases}, |
| 89 |
+ |
\label{eq:shiftingForm} |
| 90 |
|
\end{equation} |
| 91 |
< |
where $v(r_{ij})$ is the unshifted form of the potential, and $v_c$ is $v(R_\textrm{c})$ and insures the potential goes to zero at the cutoff radius.\cite{Allen87} Using the simple damped Coulomb potential as the starting point, |
| 91 |
> |
where $v(r_{ij})$ is the unshifted form of the potential, and $v_c$ is $v(R_\textrm{c})$ and insures the potential goes to zero at the cutoff radius.\cite{Allen87} If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99} |
| 92 |
|
\begin{equation} |
| 93 |
< |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 94 |
< |
\label{eq:dampCoulomb} |
| 93 |
> |
V^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 94 |
> |
\end{cases}. |
| 95 |
> |
\label{eq:WolfSP} |
| 96 |
|
\end{equation} |
| 97 |
< |
the resulting shifted force potential is |
| 97 |
> |
The forces associated with this potential are obtained by taking the derivative, resulting in the following, |
| 98 |
|
\begin{equation} |
| 99 |
< |
V^\mathrm{SF}\left(r_{ij}\right)=q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}}+\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. |
| 99 |
> |
F^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 100 |
> |
\end{cases}. |
| 101 |
> |
\label{eq:FWolfSP} |
| 102 |
> |
\end{equation} |
| 103 |
> |
These forces are identical to the forces of the standard electrostatic interaction, and this was addressed by Wolf \textit{et al.} as undesirable. They pointed out that the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. Additionally, there is a discontinuity in the forces. This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$. |
| 104 |
> |
|
| 105 |
> |
If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential, |
| 106 |
> |
\begin{equation} |
| 107 |
> |
V^\textrm{SF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left[\frac{1}{R_\textrm{c}^2}\right](r_{ij}-R_\textrm{c})\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 108 |
> |
\end{cases}. |
| 109 |
|
\label{eq:SFPot} |
| 110 |
|
\end{equation} |
| 111 |
< |
Equation \ref{eq:SFPot} is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term is simply equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in equation \ref{eq:ZahnPot}, resulting in a discontinuity in the potential as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 111 |
> |
Taking the derivative of this shifted force potential gives the following forces, |
| 112 |
> |
\begin{equation} |
| 113 |
> |
F^\textrm{SF}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 114 |
> |
\end{cases}. |
| 115 |
> |
\label{eq:SFForces} |
| 116 |
> |
\end{equation} |
| 117 |
> |
Using this formulation rather than the simple shifted potential (Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the forces in addition to the potential. This form also has the benefit that the image charges have a force presence, addressing the concerns about a missing physical component. One side effect of this treatment is a slight alteration in the shape of the potential that comes about from the derivative term. Thus, a degree of clarity about the original formulation of the potential is lost in order to gain functionality in dynamics simulations. |
| 118 |
|
|
| 119 |
< |
It is important to note that shifted force techniques have a drawback in that they alter the shape of the original potential. We thereby lose a degree of clarity about the original formulation of the potential in order to gain functionality in dynamics simulations. An alternative direction would be use the derivatives of the original potential for the forces. This was addressed by Wolf \textit{et al.} as undesirable, because the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. In Monte Carlo simulations, this argument is moot, because forces are not evaluated. We decided to consider both the Shifted-Force technique described above and this Shifted-Potential technique to determine their usability in the evaluation of both energetic and dynamic results in simulations with electrostatics. |
| 119 |
> |
Wolf \textit{et al.} originally discussed the energetics of the shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that it was still insufficient for accurate determination of the energy. The energy would fluctuate around the expected value with increasing cutoff radius, but the oscillations appeared to be converging toward the correct value.\cite{Wolf99} A damping function was incorporated to accelerate convergence; and though alternative functional forms could be used,\cite{Jones56,Heyes81} the complimentary error function was chosen to draw parallels to the Ewald summation. Incorporating damping into the simple Coulomb potential, |
| 120 |
> |
\begin{equation} |
| 121 |
> |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 122 |
> |
\label{eq:dampCoulomb} |
| 123 |
> |
\end{equation} |
| 124 |
> |
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived \textit{via} equation \ref{eq:shiftingForm}, |
| 125 |
> |
\begin{equation} |
| 126 |
> |
V^{\textrm{DSP}}(r_{ij}) = \begin{cases} q_iq_j\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right] &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 127 |
> |
\end{cases}. |
| 128 |
> |
\label{eq:DSPPot} |
| 129 |
> |
\end{equation} |
| 130 |
> |
The derivative of this Shifted-Potential can be taken to obtain forces for use in MD, |
| 131 |
> |
\begin{equation} |
| 132 |
> |
F^{\textrm{DSP}}(r_{ij}) = \begin{cases} q_iq_j\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right] &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 133 |
> |
\end{cases}. |
| 134 |
> |
\label{eq:DSPForces} |
| 135 |
> |
\end{equation} |
| 136 |
> |
Again, this Shifted-Potential suffers from a discontinuity in the forces, and a lack of an image-charge component in the forces. To remedy these concerns, a Shifted-Force variant is obtained by inclusion of the derivative term in equation \ref{eq:shiftingForm} to give, |
| 137 |
> |
\begin{equation} |
| 138 |
> |
V^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}}\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 139 |
> |
\end{cases}. |
| 140 |
> |
\label{eq:DSFPot} |
| 141 |
> |
\end{equation} |
| 142 |
> |
The derivative of the above potential gives the following forces, |
| 143 |
> |
\begin{equation} |
| 144 |
> |
F^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2R_{\textrm{c}}^2)}}{R_{\textrm{c}}}\right]\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 145 |
> |
\end{cases}. |
| 146 |
> |
\label{eq:DSFForces} |
| 147 |
> |
\end{equation} |
| 148 |
|
|
| 149 |
+ |
This new Shifted-Force potential is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation \ref{eq:shiftingForm} is equal to equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in the Zahn potential, resulting in a discontinuity as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 150 |
+ |
|
| 151 |
|
\section{Methods} |
| 152 |
|
|
| 153 |
|
\subsection{What Qualities are Important?}\label{sec:Qualities} |
