--- trunk/multipole/multipole_2/multipole2.tex	2014/06/14 23:35:39	4184
+++ trunk/multipole/multipole_2/multipole2.tex	2014/08/07 21:13:00	4207
@@ -47,45 +47,41 @@ preprint,
 
 %\preprint{AIP/123-QED}
 
-\title{Real space alternatives to the Ewald
-Sum. II. Comparison of Methods} % Force line breaks with \\
+\title{Real space electrostatics for multipoles. II. Comparisons with
+  the Ewald Sum}
 
 \author{Madan Lamichhane}
- \affiliation{Department of Physics, University
-of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\
+ \affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556}
 
 \author{Kathie E. Newman}
-\affiliation{Department of Physics, University
-of Notre Dame, Notre Dame, IN 46556}
+\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556}
 
 \author{J. Daniel Gezelter}%
  \email{gezelter@nd.edu.}
-\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash
-}%
+\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556
+}
 
-\date{\today}% It is always \today, today,
-             %  but any date may be explicitly specified
+\date{\today}
 
 \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 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 using all of these methods. The SP method shows
+  We report on tests of the shifted potential (SP), gradient shifted
+  force (GSF), and Taylor shifted force (TSF) real-space methods for
+  multipole interactions 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 designed to
+  test up to quadrupole-quadrupole 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 efficiently computing electrostatic
-  interactions in molecular dynamics simulations.
+  forces, and torques and also 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
@@ -94,61 +90,91 @@ of Notre Dame, Notre Dame, IN 46556}
 
 \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.\cite{Ewald21,deLeeuw80,Smith81,Allen87} 
 
 When carried out as originally formulated, the reciprocal-space
 portion of the Ewald sum exhibits relatively poor computational
-scaling, making it prohibitive for large systems. By utilizing
-particle meshes and three dimensional fast Fourier transforms (FFT),
-the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald
-(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease
-the computational cost from $O(N^2)$ down to $O(N \log
-N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}.
+scaling, making it prohibitive for large systems. By utilizing a
+particle mesh and three dimensional fast Fourier transforms (FFT), the
+particle-mesh Ewald (PME), particle-particle particle-mesh Ewald
+(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME)
+methods can decrease 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 method can also be problematic for molecular
-interfaces.\cite{Fennell:2006lq} 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.
+interfaces require modifications to the method.  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 (2-D) electrostatic
+interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
+These methods were originally quite computationally
+expensive.\cite{Spohr97,Yeh99} There have been several successful
+efforts that reduced the computational cost of 2-D lattice summations,
+bringing them more in line with the scaling for the full 3-D
+treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The
+inherent periodicity required by the Ewald method can also be
+problematic in a number of protein/solvent and ionic solution
+environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,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 phase 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
-application of Wolf's method can obtain accurate estimates of Madelung
-constants using relatively short cutoff radii.  Recently, Fukuda used
-neutralization of the higher order moments for calculation of the
-electrostatic interactions in point charge
-systems.\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 can provide 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.} suggest 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
+maintains 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 the 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{Zahn:2002hc,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 successfuly
+used additional neutralization of higher order moments for systems of
+point charges.\cite{Fukuda:2013sf}
+
 \begin{figure}
   \centering
-  \includegraphics[width=\linewidth]{schematic.pdf}
+  \includegraphics[width=\linewidth]{schematic.eps}
   \caption{Top: Ionic systems exhibit local clustering of dissimilar
     charges (in the smaller grey circle), so interactions are
     effectively charge-multipole at longer distances.  With hard
@@ -163,46 +189,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.}  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.\cite{LT} 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
@@ -222,12 +218,12 @@ Even at elevated temperatures, there is, on average, l
 
 The shorter effective range of electrostatic interactions is not
 limited to perfect crystals, but can also apply in disordered fluids.
-Even at elevated temperatures, there is, on average, local charge
-balance in an ionic liquid, where each positive ion has surroundings
-dominated by negaitve ions and vice versa.  The reversed-charge images
-on the cutoff sphere that are integral to the Wolf and DSF approaches
-retain the effective multipolar interactions as the charges traverse
-the cutoff boundary.
+Even at elevated temperatures, there is local charge balance in an
+ionic liquid, where each positive ion has surroundings dominated by
+negative ions and vice versa.  The reversed-charge images on the
+cutoff sphere that are integral to the Wolf and DSF approaches retain
+the effective multipolar interactions as the charges traverse the
+cutoff boundary.
 
 In multipolar fluids (see Fig. \ref{fig:schematic}) there is
 significant orientational averaging that additionally reduces the
@@ -246,72 +242,66 @@ The forces and torques acting on atomic sites are the 
 % to the non-neutralized value of the higher order moments within the
 % cutoff sphere.
 
-The forces and torques acting on atomic sites are the fundamental
-factors driving dynamics in molecular simulations. Fennell and
-Gezelter proposed the damped shifted force (DSF) energy kernel to
-obtain consistent energies and forces on the atoms within the cutoff
-sphere. Both the energy and the force go smoothly to zero as an atom
-aproaches the cutoff radius. The comparisons of the accuracy these
-quantities between the DSF kernel and SPME was surprisingly
-good.\cite{Fennell:2006lq} The DSF method has seen increasing use for
-calculating electrostatic interactions in molecular systems with
-relatively uniform charge
-densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13}
+Forces and torques acting on atomic sites are fundamental in driving
+dynamics in molecular simulations, and the damped shifted force (DSF)
+energy kernel provides consistent energies and forces on charged atoms
+within the cutoff sphere. Both the energy and the force go smoothly to
+zero as an atom aproaches the cutoff radius. The comparisons of the
+accuracy these quantities between the DSF kernel and SPME was
+surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method
+has seen increasing use in molecular systems with relatively uniform
+charge
+densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,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)
+Additionally, because electrons in a molecule are not localized at
+specific points, the assignment of partial charges to atomic centers
+is always an approximation.  For increased accuracy, atomic sites can
+also be assigned point multipoles and polarizabilities.  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
+can 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
 shifted force (DSF) and Wolf methods originally developed for
 charge-charge interactions and generalize them for higher order
