--- trunk/multipole/multipole_2/multipole2.tex	2014/08/18 14:54:19	4212
+++ trunk/multipole/multipole_2/multipole2.tex	2014/09/10 21:06:56	4214
@@ -165,7 +165,7 @@ simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, a
 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
+simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennell 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
@@ -223,26 +223,26 @@ cutoff sphere that are integral to the Wolf and DSF ap
 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.
+cutoff sphere that are integral to the Wolf and damped shifted force
+(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
 effect of long-range multipolar interactions.  The image multipoles
-that are introduced in the TSF, GSF, and SP methods mimic this effect
+that are introduced in the Taylor shifted force (TSF), gradient
+shifted force (GSF), and shifted potential (SP) methods mimic this effect
 and reduce the effective range of the multipolar interactions as
 interacting molecules traverse each other's cutoff boundaries.
 
 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 approaches 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
+dynamics in molecular simulations, and the 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
+approaches 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}
@@ -257,7 +257,7 @@ With moderate damping coefficients, $\alpha \sim 0.2$,
 produce complementary error functions when truncated at a finite
 distance.
 
-With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method
+With moderate damping coefficients, $\alpha \sim 0.2$ \AA$^{-1}$, the DSF method
 produced very good agreement with SPME for interaction energies,
 forces and torques for charge-charge
 interactions.\cite{Fennell:2006lq}
@@ -343,7 +343,7 @@ U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \
 \cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance
 dot products:
 \begin{equation}
-U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
+U_{\mathbf{D}_{a}\mathbf{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)
@@ -372,10 +372,10 @@ U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
 which have been projected onto the surface of the cutoff sphere
 without changing their relative orientation,
 \begin{equation}
-U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
-U_{D_{a}D_{b}}(r_c)
+U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r)  = U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r) -
+U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c)
    - (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot
-  \nabla U_{D_{a}D_{b}}(r_c).
+  \nabla U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c).
 \end{equation}
 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
@@ -482,11 +482,11 @@ in the test cases are given in table~\ref{tab:pars}. 
 in the test cases are given in table~\ref{tab:pars}. 
 
 \begin{table}
-\caption{The parameters used in the systems used to evaluate the new
-  real-space methods.  The most comprehensive test was a liquid
-  composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-}
-  ions).  This test exercises all orders of the multipolar
-  interactions developed in the first paper.\label{tab:pars}}
+  \caption{The parameters used in the systems used to evaluate the new
+    real-space methods.  The most comprehensive test was a liquid
+    composed of 2000 soft DQ liquid molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-}
+    ions).  This test exercises all orders of the multipolar
+    interactions developed in the first paper.\label{tab:pars}}
 \begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline
              & \multicolumn{2}{c|}{LJ parameters} &
              \multicolumn{5}{c|}{Electrostatic moments} & & & & \\
@@ -499,7 +499,7 @@ Soft Quadrupolar solid & 2.837 & 1.0   &  &  & -1&-1&-
     Soft Dipolar solid & 2.837 & 1.0   &  & 2.35 & & & & $10^4$  & 17.6 &17.6 & 0 \\
 Soft Quadrupolar fluid & 3.051 & 0.152 &  &  & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155  \\
 Soft Quadrupolar solid & 2.837 & 1.0   &  &  & -1&-1&-2.5 & $10^4$  & 17.6&17.6&0 \\
-      SSDQ water  & 3.051 & 0.152 &  & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\
+      Soft DQ liquid  & 3.051 & 0.152 &  & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\
               \ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\
               \ce{Cl-} & 4.445 & 0.1   & -1& & & & & 35.4527& & & \\ \hline
 \end{tabularx}
@@ -511,13 +511,15 @@ relatively strict translational order.  The SSDQ model
 charges in addition to the multipolar fluid.  The solid-phase
 parameters were chosen so that the systems can explore some
 orientational freedom for the multipolar sites, while maintaining
-relatively strict translational order.  The SSDQ model used here is
-not a particularly accurate water model, but it does test
-dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole
-interactions at roughly the same magnitudes. The last test case, SSDQ
-water with dissolved ions, exercises \textit{all} levels of the
-multipole-multipole interactions we have derived so far and represents
-the most complete test of the new methods.
+relatively strict translational order.  The soft DQ liquid model used
+here based loosely on the SSDQO water
+model,\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} but is not itself a
+particularly accurate water model.  However, the soft DQ model does
+test dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole
+interactions at roughly the same magnitudes. The last test case, a
+soft DQ liquid with dissolved ions, exercises \textit{all} levels of
+the multipole-multipole interactions we have derived so far and
+represents the most complete test of the new methods.
 
 In the following section, we present results for the total
 electrostatic energy, as well as the electrostatic contributions to
