--- trunk/multipole/multipole_2/multipole2.tex 2014/06/16 12:49:00 4190 +++ trunk/multipole/multipole_2/multipole2.tex 2014/08/07 20:53:27 4206 @@ -47,7 +47,8 @@ preprint, %\preprint{AIP/123-QED} -\title{Real space alternatives to the Ewald Sum. II. Comparison of Methods} +\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} @@ -117,15 +118,15 @@ interactions,\cite{Parry:1975if,Parry:1976fq,Clarke77, 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} -but these methods were originally quite computationally +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,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} +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. The inherent periodicity in the Ewald method can also -be problematic for interfacial molecular -systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} +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)$ @@ -324,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 @@ -384,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, @@ -417,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 @@ -444,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 @@ -479,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 @@ -573,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 @@ -792,7 +791,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 @@ -865,7 +864,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 @@ -879,7 +878,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 @@ -937,7 +936,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) @@ -957,13 +956,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 @@ -987,16 +987,166 @@ $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 \& 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