--- trunk/multipole/multipole_2/multipole2.tex 2014/06/12 14:58:06 4181 +++ trunk/multipole/multipole_2/multipole2.tex 2014/06/15 00:18:18 4185 @@ -69,43 +69,43 @@ of Notre Dame, Notre Dame, IN 46556} \begin{abstract} We have tested the real-space shifted potential (SP), gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods - for multipoles that were developed in the first paper in this series - against 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 sum. We have - also investigated the energy conservation properties of the new - methods in molecular dynamics simulations using all of these - methods. 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 - efficiently computing electrostatic interactions in molecular - dynamics simulations. + 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. \end{abstract} %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy % Classification Scheme. -\keywords{Electrostatics, Multipoles, Real-space} +%\keywords{Electrostatics, Multipoles, Real-space} \maketitle \section{\label{sec:intro}Introduction} Computing the interactions between electrostatic sites is one of the -most expensive aspects of molecular simulations, which is why there -have been significant efforts to develop practical, efficient and -convergent methods for handling these interactions. Ewald's method is -perhaps the best known and most accurate method for evaluating -energies, forces, and torques in explicitly-periodic simulation -cells. In this approach, 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.\cite{Clarke:1986eu,Woodcock75} +most expensive aspects of molecular simulations. There have been +significant efforts to develop practical, efficient and convergent +methods for handling these interactions. Ewald's method is perhaps the +best known and most accurate method for evaluating energies, forces, +and torques in explicitly-periodic simulation cells. In this approach, +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} When carried out as originally formulated, the reciprocal-space portion of the Ewald sum exhibits relatively poor computational @@ -116,44 +116,74 @@ Because of the artificial periodicity required for the 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, the -method may require modification to compute interactions for +Because of the artificial periodicity required for the Ewald sum, interfacial molecular systems such as membranes and liquid-vapor -interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} -To simulate interfacial systems, Parry's extension of the 3D Ewald sum -is appropriate for slab geometries.\cite{Parry:1975if} The inherent -periodicity in the Ewald’s method can also be problematic for -interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald -methods that were developed to handle two-dimensional (2D) -electrostatic interactions in interfacial systems have not had similar -particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, - Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} +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} \subsection{Real-space methods} Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ method for calculating electrostatic interactions between point -charges. They argued that the effective Coulomb interaction in -condensed systems is actually short ranged.\cite{Wolf92,Wolf95} For an -ordered lattice (e.g., when computing the Madelung constant of an -ionic solid), the material can be considered as a set of ions -interacting with neutral dipolar or quadrupolar ``molecules'' giving -an effective distance dependence for the electrostatic interactions of -$r^{-5}$ (see figure \ref{fig:schematic}). For this reason, careful -applications of Wolf's method are able to obtain accurate estimates of -Madelung constants using relatively short cutoff radii. Recently, -Fukuda used neutralization of the higher order moments for the -calculation of the electrostatic interaction of the point charges -system.\cite{Fukuda:2013sf} +charges. They argued that the effective Coulomb interaction in most +condensed phase systems is effectively short +ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when +computing the Madelung constant of an ionic solid), the material can +be considered as a set of ions interacting with neutral dipolar or +quadrupolar ``molecules'' giving an effective distance dependence for +the electrostatic interactions of $r^{-5}$ (see figure +\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple +cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the +electrostatic energy per ion converges more rapidly to the Madelung +energy than the dipolar approximation.\cite{Wolf92} To find the +correct Madelung constant, Lacman suggested that the NaCl structure +could be constructed in a way that the finite crystal terminates with +complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees +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 +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 +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 +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 +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 + 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} + \begin{figure} \centering \includegraphics[width=\linewidth]{schematic.pdf} \caption{Top: Ionic systems exhibit local clustering of dissimilar charges (in the smaller grey circle), so interactions are - effectively charge-multipole in order at longer distances. With - hard cutoffs, motion of individual charges in and out of the - cutoff sphere can break the effective multipolar ordering. - Bottom: dipolar crystals and fluids have a similar effective + effectively charge-multipole at longer distances. With hard + cutoffs, motion of individual charges in and out of the cutoff + sphere can break the effective multipolar ordering. Bottom: + dipolar crystals and fluids have a similar effective \textit{quadrupolar} ordering (in the smaller grey circles), and orientational averaging helps to reduce the effective range of the interactions in the fluid. Placement of reversed image multipoles @@ -162,46 +192,16 @@ The direct truncation of interactions at a cutoff radi \label{fig:schematic} \end{figure} -The direct truncation of interactions at a cutoff radius creates -truncation defects. Wolf \textit{et al.} further argued 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 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 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 - al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} -and by Wu and Brooks.\cite{Wu:044107} - -Considering the interaction of one central ion in an ionic crystal -with a portion of the crystal at some distance, the effective Columbic -potential is found to decrease as $r^{-5}$. If one views the \ce{NaCl} -crystal as a simple cubic (SC) structure with an octupolar -\ce{(NaCl)4} basis, the electrostatic energy per ion converges more -rapidly to the Madelung energy than the dipolar -approximation.\cite{Wolf92} To find the correct Madelung constant, -Lacman suggested that the NaCl structure could be constructed in a way -that the finite crystal terminates with complete \ce{(NaCl)4} -molecules.\cite{Lacman65} The central ion sees 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} - -One can make a similar argument for crystals of point 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 moments, while the rest have -zero net dipole and retain contributions only from higher order -multipoles. The effective interaction between a dipole at the center +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 +moments, while the rest have zero net dipole and retain contributions +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 of a crystal and a group of eight dipoles farther away is significantly shorter ranged than the $r^{-3}$ that one would expect for raw dipole-dipole interactions. Only in crystals which retain a @@ -258,48 +258,42 @@ The damping function used in our research has been dis densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} \subsection{The damping function} -The damping function used in our research has been discussed in detail -in the first paper of this series.\cite{PaperI} The radial kernel -$1/r$ for the interactions between point charges can be replaced by -the complementary error function $\mathrm{erfc}(\alpha r)/r$ to -accelerate the rate of convergence, where $\alpha$ is a damping -parameter with units of inverse distance. Altering the value of -$\alpha$ is equivalent to changing the width of Gaussian charge -distributions that replace each point charge -- Gaussian overlap -integrals yield complementary error functions when truncated at a -finite distance. +The damping function has been discussed in detail in the first paper +of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the +interactions between point charges can be replaced by the +complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate +convergence, where $\alpha$ is a damping parameter with units of +inverse distance. Altering the value of $\alpha$ is equivalent to +changing the width of Gaussian charge distributions that replace each +point charge, as Coulomb integrals with Gaussian charge distributions +produce complementary error functions when truncated at a finite +distance. -By using suitable value of damping alpha ($\alpha \sim 0.2$) for a -cutoff radius ($r_{­c}=9 A$), Fennel and Gezelter produced very good -agreement with SPME for the interaction energies, forces and torques -for charge-charge interactions.\cite{Fennell:2006lq} +With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method +produced very good agreement with SPME for interaction energies, +forces and torques for charge-charge +interactions.\cite{Fennell:2006lq} \subsection{Point multipoles in molecular modeling} Coarse-graining approaches which treat entire molecular subsystems as a single rigid body are now widely used. A common feature of many coarse-graining approaches is simplification of the electrostatic interactions between bodies so that fewer site-site interactions are -required to compute configurational energies. Many coarse-grained -molecular structures would normally consist of equal positive and -negative charges, and rather than use multiple site-site interactions, -the interaction between higher order multipoles can also be used to -evaluate a single molecule-molecule -interaction.\cite{Ren06,Essex10,Essex11} +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 a relatively -rough 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) +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) 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} But -using point multipole with the real space truncation without -accounting for multipolar neutrality will create energy conservation -issues in molecular dynamics (MD) simulations. +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. In this paper we test a set of real-space methods that were developed for point multipolar interactions. These methods extend the damped @@ -339,7 +333,7 @@ $\bf a$. 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$. +$\bf a$, etc. % Interactions between multipoles can be expressed as higher derivatives % of the bare Coulomb potential, so one way of ensuring that the forces @@ -412,7 +406,7 @@ to another site within cutoff sphere are derived from connection to unmodified electrostatics as well as the smooth transition to zero in both these functions as $r\rightarrow r_c$. The electrostatic forces and torques acting on the central multipole due -to another site within cutoff sphere are derived from +to another site within the cutoff sphere are derived from Eq.~\ref{generic}, accounting for the appropriate number of derivatives. Complete energy, force, and torque expressions are presented in the first paper in this series (Reference @@ -430,7 +424,7 @@ U_{D_{\bf a} D_{\bf b}}(r_c) 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 - \vec{\nabla} U_{D_{\bf a}D_{\bf b}}(r) \Big \lvert _{r_c} + \nabla U_{D_{\bf a}D_{\bf 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 @@ -454,18 +448,21 @@ U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathb 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(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} -\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] + 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] \label{generic2} \end{equation} where the sum describes a separate force-shifting that is applied to -each orientational contribution to the energy. +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}}$ +represent the orientations the multipoles. The third term converges more rapidly than the first two terms as a function of radius, hence the contribution of the third term is very small for large cutoff radii. The force and torque derived from -equation \ref{generic2} are consistent with the energy expression and +Eq. \ref{generic2} are consistent with the energy expression and approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods can be considered generalizations of the original DSF method for higher order multipole interactions. GSF and TSF are also identical up @@ -473,7 +470,7 @@ GSF potential are presented in the first paper in this the energy, force and torque for higher order multipole-multipole interactions. Complete energy, force, and torque expressions for the GSF potential are presented in the first paper in this series -(Reference~\onlinecite{PaperI}) +(Reference~\onlinecite{PaperI}). \subsection{Shifted potential (SP) } @@ -487,7 +484,7 @@ U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \ri 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(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] +U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] \label{eq:SP} \end{equation} where the sum describes separate potential shifting that is done for