@@ -595,15 +597,18 @@ ensemble.  We collected 250 different configurations a
 simulations, each system was created with 2,048 randomly-oriented
 molecules.  These were equilibrated at a temperature of 300K for 1 ns.
 Each system was then simulated for 1 ns in the microcanonical (NVE)
-ensemble.  We collected 250 different configurations at equal time
-intervals. For the liquid system that included ionic species, we
-converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24
-\ce{Cl-} ions and re-equilibrated. After equilibration, the system was
-run under the same conditions for 1 ns. A total of 250 configurations
-were collected. In the following comparisons of energies, forces, and
-torques, the Lennard-Jones potentials were turned off and only the
-purely electrostatic quantities were compared with the same values
-obtained via the Ewald sum.
+ensemble with the Dullweber, Leimkuhler, and McLachlan (DLM)
+symplectic splitting integrator using 1 fs
+timesteps.\cite{Dullweber1997} We collected 250 different
+configurations at equal time intervals. For the liquid system that
+included ionic species, we converted 48 randomly-distributed molecules
+into 24 \ce{Na+} and 24 \ce{Cl-} ions and re-equilibrated. After
+equilibration, the system was run under the same conditions for 1
+ns. A total of 250 configurations were collected. In the following
+comparisons of energies, forces, and torques, the Lennard-Jones
+potentials were turned off and only the purely electrostatic
+quantities were compared with the same values obtained via the Ewald
+sum.
 
 \subsection{Accuracy of Energy Differences, Forces and Torques}
 The pairwise summation techniques (outlined above) were evaluated for
@@ -619,10 +624,10 @@ with the multipolar Ewald reference method.  Unitary r
 Since none of the real-space methods provide exact energy differences,
 we used least square regressions analysis for the six different
 molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF
-with the multipolar Ewald reference method.  Unitary results for both
-the correlation (slope) and correlation coefficient for these
-regressions indicate perfect agreement between the real-space method
-and the multipolar Ewald sum.
+with the multipolar Ewald reference method.  A result of unity for
+both the correlation (slope) and coefficient of determination ($R^2$)
+for these regressions would indicate perfect agreement between the
+real-space method and the multipolar Ewald sum.
 
 Molecular systems were run long enough to explore independent
 configurations and 250 configurations were recorded for comparison.
@@ -655,8 +660,8 @@ evaluated,
 the forces obtained via the Ewald sum and the real-space methods were
 evaluated,
 \begin{equation}
-\cos\theta_i =  \frac{\vec{f}_i^\mathrm{~Ewald} \cdot
-  \vec{f}_i^\mathrm{~GSF}}{\left|\vec{f}_i^\mathrm{~Ewald}\right| \left|\vec{f}_i^\mathrm{~GSF}\right|}
+  \cos\theta_i =  \frac{\mathbf{f}_i^\mathrm{~Ewald} \cdot
+    \mathbf{f}_i^\mathrm{~GSF}}{\left|\mathbf{f}_i^\mathrm{~Ewald}\right| \left|\mathbf{f}_i^\mathrm{~GSF}\right|}
 \end{equation}
 The total angular displacement of the vectors was calculated as,
 \begin{equation}
@@ -679,15 +684,15 @@ system of 2000 SSDQ water molecules with 24 \ce{Na+} a
 
 \subsection{Energy conservation}
 To test conservation the energy for the methods, the mixed molecular
