--- trunk/multipole/multipole_2/multipole2.tex 2014/06/15 00:18:18 4185 +++ trunk/multipole/multipole_2/multipole2.tex 2014/08/06 19:10:04 4203 @@ -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 multipole interactions that were 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 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. 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 efficient computation of electrostatic interactions in - molecular dynamics simulations. + 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 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,7 +90,6 @@ 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. There have been @@ -105,30 +100,33 @@ space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} 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. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} +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, interfacial molecular systems such as membranes and liquid-vapor -interfaces require modifications to the -method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} -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 (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. The inherent periodicity in the Ewald’s method can -also be problematic for interfacial molecular -systems.\cite{Fennell:2006lq} +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)$ @@ -150,34 +148,33 @@ application of Wolf's method are able to obtain accura 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 are able to obtain accurate estimates of +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.} argued that truncation errors are due +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 -maintained a discontinuity at the cutoff sphere. In the original Wolf +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 their force expression, and as a result, the total +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{Fennell:2006lq} Related real-space methods were also -proposed by Chen \textit{et +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 used -neutralization of the higher order moments for the calculation of the -electrostatic interaction of the point charge -systems.\cite{Fukuda:2013sf} +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 @@ -195,10 +192,10 @@ dipoles at the eight corners of a unit cube. Only thr 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. Only three of these -basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole +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 +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 @@ -221,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 +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. In multipolar fluids (see Fig. \ref{fig:schematic}) there is significant orientational averaging that additionally reduces the @@ -245,17 +242,16 @@ 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 has been discussed in detail in the first paper @@ -282,29 +278,30 @@ Because electrons in a molecule are not localized at s 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 always an -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} However, truncating point multipoles without smoothing the forces and torques -will create energy conservation issues in molecular dynamics simulations. +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 } @@ -577,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 @@ -589,7 +586,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 @@ -796,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 @@ -869,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 @@ -883,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 @@ -941,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) @@ -988,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 @@ -1000,7 +997,112 @@ $k$-space cutoff values. 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