-multipoles. The detailed mathematical development of these methods has
-been presented in the first paper in this series, while this work
-covers the testing the energies, forces, torques, and energy
+multipoles.  The detailed mathematical development of these methods
+has been presented in the first paper in this series, while this work
+covers the testing of energies, forces, torques, and energy
 conservation properties of the methods in realistic simulation
 environments.  In all cases, the methods are compared with the
-reference method, a full multipolar Ewald treatment.
+reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98}
 
 
 %\subsection{Conservation of total energy } 
@@ -335,12 +325,11 @@ U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1
 expressed as the product of two multipole operators and a Coulombic
 kernel,
 \begin{equation}
-U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
+U_{ab}(r)= M_{a} M_{b} \frac{1}{r}  \label{kernel}.
 \end{equation}
-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$, etc.
+where the multipole operator for site $a$, $M_{a}$, is
+expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object
+$a$, etc.
 
 % Interactions between multipoles can be expressed as higher derivatives
 % of the bare Coulomb potential, so one way of ensuring that the forces
@@ -395,15 +384,15 @@ radius.  For example, the direct dipole dot product
 contributions to the potential, and ensures that the forces and
 torques from each of these contributions will vanish at the cutoff
 radius.  For example, the direct dipole dot product
-($\mathbf{D}_{\bf a}
-\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance
+($\mathbf{D}_{a}
+\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance
 dot products:
 \begin{equation}
-U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
-  \mathbf{D}_{\bf a} \cdot
-\mathbf{D}_{\bf b} \right) v_{21}(r) +
-\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right)
-\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right]
+U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
+  \mathbf{D}_{a} \cdot
+\mathbf{D}_{b} \right) v_{21}(r) +
+\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right)
+\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right]
 \end{equation}
 
 For the Taylor shifted (TSF) method with the undamped kernel,
@@ -428,13 +417,12 @@ U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r)
 which have been projected onto the surface of the cutoff sphere
 without changing their relative orientation,
 \begin{equation}