-system of 2000 SSDQ water molecules with 24 \ce{Na+} and 24 \ce{Cl-}
-ions was run for 1 ns in the microcanonical ensemble at an average
-temperature of 300K.  Each of the different electrostatic methods
-(Ewald, Hard, SP, GSF, and TSF) was tested for a range of different
-damping values. The molecular system was started with same initial
-positions and velocities for all cutoff methods. The energy drift
-($\delta E_1$) and standard deviation of the energy about the slope
-($\delta E_0$) were evaluated from the total energy of the system as a
-function of time.  Although both measures are valuable at
+system of 2000 soft DQ liquid molecules with 24 \ce{Na+} and 24
+\ce{Cl-} ions was run for 1 ns in the microcanonical ensemble at an
+average temperature of 300K.  Each of the different electrostatic
+methods (Ewald, Hard, SP, GSF, and TSF) was tested for a range of
+different damping values. The molecular system was started with same
+initial positions and velocities for all cutoff methods. The energy
+drift ($\delta E_1$) and standard deviation of the energy about the
+slope ($\delta E_0$) were evaluated from the total energy of the
+system as a function of time.  Although both measures are valuable at
 investigating new methods for molecular dynamics, a useful interaction
 model must allow for long simulation times with minimal energy drift.
 
@@ -696,29 +701,29 @@ model must allow for long simulation times with minima
  
 \begin{figure}
   \centering
-  \includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps}
+  \includegraphics[width=0.6\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
     to 1 (dashed line) indicate $\Delta E$ values indistinguishable
     from those obtained using the multipolar Ewald sum.  Different
     values of the cutoff radius are indicated with different symbols
-    (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted
+    (9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = inverted
     triangles).\label{fig:slopeCorr_energy}}
 \end{figure} 
 
-The combined correlation coefficient and slope for all six systems is
-shown in Figure ~\ref{fig:slopeCorr_energy}.  Most of the methods
-reproduce the Ewald configurational energy differences with remarkable
-fidelity.  Undamped hard cutoffs introduce a significant amount of
-random scatter in the energy differences which is apparent in the
-reduced value of the correlation coefficient for this method.  This
-can be easily understood as configurations which exhibit small
-traversals of a few dipoles or quadrupoles out of the cutoff sphere
-will see large energy jumps when hard cutoffs are used.  The
-orientations of the multipoles (particularly in the ordered crystals)
-mean that these energy jumps can go in either direction, producing a
-significant amount of random scatter, but no systematic error.
+The combined coefficient of determination and slope for all six
+systems is shown in Figure ~\ref{fig:slopeCorr_energy}.  Most of the
+methods reproduce the Ewald configurational energy differences with
+remarkable fidelity.  Undamped hard cutoffs introduce a significant
+amount of random scatter in the energy differences which is apparent
+in the reduced value of $R^2$ for this method.  This can be easily
+understood as configurations which exhibit small traversals of a few
+dipoles or quadrupoles out of the cutoff sphere will see large energy
+jumps when hard cutoffs are used.  The orientations of the multipoles
+(particularly in the ordered crystals) mean that these energy jumps
+can go in either direction, producing a significant amount of random
+scatter, but no systematic error.
 
 The TSF method produces energy differences that are highly correlated
 with the Ewald results, but it also introduces a significant
@@ -729,9 +734,9 @@ excellent fidelity, particularly for moderate damping 
 effect, particularly for the crystalline systems.
 
 Both the SP and GSF methods appear to reproduce the Ewald results with
-excellent fidelity, particularly for moderate damping ($\alpha =
-0.1-0.2$\AA$^{-1}$) and with a commonly-used cutoff value ($r_c =
-12$\AA).  With the exception of the undamped hard cutoff, and the TSF
+excellent fidelity, particularly for moderate damping ($\alpha \approx
+0.2$~\AA$^{-1}$) and with a commonly-used cutoff value ($r_c =
+12$~\AA).  With the exception of the undamped hard cutoff, and the TSF
 method with short cutoffs, all of the methods would be appropriate for
 use in Monte Carlo simulations.
  
@@ -761,10 +766,9 @@ commonly-used cutoff values ($r_c = 12$\AA).  The TSF 
 energy conservation issues, and this perturbation is evident in the
 statistics accumulated for the molecular forces.  The GSF
 perturbations are minimal, particularly for moderate damping and
-commonly-used cutoff values ($r_c = 12$\AA).  The TSF method shows
-reasonable agreement in the correlation coefficient but again the
-systematic error in the forces is concerning if replication of Ewald
-forces is desired.
+commonly-used cutoff values ($r_c = 12$~\AA).  The TSF method shows
+reasonable agreement in $R^2$, but again the systematic error in the
+forces is concerning if replication of Ewald forces is desired.
 
 It is important to note that the forces and torques from the SP and
 the Hard cutoffs are not identical. The SP method shifts each
@@ -782,8 +786,8 @@ for multipoles even though the forces for point charge
     the reference Ewald sum. Results with a value equal to 1 (dashed
     line) indicate force magnitude values indistinguishable from those
     obtained using the multipolar Ewald sum.  Different values of the
-    cutoff radius are indicated with different symbols (9\AA\ =
-    circles, 12\AA\ = squares, and 15\AA\ = inverted
+    cutoff radius are indicated with different symbols (9~\AA\ =
+    circles, 12~\AA\ = squares, and 15~\AA\ = inverted
     triangles).\label{fig:slopeCorr_force}}
 \end{figure} 
 
@@ -796,8 +800,8 @@ for multipoles even though the forces for point charge
     the reference Ewald sum. Results with a value equal to 1 (dashed
     line) indicate force magnitude values indistinguishable from those
     obtained using the multipolar Ewald sum.  Different values of the
-    cutoff radius are indicated with different symbols (9\AA\ =
-    circles, 12\AA\ = squares, and 15\AA\ = inverted
+    cutoff radius are indicated with different symbols (9~\AA\ =
+    circles, 12~\AA\ = squares, and 15~\AA\ = inverted
     triangles).\label{fig:slopeCorr_torque}}
 \end{figure}
 
@@ -812,8 +816,8 @@ of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1
 reproduces the torques in quite good agreement with the Ewald sum.
 The other real-space methods can cause some deviations, but excellent
 agreement with the Ewald sum torques is recovered at moderate values
-of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1}$) and cutoff
-radius ($r_c \ge 12$\AA).  The TSF method exhibits only fair agreement
+of the damping coefficient ($\alpha \approx 0.2$~\AA$^{-1}$) and cutoff
+radius ($r_c \ge 12$~\AA).  The TSF method exhibits only fair agreement
 in the slope when compared with the Ewald torques even for larger
 cutoff radii.  It appears that the severity of the perturbations in
 the TSF method are most in evidence for the torques.
