--- trunk/multipole/multipole_2/multipole2.tex	2014/08/06 19:10:04	4203
+++ trunk/multipole/multipole_2/multipole2.tex	2014/08/07 21:13:00	4207
@@ -220,7 +220,7 @@ negaitve ions and vice versa.  The reversed-charge ima
 limited to perfect crystals, but can also apply in disordered fluids.
 Even at elevated temperatures, there is 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
+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.
@@ -325,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
@@ -385,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,
@@ -418,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
@@ -445,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
@@ -480,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
@@ -526,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}
@@ -596,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
@@ -958,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
@@ -988,69 +990,73 @@ $k$-space cutoff values.
   \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 Features\label{sec:structure}}
-One of the best tests of modified interaction potentials is the
-fidelity with which they can reproduce structural features 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 any 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).  Interested readers may consult the
-supplementary information for plots of these pair distribution
-functions.
+\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}.
 
-A direct measure of the structural features that is a more
-enlightening test of the modified electrostatic methods is the average
-value of the electrostatic energy $\langle U_\mathrm{elect} \rangle$
-which is obtained by sampling the liquid-state configurations
-experienced by a liquid evolving entirely under the influence of the
-methods being investigated.  In figure \ref{fig:Uelect} we show how
-$\langle U_\mathrm{elect} \rangle$ for varies with the damping parameter,
-$\alpha$, for each of the methods.
-
 \begin{figure}
   \centering
-  \includegraphics[width=\textwidth]{averagePotentialEnergy_r9_12.eps}
-\label{fig:Uelect}        
-\caption{The average electrostatic potential energy, 
-  $\langle U_\mathrm{elect} \rangle$ for the SSDQ water with ions as a function
-  of the damping parameter, $\alpha$, for each of the real-space
-  electrostatic methods. Top panel: simulations run with a real-space
-  cutoff, $r_c = 9$\AA.  Bottom panel: the same quantity, but with a
-  larger cutoff, $r_c = 12$\AA.}
+  \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} 
 
-It is clear that moderate damping is important for converging the mean
-potential energy values, particularly for the two shifted force
-methods (GSF and TSF).  A value of $\alpha \approx 0.18$ \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.36$ \AA$^{-1}$ for $r_c = 9$ \AA).  It is also clear from
-fig. \ref{fig:Uelect} that it is possible to overdamp the real-space
-electrostatic methods, causing the estimate of the energy to drop
-below the Ewald results.
+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.
 
-\subsection{Reproduction of Dynamic Properties\label{sec:structure}}
 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,
@@ -1060,49 +1066,90 @@ the multipoles.  For the soft dipolar fluid, and the S
 \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
+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}
+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 $NVE$
-trajectories used for translational diffusion by calculating the
-orientational time correlation function,
+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^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t)
-                \cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle,
+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{u}}_i^\alpha$ is the unit vector of molecule $i$ along
-axis $\gamma$.  The body-fixed reference frame used for our
-orientational correlation functions has the $z$-axis running along the
-dipoles, and for the SSDQ water model, the $y$-axis connects the two
-implied hydrogen atoms.
+$\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}
 
-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.  These
-decay times are directly 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{fig:dynamics}. From this data, it is apparent
+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 high fidelity by the GSF method. 
+reproduced with reasonable fidelity by the GSF method.
 
-The $\tau_2$ results in \ref{fig:dynamics} show a much greater
-difference between the real-space and the Ewald results.
+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
