--- trunk/multipole/multipole_2/multipole2.tex 2014/06/12 14:58:06 4181 +++ trunk/multipole/multipole_2/multipole2.tex 2014/06/14 23:35:39 4184 @@ -70,27 +70,27 @@ of Notre Dame, Notre Dame, IN 46556} 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. + against 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 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. \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 @@ -122,38 +122,39 @@ periodicity in the Ewald’s method can also be proble 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} +periodicity in the Ewald method can also be problematic for molecular +interfaces.\cite{Fennell:2006lq} 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. \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 +condensed phase 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} +application of Wolf's method can obtain accurate estimates of Madelung +constants using relatively short cutoff radii. Recently, Fukuda used +neutralization of the higher order moments for calculation of the +electrostatic interactions in 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 @@ -163,7 +164,7 @@ truncation defects. Wolf \textit{et al.} further argue \end{figure} The direct truncation of interactions at a cutoff radius creates -truncation defects. Wolf \textit{et al.} further argued that +truncation defects. 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 @@ -339,7 +340,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 +413,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 +431,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 +455,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 +477,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 +491,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