@@ -825,26 +829,27 @@ directionality is shown in terms of circular variance
 these quantities. Force and torque vectors for all six systems were
 analyzed using Fisher statistics, and the quality of the vector
 directionality is shown in terms of circular variance
-($\mathrm{Var}(\theta)$) in figure
-\ref{fig:slopeCorr_circularVariance}. The force and torque vectors
-from the new real-space methods exhibit nearly-ideal Fisher probability
-distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods
-exhibit the best vectorial agreement with the Ewald sum. The force and
-torque vectors from GSF method also show good agreement with the Ewald
-method, which can also be systematically improved by using moderate
-damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c =
-12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds
-to a distribution with 95\% of force vectors within $6.37^\circ$ of
-the corresponding Ewald forces. The TSF method produces the poorest
-agreement with the Ewald force directions.
+($\mathrm{Var}(\theta)$) in
+Fig. \ref{fig:slopeCorr_circularVariance}. The force and torque
+vectors from the new real-space methods exhibit nearly-ideal Fisher
+probability distributions (Eq.~\ref{eq:pdf}). Both the hard and SP
+cutoff methods exhibit the best vectorial agreement with the Ewald
+sum. The force and torque vectors from GSF method also show good
+agreement with the Ewald method, which can also be systematically
+improved by using moderate damping and a reasonable cutoff radius. For
+$\alpha = 0.2$~\AA$^{-1}$ and $r_c = 12$~\AA, we observe
+$\mathrm{Var}(\theta) = 0.00206$, which corresponds to a distribution
+with 95\% of force vectors within $6.37^\circ$ of the corresponding
+Ewald forces. The TSF method produces the poorest agreement with the
+Ewald force directions.
 
 Torques are again more perturbed than the forces by the new real-space
 methods, but even here the variance is reasonably small.  For the same
-method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA),
-the circular variance was 0.01415, corresponds to a distribution which
-has 95\% of torque vectors are within $16.75^\circ$ of the Ewald
-results. Again, the direction of the force and torque vectors can be
-systematically improved by varying $\alpha$ and $r_c$.
+method (GSF) with the same parameters ($\alpha = 0.2$~\AA$^{-1}$, $r_c
+= 12$~\AA), the circular variance was 0.01415, corresponds to a
+distribution which has 95\% of torque vectors are within $16.75^\circ$
+of the Ewald results. Again, the direction of the force and torque
+vectors can be systematically improved by varying $\alpha$ and $r_c$.
 
 \begin{figure}
   \centering