-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
-  \nabla U_{D_{\bf a}D_{\bf b}}(r_c).
+U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
+U_{D_{a}D_{b}}(r_c)
+   - (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot
+  \nabla U_{D_{a}D_{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
+Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance
 (although the signs are reversed for the dipole that has been
 projected onto the cutoff sphere).  In many ways, this simpler
 approach is closer in spirit to the original shifted force method, in
@@ -455,15 +443,15 @@ and any multipolar site inside the cutoff radius is gi
 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(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]
+U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
+U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c)
+\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
 \label{generic2}
 \end{equation}
 where the sum describes a separate force-shifting that is applied to
 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}}$
+($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$
 represent the orientations the multipoles.
 
 The third term converges more rapidly than the first two terms as a
@@ -490,8 +478,8 @@ U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf
 interactions with the central multipole and the image. This
 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(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
+U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
+U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
 \label{eq:SP}
 \end{equation}      	
 where the sum describes separate potential shifting that is done for
@@ -536,7 +524,7 @@ in the test cases are given in table~\ref{tab:pars}.
 used the multipolar Ewald sum as a reference method for comparing
 energies, forces, and torques for molecular models that mimic
 disordered and ordered condensed-phase systems.  The parameters used
-in the test cases are given in table~\ref{tab:pars}.
+in the test cases are given in table~\ref{tab:pars}. 
 
 \begin{table}
 \label{tab:pars}
@@ -584,7 +572,7 @@ the simulation proceeds. These differences are the mos
 and have been compared with the values obtained from the multipolar
 Ewald sum.  In Monte Carlo (MC) simulations, the energy differences
 between two configurations is the primary quantity that governs how
-the simulation proceeds. These differences are the most imporant
+the simulation proceeds. These differences are the most important
 indicators of the reliability of a method even if the absolute
 energies are not exact.  For each of the multipolar systems listed
 above, we have compared the change in electrostatic potential energy
@@ -596,7 +584,7 @@ program, OpenMD,\cite{openmd} which was used for all c
 \subsection{Implementation}
 The real-space methods developed in the first paper in this series
 have been implemented in our group's open source molecular simulation
-program, OpenMD,\cite{openmd} which was used for all calculations in
+program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in
 this work.  The complementary error function can be a relatively slow
 function on some processors, so all of the radial functions are
 precomputed on a fine grid and are spline-interpolated to provide
@@ -606,10 +594,13 @@ approximations.\cite{Smith82,Smith98} In all cases, th
 with a reciprocal space cutoff, $k_\mathrm{max} = 7$.  Our version of
 the Ewald sum is a re-implementation of the algorithm originally
 proposed by Smith that does not use the particle mesh or smoothing
-approximations.\cite{Smith82,Smith98} In all cases, the quantities
-being compared are the electrostatic contributions to energies, force,
-and torques.  All other contributions to these quantities (i.e. from
-Lennard-Jones interactions) are removed prior to the comparisons.
+approximations.\cite{Smith82,Smith98} This implementation was tested
+extensively against the analytic energy constants for the multipolar
+lattices that are discussed in reference \onlinecite{PaperI}.  In all
+cases discussed below, the quantities being compared are the
+electrostatic contributions to energies, force, and torques.  All
+other contributions to these quantities (i.e. from Lennard-Jones
+interactions) are removed prior to the comparisons.
 
 The convergence parameter ($\alpha$) also plays a role in the balance
 of the real-space and reciprocal-space portions of the Ewald
@@ -803,7 +794,7 @@ model must allow for long simulation times with minima
  
 \begin{figure}
   \centering
-  \includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf}
+  \includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps}
   \caption{Statistical analysis of the quality of configurational
     energy differences for the real-space electrostatic methods
     compared with the reference Ewald sum.  Results with a value equal
@@ -876,7 +867,7 @@ forces is desired.
 
 \begin{figure}
   \centering
-  \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} 
+  \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} 
   \caption{Statistical analysis of the quality of the force vector
     magnitudes for the real-space electrostatic methods compared with
     the reference Ewald sum. Results with a value equal to 1 (dashed
@@ -890,7 +881,7 @@ forces is desired.
 
 \begin{figure}
   \centering
-  \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf}
+  \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps}
   \caption{Statistical analysis of the quality of the torque vector
     magnitudes for the real-space electrostatic methods compared with
     the reference Ewald sum. Results with a value equal to 1 (dashed
@@ -948,7 +939,7 @@ systematically improved by varying $\alpha$ and $r_c$.
 
 \begin{figure}
   \centering
-  \includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf}
+  \includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps}
   \caption{The circular variance of the direction of the force and
     torque vectors obtained from the real-space methods around the
     reference Ewald vectors. A variance equal to 0 (dashed line)
@@ -968,13 +959,14 @@ temperature of 300K.  After equilibration, this liquid
 in this series and provides the most comprehensive test of the new
 methods.  A liquid-phase system was created with 2000 water molecules
 and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a
-temperature of 300K.  After equilibration, this liquid-phase system
-was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with
-a cutoff radius of 12\AA.  The value of the damping coefficient was
-also varied from the undamped case ($\alpha = 0$) to a heavily damped
-case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods.  A
-sample was also run using the multipolar Ewald sum with the same
-real-space cutoff.
+temperature of 300K.  After equilibration in the canonical (NVT)
+ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system
+was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald,
+Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA.  The
+value of the damping coefficient was also varied from the undamped
+case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$
+\AA$^{-1}$) for all of the real space methods.  A sample was also run
+using the multipolar Ewald sum with the same real-space cutoff.
 
 In figure~\ref{fig:energyDrift} we show the both the linear drift in
 energy over time, $\delta E_1$, and the standard deviation of energy
