--- trunk/multipole/multipole_2/multipole2.tex 2014/06/15 16:27:42 4188 +++ trunk/multipole/multipole_2/multipole2.tex 2014/08/06 19:10:04 4203 @@ -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)$ @@ -173,7 +174,7 @@ 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 @@ -573,7 +574,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 +793,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 +866,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 +880,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 +938,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) @@ -984,7 +985,7 @@ $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 @@ -994,9 +995,114 @@ $k$-space cutoff values. 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. + +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.} \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. +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, +\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 $NVE$ +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, +\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. + +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 +that the values for both $D$ and $\tau_2$ using the Ewald sum are +reproduced with high 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. + + \section{CONCLUSION} In the first paper in this series, we generalized the charge-neutralized electrostatic energy originally developed by Wolf