@@ -855,26 +860,26 @@ systematically improved by varying $\alpha$ and $r_c$.
     indicates direction of the force or torque vectors are
     indistinguishable from those obtained from the Ewald sum. Here
     different symbols represent different values of the cutoff radius
-    (9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}}
+    (9~\AA\ = circle, 12~\AA\ = square, 15~\AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}}
 \end{figure}
 
 \subsection{Energy conservation\label{sec:conservation}}
 
 We have tested the conservation of energy one can expect to see with
-the new real-space methods using the SSDQ water model with a small
+the new real-space methods using the soft DQ liquid model with a small
 fraction of solvated ions. This is a test system which exercises all
 orders of multipole-multipole interactions derived in the first paper
 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
+methods.  A liquid-phase system was created with 2000 liquid-phase
+molecules and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a
 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
+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.
+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
@@ -896,9 +901,9 @@ cutoff values are utilized.
 \begin{figure}
   \centering
   \includegraphics[width=\textwidth]{newDrift_12.eps}
-  \caption{Energy conservation of the real-space methods for the SSDQ
-    water/ion system. $\delta \mathrm{E}_1$ is the linear drift in
-    energy over time (in kcal/mol/particle/ns) and $\delta
+  \caption{Energy conservation of the real-space methods for the soft
+    DQ liauid / ion system. $\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 in
     the green region at the bottom exhibit better energy conservation
@@ -918,19 +923,20 @@ the different systems under investigation).  An exampl
 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}.
+the different systems under investigation). 
 
-\begin{figure}
-  \centering
-  \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
-\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.\label{fig:gofr}}
-\end{figure} 
+% An example of this agreement for the soft DQ liquid/ion system is
+% shown in Fig. \ref{fig:gofr}.
 
+% \begin{figure}
+%   \centering
+%   \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
+% \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.\label{fig:gofr}}
+% \end{figure} 
+
 There is a minor over-structuring of the first solvation shell when
 using TSF or when overdamping with any of the real-space methods.
 With moderate damping, GSF and SP produce pair distributions that are