@@ -995,19 +987,169 @@ $k$-space cutoff values.
 
 \begin{figure}
   \centering
-  \includegraphics[width=\textwidth]{newDrift_12.pdf}
+  \includegraphics[width=\textwidth]{newDrift_12.eps}
 \label{fig:energyDrift}        
 \caption{Analysis of the energy conservation of the real-space
-  electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in
-  energy over time (in kcal / mol / particle / ns) and $\delta
-  \mathrm{E}_0$ is the standard deviation of energy fluctuations
-  around this drift (in kcal / mol / particle).  All simulations were
-  of a 2000-molecule simulation of SSDQ water with 48 ionic charges at
+  methods. $\delta \mathrm{E}_1$ is the linear drift in energy over
+  time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is
+  the standard deviation of energy fluctuations around this drift (in
+  kcal / mol / particle).  Points that appear below the dashed grey
+  (Ewald) lines exhibit better energy conservation than commonly-used
+  parameters for Ewald-based electrostatics.  All simulations were of
+  a 2000-molecule simulation of SSDQ water with 48 ionic charges at
   300 K starting from the same initial configuration. All runs
   utilized the same real-space cutoff, $r_c = 12$\AA.}
+\end{figure} 
+
+\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}}
+The most important test of the modified interaction potentials is the
+fidelity with which they can reproduce structural features and
+dynamical properties in a liquid.  One commonly-utilized measure of
+structural ordering is the pair distribution function, $g(r)$, which
+measures local density deviations in relation to the bulk density.  In
+the electrostatic approaches studied here, the short-range repulsion
+from the Lennard-Jones potential is identical for the various
+electrostatic methods, and since short range repulsion determines much
+of the local liquid ordering, one would not expect to see many
+differences in $g(r)$.  Indeed, the pair distributions are essentially
+identical for all of the electrostatic methods studied (for each of
+the different systems under investigation).  An example of this
+agreement for the SSDQ water/ion system is shown in
+Fig. \ref{fig:gofr}.
+
+\begin{figure}
+  \centering
+  \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
+\label{fig:gofr}        
+\caption{The pair distribution functions, $g(r)$, for the SSDQ
+  water/ion system obtained using the different real-space methods are
+  essentially identical with the result from the Ewald
+  treatment.}
 \end{figure} 
 
+There is a very slight overstructuring of the first solvation shell
+when using when using TSF at lower values of the damping coefficient
+($\alpha \le 0.1$) or when using undamped GSF.  With moderate damping,
+GSF and SP produce pair distributions that are identical (within
+numerical noise) to their Ewald counterparts.
 
