--- trunk/multipole/multipole_2/multipole2.tex	2014/06/12 14:58:06	4181
+++ trunk/multipole/multipole_2/multipole2.tex	2014/06/15 00:18:18	4185
@@ -69,43 +69,43 @@ of Notre Dame, Notre Dame, IN 46556}
 \begin{abstract}
   We have tested the real-space shifted potential (SP),
   gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods
-  for multipoles that were developed in the first paper in this series
-  against a reference method. The tests were carried out in a variety
-  of condensed-phase environments which were designed to test all
-  levels of the multipole-multipole interactions.  Comparisons of the
-  energy differences between configurations, molecular forces, and
-  torques were used to analyze how well the real-space models perform
-  relative to the more computationally expensive Ewald sum.  We have
-  also investigated the energy conservation properties of the new
-  methods in molecular dynamics simulations using all of these
-  methods. The SP method shows excellent agreement with
-  configurational energy differences, forces, and torques, and would
-  be suitable for use in Monte Carlo calculations.  Of the two new
-  shifted-force methods, the GSF approach shows the best agreement
-  with Ewald-derived energies, forces, and torques and exhibits energy
-  conservation properties that make it an excellent choice for
-  efficiently computing electrostatic interactions in molecular
-  dynamics simulations.
+  for multipole interactions that were developed in the first paper in
+  this series, using the multipolar Ewald sum as a reference
+  method. The tests were carried out in a variety of condensed-phase
+  environments which were designed to test all levels of the
+  multipole-multipole interactions.  Comparisons of the energy
+  differences between configurations, molecular forces, and torques
+  were used to analyze how well the real-space models perform relative
+  to the more computationally expensive Ewald treatment.  We have also
+  investigated the energy conservation properties of the new methods
+  in molecular dynamics simulations. The SP method shows excellent
+  agreement with configurational energy differences, forces, and
+  torques, and would be suitable for use in Monte Carlo calculations.
+  Of the two new shifted-force methods, the GSF approach shows the
+  best agreement with Ewald-derived energies, forces, and torques and
+  exhibits energy conservation properties that make it an excellent
+  choice for efficient computation of electrostatic interactions in
+  molecular dynamics simulations.
 \end{abstract}
 
 %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                              % Classification Scheme.
-\keywords{Electrostatics, Multipoles, Real-space}
+%\keywords{Electrostatics, Multipoles, Real-space}
 
 \maketitle
 
 
 \section{\label{sec:intro}Introduction}
 Computing the interactions between electrostatic sites is one of the
-most expensive aspects of molecular simulations, which is why there
-have been significant efforts to develop practical, efficient and
-convergent methods for handling these interactions. Ewald's method is
-perhaps the best known and most accurate method for evaluating
-energies, forces, and torques in explicitly-periodic simulation
-cells. In this approach, the conditionally convergent electrostatic
-energy is converted into two absolutely convergent contributions, one
-which is carried out in real space with a cutoff radius, and one in
-reciprocal space.\cite{Clarke:1986eu,Woodcock75}
+most expensive aspects of molecular simulations. There have been
+significant efforts to develop practical, efficient and convergent
+methods for handling these interactions. Ewald's method is perhaps the
+best known and most accurate method for evaluating energies, forces,
+and torques in explicitly-periodic simulation cells. In this approach,
+the conditionally convergent electrostatic energy is converted into
+two absolutely convergent contributions, one which is carried out in
+real space with a cutoff radius, and one in reciprocal
+space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75}
 
 When carried out as originally formulated, the reciprocal-space
 portion of the Ewald sum exhibits relatively poor computational
@@ -116,44 +116,74 @@ Because of the artificial periodicity required for the
 the computational cost from $O(N^2)$ down to $O(N \log
 N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}.
 
-Because of the artificial periodicity required for the Ewald sum, the
-method may require modification to compute interactions for
+Because of the artificial periodicity required for the Ewald sum,
 interfacial molecular systems such as membranes and liquid-vapor
-interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
-To simulate interfacial systems, Parry's extension of the 3D Ewald sum
-is appropriate for slab geometries.\cite{Parry:1975if} The inherent
-periodicity in the Ewald’s method can also be problematic for
-interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald
-methods that were developed to handle two-dimensional (2D)
-electrostatic interactions in interfacial systems have not had similar
-particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77,
-  Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq}
+interfaces require modifications to the
+method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
+Parry's extension of the three dimensional Ewald sum is appropriate
+for slab geometries.\cite{Parry:1975if} Modified Ewald methods that
+were developed to handle two-dimensional (2D) electrostatic
+interactions in interfacial systems have not seen similar
+particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77,
+  Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly
+with system size. The inherent periodicity in the Ewald’s method can
+also be problematic for interfacial molecular
+systems.\cite{Fennell:2006lq}
 
 \subsection{Real-space methods}
 Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$
 method for calculating electrostatic interactions between point
-charges. They argued that the effective Coulomb interaction in
-condensed systems is actually short ranged.\cite{Wolf92,Wolf95} For an
-ordered lattice (e.g., when computing the Madelung constant of an
-ionic solid), the material can be considered as a set of ions
-interacting with neutral dipolar or quadrupolar ``molecules'' giving
-an effective distance dependence for the electrostatic interactions of
-$r^{-5}$ (see figure \ref{fig:schematic}).  For this reason, careful
-applications of Wolf's method are able to obtain accurate estimates of
-Madelung constants using relatively short cutoff radii.  Recently,
-Fukuda used neutralization of the higher order moments for the
-calculation of the electrostatic interaction of the point charges
-system.\cite{Fukuda:2013sf}
+charges. They argued that the effective Coulomb interaction in most
+condensed phase systems is effectively short
+ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when
+computing the Madelung constant of an ionic solid), the material can
+be considered as a set of ions interacting with neutral dipolar or
+quadrupolar ``molecules'' giving an effective distance dependence for
+the electrostatic interactions of $r^{-5}$ (see figure
+\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple
+cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the
+electrostatic energy per ion converges more rapidly to the Madelung
+energy than the dipolar approximation.\cite{Wolf92} To find the
+correct Madelung constant, Lacman suggested that the NaCl structure
+could be constructed in a way that the finite crystal terminates with
+complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees
+what is effectively a set of octupoles at large distances. These facts
+suggest that the Madelung constants are relatively short ranged for
+perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful
+application of Wolf's method are able to obtain accurate estimates of
+Madelung constants using relatively short cutoff radii.
 
+Direct truncation of interactions at a cutoff radius creates numerical
+errors.  Wolf \textit{et al.}  argued that truncation errors are due
+to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To
+neutralize this charge they proposed placing an image charge on the
+surface of the cutoff sphere for every real charge inside the cutoff.
+These charges are present for the evaluation of both the pair
+interaction energy and the force, although the force expression
+maintained a discontinuity at the cutoff sphere.  In the original Wolf
+formulation, the total energy for the charge and image were not equal
+to the integral of their force expression, and as a result, the total
+energy would not be conserved in molecular dynamics (MD)
+simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and
+Gezelter later proposed shifted force variants of the Wolf method with
+commensurate force and energy expressions that do not exhibit this
+problem.\cite{Fennell:2006lq} Related real-space methods were also
+proposed by Chen \textit{et
+  al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw}
+and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has used
+neutralization of the higher order moments for the calculation of the
+electrostatic interaction of the point charge
+systems.\cite{Fukuda:2013sf}
+
 \begin{figure}
   \centering
   \includegraphics[width=\linewidth]{schematic.pdf}
   \caption{Top: Ionic systems exhibit local clustering of dissimilar
     charges (in the smaller grey circle), so interactions are
-    effectively charge-multipole in order at longer distances.  With
-    hard cutoffs, motion of individual charges in and out of the
-    cutoff sphere can break the effective multipolar ordering.
-    Bottom: dipolar crystals and fluids have a similar effective
+    effectively charge-multipole at longer distances.  With hard
+    cutoffs, motion of individual charges in and out of the cutoff
+    sphere can break the effective multipolar ordering.  Bottom:
+    dipolar crystals and fluids have a similar effective
     \textit{quadrupolar} ordering (in the smaller grey circles), and
     orientational averaging helps to reduce the effective range of the
     interactions in the fluid.  Placement of reversed image multipoles
@@ -162,46 +192,16 @@ The direct truncation of interactions at a cutoff radi
   \label{fig:schematic}
 \end{figure}
 
-The direct truncation of interactions at a cutoff radius creates
-truncation defects. Wolf \textit{et al.} further argued that
-truncation errors are due to net charge remaining inside the cutoff
-sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed
-placing an image charge on the surface of the cutoff sphere for every
-real charge inside the cutoff.  These charges are present for the
-evaluation of both the pair interaction energy and the force, although
-the force expression maintained a discontinuity at the cutoff sphere.
-In the original Wolf formulation, the total energy for the charge and
-image were not equal to the integral of their force expression, and as
-a result, the total energy would not be conserved in molecular
-dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and
-Fennel and Gezelter later proposed shifted force variants of the Wolf
-method with commensurate force and energy expressions that do not
-exhibit this problem.\cite{Fennell:2006lq}   Related real-space
-methods were also proposed by Chen \textit{et
-  al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw}
-and by Wu and Brooks.\cite{Wu:044107}
-
-Considering the interaction of one central ion in an ionic crystal
-with a portion of the crystal at some distance, the effective Columbic
-potential is found to decrease as $r^{-5}$. If one views the \ce{NaCl}
-crystal as a simple cubic (SC) structure with an octupolar
-\ce{(NaCl)4} basis, the electrostatic energy per ion converges more
-rapidly to the Madelung energy than the dipolar
-approximation.\cite{Wolf92} To find the correct Madelung constant,
-Lacman suggested that the NaCl structure could be constructed in a way
-that the finite crystal terminates with complete \ce{(NaCl)4}
-molecules.\cite{Lacman65} The central ion sees what is effectively a
-set of octupoles at large distances. These facts suggest that the
-Madelung constants are relatively short ranged for perfect ionic
-crystals.\cite{Wolf:1999dn}
-
-One can make a similar argument for crystals of point multipoles. The
-Luttinger and Tisza treatment of energy constants for dipolar lattices
-utilizes 24 basis vectors that contain dipoles at the eight corners of
-a unit cube.  Only three of these basis vectors, $X_1, Y_1,
-\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have
-zero net dipole and retain contributions only from higher order
-multipoles.  The effective interaction between a dipole at the center
+One can make a similar effective range argument for crystals of point
+\textit{multipoles}. The Luttinger and Tisza treatment of energy
+constants for dipolar lattices utilizes 24 basis vectors that contain
+dipoles at the eight corners of a unit cube.  Only three of these
+basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole
+moments, while the rest have zero net dipole and retain contributions
+only from higher order multipoles.  The lowest energy crystalline
+structures are built out of basis vectors that have only residual
+quadrupolar moments (e.g. the $Z_5$ array). In these low energy
+structures, the effective interaction between a dipole at the center
 of a crystal and a group of eight dipoles farther away is
 significantly shorter ranged than the $r^{-3}$ that one would expect
 for raw dipole-dipole interactions.  Only in crystals which retain a
@@ -258,48 +258,42 @@ The damping function used in our research has been dis
 densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13}
 
 \subsection{The damping function}
-The damping function used in our research has been discussed in detail
-in the first paper of this series.\cite{PaperI} The radial kernel
-$1/r$ for the interactions between point charges can be replaced by
-the complementary error function $\mathrm{erfc}(\alpha r)/r$ to
-accelerate the rate of convergence, where $\alpha$ is a damping
-parameter with units of inverse distance.  Altering the value of
-$\alpha$ is equivalent to changing the width of Gaussian charge
-distributions that replace each point charge -- Gaussian overlap
-integrals yield complementary error functions when truncated at a
-finite distance.
+The damping function has been discussed in detail in the first paper
+of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the
+interactions between point charges can be replaced by the
+complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate
+convergence, where $\alpha$ is a damping parameter with units of
+inverse distance.  Altering the value of $\alpha$ is equivalent to
+changing the width of Gaussian charge distributions that replace each
+point charge, as Coulomb integrals with Gaussian charge distributions
+produce complementary error functions when truncated at a finite
+distance.
 
-By using suitable value of damping alpha ($\alpha \sim 0.2$) for a
-cutoff radius ($r_{­c}=9 A$), Fennel and Gezelter produced very good
-agreement with SPME for the interaction energies, forces and torques
-for charge-charge interactions.\cite{Fennell:2006lq}
+With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method
+produced very good agreement with SPME for interaction energies,
+forces and torques for charge-charge
+interactions.\cite{Fennell:2006lq}
 
 \subsection{Point multipoles in molecular modeling}
 Coarse-graining approaches which treat entire molecular subsystems as
 a single rigid body are now widely used. A common feature of many
 coarse-graining approaches is simplification of the electrostatic
 interactions between bodies so that fewer site-site interactions are
-required to compute configurational energies.  Many coarse-grained
-molecular structures would normally consist of equal positive and
-negative charges, and rather than use multiple site-site interactions,
-the interaction between higher order multipoles can also be used to
-evaluate a single molecule-molecule
-interaction.\cite{Ren06,Essex10,Essex11}
+required to compute configurational
+energies.\cite{Ren06,Essex10,Essex11}
 
 Because electrons in a molecule are not localized at specific points,
-the assignment of partial charges to atomic centers is a relatively
-rough approximation.  Atomic sites can also be assigned point
-multipoles and polarizabilities to increase the accuracy of the
-molecular model.  Recently, water has been modeled with point
-multipoles up to octupolar order using the soft sticky
-dipole-quadrupole-octupole (SSDQO)
+the assignment of partial charges to atomic centers is always an
+approximation.  Atomic sites can also be assigned point multipoles and
+polarizabilities to increase the accuracy of the molecular model.
+Recently, water has been modeled with point multipoles up to octupolar
+order using the soft sticky dipole-quadrupole-octupole (SSDQO)
 model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point
 multipoles up to quadrupolar order have also been coupled with point
 polarizabilities in the high-quality AMOEBA and iAMOEBA water
-models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But
-using point multipole with the real space truncation without
-accounting for multipolar neutrality will create energy conservation
-issues in molecular dynamics (MD) simulations.
+models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However,
+truncating point multipoles without smoothing the forces and torques
+will create energy conservation issues in molecular dynamics simulations.
 
 In this paper we test a set of real-space methods that were developed
 for point multipolar interactions.  These methods extend the damped
@@ -339,7 +333,7 @@ $\bf a$.
 where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is
 expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf
     a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object
-$\bf a$.
+$\bf a$, etc.
 
 % Interactions between multipoles can be expressed as higher derivatives
 % of the bare Coulomb potential, so one way of ensuring that the forces
@@ -412,7 +406,7 @@ to another site within cutoff sphere are derived from
 connection to unmodified electrostatics as well as the smooth
 transition to zero in both these functions as $r\rightarrow r_c$.  The
 electrostatic forces and torques acting on the central multipole due
-to another site within cutoff sphere are derived from
+to another site within the cutoff sphere are derived from
 Eq.~\ref{generic}, accounting for the appropriate number of
 derivatives. Complete energy, force, and torque expressions are
 presented in the first paper in this series (Reference
@@ -430,7 +424,7 @@ U_{D_{\bf a} D_{\bf b}}(r_c)
 U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r) -
 U_{D_{\bf a} D_{\bf b}}(r_c)
    - (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot
-  \vec{\nabla} U_{D_{\bf a}D_{\bf b}}(r) \Big \lvert _{r_c}
+  \nabla U_{D_{\bf a}D_{\bf b}}(r_c).
 \end{equation}
 Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf
   a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance
@@ -454,18 +448,21 @@ U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathb
 In general, the gradient shifted potential between a central multipole
 and any multipolar site inside the cutoff radius is given by,
 \begin{equation}
-U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
-U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r}
-\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert  _{r_c} \right]
+  U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
+    U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}}
+    \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
 \label{generic2}
 \end{equation}
 where the sum describes a separate force-shifting that is applied to
-each orientational contribution to the energy. 
+each orientational contribution to the energy.  In this expression,
+$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles
+($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$
+represent the orientations the multipoles.
 
 The third term converges more rapidly than the first two terms as a
 function of radius, hence the contribution of the third term is very
 small for large cutoff radii.  The force and torque derived from
-equation \ref{generic2} are consistent with the energy expression and
+Eq. \ref{generic2} are consistent with the energy expression and
 approach zero as $r \rightarrow r_c$.  Both the GSF and TSF methods
 can be considered generalizations of the original DSF method for
 higher order multipole interactions. GSF and TSF are also identical up
@@ -473,7 +470,7 @@ GSF potential are presented in the first paper in this
 the energy, force and torque for higher order multipole-multipole
 interactions. Complete energy, force, and torque expressions for the
 GSF potential are presented in the first paper in this series
-(Reference~\onlinecite{PaperI})
+(Reference~\onlinecite{PaperI}).
 
 
 \subsection{Shifted potential (SP) }
@@ -487,7 +484,7 @@ U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \ri
 effectively shifts the total potential to zero at the cutoff radius,
 \begin{equation}
 U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
-U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
+U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
 \label{eq:SP}
 \end{equation}      	
 where the sum describes separate potential shifting that is done for