@@ -942,7 +948,7 @@ in $g(r)$ ($a = 4.2$ \AA\  for the SSDQ water/ion syst
 \end{equation}
 where $\rho$ is the number density of the site-site pair interactions,
 and $a$ is the radial location of the minima following the first peak
-in $g(r)$ ($a = 4.2$ \AA\  for the SSDQ water/ion system).  The
+in $g(r)$ ($a = 4.2$~\AA\  for the soft DQ liquid / ion system).  The
 coordination number is shown as a function of the damping coefficient
 for all of the real space methods in Fig. \ref{fig:Props}.
 
@@ -950,24 +956,24 @@ of the methods.  In fig \ref{fig:Props} we demonstrate
 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 fig \ref{fig:Props} we demonstrate how $\langle
+of the methods.  In Fig. \ref{fig: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
+\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.28$ \AA$^{-1}$ for $r_c = 9$ \AA).
+dramatic damping ($\alpha \approx 0.28$~\AA$^{-1}$ for $r_c = 9$~\AA).
 Overdamping the real-space electrostatic methods occurs with $\alpha >
-0.3$, causing the estimate of the electrostatic energy to drop below
-the Ewald results.
+0.3$~\AA$^{-1}$, causing the estimate of the electrostatic energy to
+drop below the Ewald results.
 
 These ``optimal'' values of the damping coefficient for structural
 features are similar to those observed for DSF electrostatics for
 purely point-charge systems, and the range $\alpha= 0.175 \rightarrow
-0.225$ \AA$^{-1}$ for $r_c = 12$\AA\ appears to be an excellent
+0.225$~\AA$^{-1}$ for $r_c = 12$~\AA\ appears to be an excellent
 compromise for mixed charge/multipolar systems.
 
 To test the fidelity of the electrostatic methods at reproducing
@@ -981,7 +987,7 @@ $\langle r^{2}(t) \rangle$.\cite{Allen87} In fig. \ref
 which measures long-time behavior and is sensitive to the forces on
 the multipoles. 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 fig. \ref{fig:Props} we
+$\langle r^{2}(t) \rangle$.\cite{Allen87} In Fig. \ref{fig:Props} we
 demonstrate how the diffusion constant depends on the choice of
 real-space methods and the damping coefficient.  Both the SP and GSF
 methods can obtain excellent agreement with Ewald again using moderate
@@ -1000,14 +1006,14 @@ systems has the $z$-axis running along the dipoles, an
 translational diffusion.  Here, $P_l$ is the Legendre polynomial of
 order $l$ and $\hat{\mathbf{A}}_\gamma$ is the unit vector for body
 axis $\gamma$.  The reference frame used for our sample dipolar
-systems has the $z$-axis running along the dipoles, and for the SSDQ
-water model, the $y$-axis connects the two implied hydrogen atom
+systems has the $z$-axis running along the dipoles, and for the soft
+DQ liquid model, the $y$-axis connects the two implied hydrogen-like
 positions.  From the orientation autocorrelation functions, we can
 obtain time constants for rotational relaxation either by fitting to a
 multi-exponential model for the orientational relaxation, or by
 integrating the correlation functions.
 
-In a good water model, the orientational decay times would be
+In a good model for water, the orientational 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
@@ -1015,23 +1021,23 @@ In fig. \ref{fig:Props} we compare the $\tau_2^y$ and 
 atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
 of the dipole-dipole relaxation from a proton NMR signal and can
 provide an estimate of the NMR relaxation time constant.\cite{Impey82}
-In fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$
+In Fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$
 values for the various real-space methods over a range of different
 damping coefficients.  The rotational relaxation for the $z$ axis
 primarily probes the torques on the dipoles, while the relaxation for
 the $y$ axis is sensitive primarily to the quadrupolar torques.
 
 \begin{figure}
+  \includegraphics[width=\textwidth]{properties.eps}
   \caption{Comparison of the structural and dynamic properties for the
-    combined multipolar liquid (SSDQ water + ions) for all of the
-    real-space methods with $r_c = 12$\AA. Electrostatic energies,
+    combined multipolar liquid (soft DQ liquid + ions) for all of the
+    real-space methods with $r_c = 12$~\AA. Electrostatic energies,
     $\langle U_\mathrm{elect} \rangle / N$ (in kcal / mol),
     coordination numbers, $n_C$, diffusion constants (in $10^{-5}
     \mathrm{cm}^2\mathrm{s}^{-1}$), and rotational correlation times
     (in ps) all show excellent agreement with Ewald results for
     damping coefficients in the range $\alpha= 0.175 \rightarrow
-    0.225$ \AA$^{-1}$. \label{fig:Props}}
-  \includegraphics[width=\textwidth]{properties.eps}
+    0.225$~\AA$^{-1}$. \label{fig:Props}}
 \end{figure}
 
 In Fig. \ref{fig:Props} it appears that values for $D$, $\tau_2^y$,
@@ -1046,8 +1052,8 @@ parameter appears to be $\alpha= 0.175 \rightarrow 0.2
 the SP and GSF methods, as their structural and dynamical properties
 still reproduce the Ewald results even in the completely undamped
 ($\alpha = 0$) case.  An optimal range for the electrostatic damping
-parameter appears to be $\alpha= 0.175 \rightarrow 0.225$ \AA$^{-1}$
-for $r_c = 12$\AA, which similar to the optimal range found for the
+parameter appears to be $\alpha= 0.175 \rightarrow 0.225$~\AA$^{-1}$
+for $r_c = 12$~\AA, which similar to the optimal range found for the
 damped shifted force potential for point charges.\cite{Fennell:2006lq}
 
 \section{CONCLUSION}
@@ -1063,7 +1069,7 @@ Fennel and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} 
 
 We also developed two natural extensions of the damped shifted-force
 (DSF) model originally proposed by Zahn {\it et al.} and extended by
-Fennel and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF
+Fennell and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF
 approaches provide smooth truncation of energies, forces, and torques
 at the real-space cutoff, and both converge to DSF electrostatics for
 point-charge interactions.  The TSF model is based on a high-order
@@ -1085,8 +1091,8 @@ useful in past treatments of the corrections to the fl
 series of the electrostatic kernel at the cutoff radius. The TSF
 method also has the unique property that a large number of derivatives
 can be made to vanish at the cutoff radius.  This property has proven
-useful in past treatments of the corrections to the fluctuation
-formula for dielectric constants.\cite{Izvekov:2008wo}
+useful in past treatments of the corrections to the Clausius-Mossotti
+fluctuation formula for dielectric constants.\cite{Izvekov:2008wo}
 
 Reproduction of both structural and dynamical features in the liquid
 systems is remarkably good for both the SP and GSF models.  Pair