+A structural property that is a more demanding test of modified
+electrostatics is the mean value of the electrostatic energy $\langle
+U_\mathrm{elect} \rangle / N$ which is obtained by sampling the
+liquid-state configurations experienced by a liquid evolving entirely
+under the influence of each of the methods.  In table \ref{tab:Props}
+we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with
+the damping parameter, $\alpha$, for each of the methods.
+
+As in the crystals studied in the first paper, damping is important
+for converging the mean electrostatic energy values, particularly for
+the two shifted force methods (GSF and TSF).  A value of $\alpha
+\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF
+energies with a cutoff of 12 \AA, while shorter cutoffs require more
+dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA).
+Overdamping the real-space electrostatic methods occurs with $\alpha >
+0.4$, causing the estimate of the energy to drop below the Ewald
+results.
+
+These ``optimal'' values of the damping coefficient are slightly
+larger than what were observed for DSF electrostatics for purely
+point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for
+$r_c = 12$\AA appears to be an excellent compromise for mixed charge
+multipole systems.
+
+To test the fidelity of the electrostatic methods at reproducing
+dynamics in a multipolar liquid, it is also useful to look at
+transport properties, particularly the diffusion constant,
+\begin{equation}
+D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left|
+  \mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle
+\label{eq:diff}
+\end{equation}
+which measures long-time behavior and is sensitive to the forces on
+the multipoles.  For the soft dipolar fluid and the SSDQ liquid
+systems, the self-diffusion constants (D) were calculated from linear
+fits to the long-time portion of the mean square displacement,
+$\langle r^{2}(t) \rangle$.\cite{Allen87}
+
+In addition to translational diffusion, orientational relaxation times
+were calculated for comparisons with the Ewald simulations and with
+experiments. These values were determined from the same 1~ns
+microcanonical trajectories used for translational diffusion by
+calculating the orientational time correlation function,
+\begin{equation}
+C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t)
+                \cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle,
+\label{eq:OrientCorr}
+\end{equation}
+where $P_l$ is the Legendre polynomial of order $l$ and
+$\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis
+$\gamma$ on a molecule..  Th body-fixed reference frame used for our
+models has the $z$-axis running along the dipoles, and for the SSDQ
+water model, the $y$-axis connects the two implied hydrogen atom
+positions.  From the orientation autocorrelation functions, we can
+obtain time constants for rotational relaxation either by fitting an
+exponential function or by integrating the entire correlation
+function.  In a good water model, these decay times would be
+comparable to water orientational relaxation times from nuclear
+magnetic resonance (NMR). The relaxation constant obtained from
+$C_2^y(t)$ is normally of experimental interest because it describes
+the relaxation of the principle axis connecting the hydrogen
+atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
+of the dipole-dipole relaxation from a proton NMR signal and should
+provide an estimate of the NMR relaxation time constant.\cite{Impey82}
+
+Results for the diffusion constants and orientational relaxation times
+are shown in figure \ref{tab:Props}. From this data, it is apparent
+that the values for both $D$ and $\tau_2$ using the Ewald sum are
+reproduced with reasonable fidelity by the GSF method.
+
+The $\tau_2$ results in \ref{tab:Props} show a much greater difference
+between the real-space and the Ewald results.
+
+\begin{table}
+\label{tab:Props}
+\caption{Comparison of the structural and dynamic properties for the
+  soft dipolar liquid test for all of the real-space methods.}
+\begin{tabular}{l|c|cccc|cccc|cccc}
+         & Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\
+$\alpha$ (\AA$^{-1}$) & &      
+ 0.0 & 0.1 & 0.2 & 0.3 &
+ 0.0 & 0.1 & 0.2 & 0.3 &
+ 0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14}
+
+$\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\
+D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)&
+470.2(6) &
+416.6(5) &
+379.6(5) &
+438.6(5) &
+476.0(6) &
+412.8(5) &
+421.1(5) &
+400.5(5) &
+437.5(6) &
+434.6(5) &
+411.4(5) &
+545.3(7) &
+459.6(6) \\
+$\tau_2$ (fs) &
+1.136 &
+1.041 &
+1.064 &
+1.109 &
+1.211 &
+1.119 &
+1.039 &
+1.058 &
+1.21  &
+1.15  &
+1.172 &
+1.153 &
+1.125 \\
+\end{tabular}
+\end{table}
+
+
 \section{CONCLUSION}
 In the first paper in this series, we generalized the
 charge-neutralized electrostatic energy originally developed by Wolf