--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/16 15:03:20	2630
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/18 05:31:17	2636
@@ -2,7 +2,7 @@
 %\documentclass[aps,prb,preprint]{revtex4}
 \documentclass[11pt]{article}
 \usepackage{endfloat}
-\usepackage{amsmath}
+\usepackage{amsmath,bm}
 \usepackage{amssymb}
 \usepackage{epsf}
 \usepackage{times}
@@ -81,8 +81,49 @@ impractical task to perform these calculations.
 impractical task to perform these calculations.
 
 \subsection{The Ewald Sum}
-blah blah blah Ewald Sum Important blah blah blah
+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,
+\begin{equation}
+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],
+\label{eq:PBCSum}
+\end{equation}
+where the sum over $\mathbf{n}$ is a sum over all periodic box replicas
+with integer coordinates $\mathbf{n} = (l,m,n)$, and the prime indicates
+$i = j$ are neglected for $\mathbf{n} = 0$.\cite{deLeeuw80} Within the
+sum, $N$ is the number of electrostatic particles, $\mathbf{r}_{ij}$ is
+$\mathbf{r}_j - \mathbf{r}_i$, $L$ is the cell length, $\bm{\Omega}_{i,j}$ are
+the Euler angles for $i$ and $j$, and $\phi$ is Poisson's equation
+($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge
+interactions). In the case of monopole electrostatics,
+eq. (\ref{eq:PBCSum}) is conditionally convergent and is discontiuous
+for non-neutral systems.
 
+This electrostatic summation problem was originally studied by Ewald
+for the case of an infinite crystal.\cite{Ewald21}. The approach he
+took was to convert this conditionally convergent sum into two
+absolutely convergent summations: a short-ranged real-space summation
+and a long-ranged reciprocal-space summation,
+\begin{equation}
+\begin{split}
+V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{3L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2,
+\end{split}
+\label{eq:EwaldSum}
+\end{equation}
+where $\alpha$ is a damping parameter, or separation constant, with
+units of \AA$^{-1}$, and $\mathbf{k}$ are the reciprocal vectors and
+equal $2\pi\mathbf{n}/L^2$. The final two terms of
+eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
+for interacting with a surrounding dielectric.\cite{Allen87} This
+dipolar term was neglected in early applications in molecular
+simulations,\cite{Brush66,Woodcock71} until it was introduced by de
+Leeuw {\it et al.} to address situations where the unit cell has a
+dipole moment and this dipole moment gets magnified through
+replication of the periodic images.\cite{deLeeuw80} This term is zero
+for systems where $\epsilon_{\rm S} = \infty$. Figure
+\ref{fig:ewaldTime} shows how the Ewald sum has been applied over
+time.  Initially, due to the small size of systems, the entire
+simulation box was replicated to convergence.  Currently, we balance a
+spherical real-space cutoff with the reciprocal sum and consider the
+surrounding dielectric.
 \begin{figure}
 \centering
 \includegraphics[width = \linewidth]{./ewaldProgression.pdf}
@@ -96,6 +137,18 @@ a surrounding dielectric term is included.}
 \label{fig:ewaldTime}
 \end{figure}
 
+The Ewald summation in the straight-forward form is an
+$\mathscr{O}(N^2)$ algorithm.  The separation constant $(\alpha)$
+plays an important role in the computational cost balance between the
+direct and reciprocal-space portions of the summation.  The choice of
+the magnitude of this value allows one to whether the real-space or
+reciprocal space portion of the summation is an $\mathscr{O}(N^2)$
+calcualtion, with the other being $\mathscr{O}(N)$.\cite{Sagui99} With
+appropriate choice of $\alpha$ and thoughtful algorithm development,
+this cost can be brought down to 
+$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to
+accelerate the Ewald summation is to se
+
 \subsection{The Wolf and Zahn Methods}
 In a recent paper by Wolf \textit{et al.}, a procedure was outlined
 for the accurate accumulation of electrostatic interactions in an
@@ -126,7 +179,7 @@ procedure gives an expression for the forces,
 derivative of this potential prior to evaluation of the limit.  This
 procedure gives an expression for the forces,
 \begin{equation}
-F_{\textrm{Wolf}}(r_{ij}) = q_i q_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\},
+F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{\left[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}\right]-\left[\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right]\right\},
 \label{eq:WolfForces}
 \end{equation}
 that incorporates both image charges and damping of the electrostatic
@@ -206,12 +259,12 @@ of the unshifted potential itself (when inside the cut
 The forces associated with the shifted potential are simply the forces
 of the unshifted potential itself (when inside the cutoff sphere),
 \begin{equation}
-F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right),
+f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
 \end{equation}
 and are zero outside.  Inside the cutoff sphere, the forces associated
 with the shifted force form can be written,
 \begin{equation}
-F_{\textrm{SF}} = \left( \frac{d v(r)}{dr} \right) - \left(\frac{d
+f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
 v(r)}{dr} \right)_{r=R_\textrm{c}}.
 \end{equation}
  
@@ -223,15 +276,15 @@ al.}'s undamped prescription:
 then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
 al.}'s undamped prescription:
 \begin{equation}
-V_\textrm{SP}(r) =
+v_\textrm{SP}(r) =
 q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
 r\leqslant R_\textrm{c},
-\label{eq:WolfSP}
+\label{eq:SPPot}
 \end{equation}
 with associated forces,
 \begin{equation}
-F_\textrm{SP}(r) = q_iq_j\left(-\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}.
-\label{eq:FWolfSP}
+f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}.
+\label{eq:SPForces}
 \end{equation}
 These forces are identical to the forces of the standard Coulomb
 interaction, and cutting these off at $R_c$ was addressed by Wolf
@@ -245,12 +298,12 @@ will give,
 The shifted force ({\sc sf}) form using the normal Coulomb potential
 will give,
 \begin{equation}
-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}.
+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}.
 \label{eq:SFPot}
 \end{equation}
 with associated forces,
 \begin{equation}
-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}.
+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}.
 \label{eq:SFForces}
 \end{equation}
 This formulation has the benefits that there are no discontinuities at
@@ -264,7 +317,7 @@ Wolf \textit{et al.} originally discussed the energeti
 to gain functionality in dynamics simulations.
 
 Wolf \textit{et al.} originally discussed the energetics of the
-shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that
+shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that
 it was still insufficient for accurate determination of the energy
 with reasonable cutoff distances.  The calculated Madelung energies
 fluctuate around the expected value with increasing cutoff radius, but
@@ -279,15 +332,15 @@ v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
 v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
 \label{eq:dampCoulomb}
 \end{equation}
-the shifted potential (Eq. (\ref{eq:WolfSP})) can be recovered
-using eq. (\ref{eq:shiftingForm}),
+the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using
+eq. (\ref{eq:shiftingForm}),
 \begin{equation}
-v_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r)}{r}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c},
+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},
 \label{eq:DSPPot}
 \end{equation}
 with associated forces,
 \begin{equation}
-f_{\textrm{DSP}}(r) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r)}{r^2}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \quad r\leqslant R_\textrm{c}.
+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}.
 \label{eq:DSPForces}
 \end{equation}
 Again, this damped shifted potential suffers from a discontinuity and
@@ -300,29 +353,31 @@ v_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc
 \label{eq:DSFPot}
 \end{split}
 \end{equation}
-The derivative of the above potential gives the following forces,
+The derivative of the above potential will lead to the following forces,
 \begin{equation}
 \begin{split}
-f_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \\ &\left.+\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] \quad r\leqslant R_\textrm{c}.
+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}.
 \label{eq:DSFForces}
 \end{split}
 \end{equation}
+If the damping parameter $(\alpha)$ is chosen to be zero, the undamped
+case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered
+from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}).
 
-This new {\sc sf} 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 eq. (\ref{eq:shiftingForm}) is equal to
-eq. (\ref{eq:dampCoulomb}) with $r$ replaced by $R_\textrm{c}$.  This
-term is {\it not} present in the Zahn potential, resulting in a
-potential discontinuity as particles cross $R_\textrm{c}$.  Second,
-the sign of the derivative portion is different.  The missing
-$v_\textrm{c}$ term would not affect molecular dynamics simulations
-(although the computed energy would be expected to have sudden jumps
-as particle distances crossed $R_c$).  The sign problem would be a
-potential source of errors, however.  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}$.  
+This new {\sc sf} 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
+eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
+with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
+in the Zahn potential, resulting in a potential discontinuity as
+particles cross $R_\textrm{c}$.  Second, the sign of the derivative
+portion is different.  The missing $v_\textrm{c}$ term would not
+affect molecular dynamics simulations (although the computed energy
+would be expected to have sudden jumps as particle distances crossed
+$R_c$).  The sign problem would be a potential source of errors,
+however.  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}$.
 
 Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
 electrostatic summation method that is continuous in both the
@@ -570,9 +625,10 @@ tolerance (typically less than $1 \times 10^{-4}$ kcal
 the energies and forces calculated.  Typical molecular mechanics
 packages default this to a value dependent on the cutoff radius and a
 tolerance (typically less than $1 \times 10^{-4}$ kcal/mol).  Smaller
-tolerances are typically associated with increased accuracy in the
-real-space portion of the summation.\cite{Essmann95} The default
-TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
+tolerances are typically associated with increased accuracy, but this
+usually means more time spent calculating the reciprocal-space portion
+of the summation.\cite{Perram88,Essmann95} The default TINKER
+tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
 calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and
 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.
 
@@ -851,7 +907,7 @@ required for accurate reproduction of crystal dynamics
 \begin{figure}
 \centering
 \includegraphics[width = \linewidth]{./comboSquare.pdf}
-\caption{Upper: Zoomed inset from figure \ref{fig:methodPS}.  As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift.  Lower: Low-frequency correlated motions for NaCl crystals at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$.  As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME.  The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.}
+\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods.  The upper plot is a zoomed inset from figure \ref{fig:methodPS}.  As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift.  The lower plot is of spectra when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$.  As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME.  The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.}
 \label{fig:dampInc}
 \end{figure}
 
