--- trunk/multipole/multipole1.tex 2014/01/02 15:46:58 3988 +++ trunk/multipole/multipole1.tex 2014/06/05 17:09:34 4173 @@ -152,22 +152,23 @@ V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i An efficient real-space electrostatic method involves the use of a pair-wise functional form, \begin{equation} -V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) + -\sum_i V_i^\mathrm{correction} +V = \sum_i \sum_{j>i} V_\mathrm{pair}(\mathbf{r}_{ij}, \Omega_i, \Omega_j) + +\sum_i V_i^\mathrm{self} \end{equation} that is short-ranged and easily truncated at a cutoff radius, \begin{equation} - V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) = \left\{ + V_\mathrm{pair}(\mathbf{r}_{ij},\Omega_i, \Omega_j) = \left\{ \begin{array}{ll} -V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\ -0 & \quad r > r_c , +V_\mathrm{approx} (\mathbf{r}_{ij}, \Omega_i, \Omega_j) & \quad \left| \mathbf{r}_{ij} \right| \le r_c \\ +0 & \quad \left| \mathbf{r}_{ij} \right| > r_c , \end{array} \right. \end{equation} -along with an easily computed correction term ($\sum_i -V_i^\mathrm{correction}$) which has linear-scaling with the number of +along with an easily computed self-interaction term ($\sum_i +V_i^\mathrm{self}$) which scales linearly with the number of particles. Here $\Omega_i$ and $\Omega_j$ represent orientational -coordinates of the two sites. The computational efficiency, energy +coordinates of the two sites, and $\mathbf{r}_{ij}$ is the vector +between the two sites. The computational efficiency, energy conservation, and even some physical properties of a simulation can depend dramatically on how the $V_\mathrm{approx}$ function behaves at the cutoff radius. The goal of any approximation method should be to @@ -180,10 +181,9 @@ computed within $r_c$. Damping using a complementary e contained within the cutoff sphere surrounding each particle. This is accomplished by shifting the potential functions to generate image charges on the surface of the cutoff sphere for each pair interaction -computed within $r_c$. Damping using a complementary error -function is applied to the potential to accelerate convergence. The -potential for the DSF method, where $\alpha$ is the adjustable damping -parameter, is given by +computed within $r_c$. Damping using a complementary error function is +applied to the potential to accelerate convergence. The interaction +for a pair of charges ($C_i$ and $C_j$) in the DSF method, \begin{equation*} V_\mathrm{DSF}(r) = C_i C_j \Biggr{[} \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} @@ -193,9 +193,9 @@ Note that in this potential and in all electrostatic q \right)\left(r_{ij}-r_c\right)\ \Biggr{]} \label{eq:DSFPot} \end{equation*} -Note that in this potential and in all electrostatic quantities that -follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for -clarity. +where $\alpha$ is the adjustable damping parameter. Note that in this +potential and in all electrostatic quantities that follow, the +standard $1/4 \pi \epsilon_{0}$ has been omitted for clarity. To insure net charge neutrality within each cutoff sphere, an additional ``self'' term is added to the potential. This term is @@ -250,7 +250,7 @@ labelling specific charges in $\bf a$ and $\bf b$ resp The Taylor expansion in $r$ can be written using an implied summation notation. Here Greek indices are used to indicate space coordinates ($x$, $y$, $z$) and the subscripts $k$ and $j$ are reserved for -labelling specific charges in $\bf a$ and $\bf b$ respectively. The +labeling specific charges in $\bf a$ and $\bf b$ respectively. The Taylor expansion, \begin{equation} \frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = @@ -276,13 +276,20 @@ C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\ a}\alpha\beta}$, respectively. These are the primitive multipoles which can be expressed as a distribution of charges, \begin{align} -C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\ -D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha} ,\\ -Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k r_{k\alpha} r_{k\beta} . +C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \label{eq:charge} \\ +D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha}, \label{eq:dipole}\\ +Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k +r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} \end{align} Note that the definition of the primitive quadrupole here differs from the standard traceless form, and contains an additional Taylor-series -based factor of $1/2$. +based factor of $1/2$. We are essentially treating the mass +distribution with higher priority; the moment of inertia tensor, +$\overleftrightarrow{\mathsf I}$, is diagonalized to obtain body-fixed +axes, and the charge distribution may result in a quadrupole tensor +that is not necessarily diagonal in the body frame. Additional +reasons for utilizing the primitive quadrupole are discussed in +section \ref{sec:damped}. It is convenient to locate charges $q_j$ relative to the center of mass of $\bf b$. Then with $\bf{r}$ pointing from $\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by @@ -306,6 +313,7 @@ In the standard multipole expansion, one typically use of $\bf a$ interacting with the same multipoles on $\bf b$. \subsection{Damped Coulomb interactions} +\label{sec:damped} In the standard multipole expansion, one typically uses the bare Coulomb potential, with radial dependence $1/r$, as shown in Eq.~(\ref{kernel}). It is also quite common to use a damped Coulomb @@ -321,7 +329,13 @@ functions $B_l(r)$ are summarized in Appendix A. either $1/r$ or $B_0(r)$, and all of the techniques can be applied to bare or damped Coulomb kernels (or any other function) as long as derivatives of these functions are known. Smith's convenient -functions $B_l(r)$ are summarized in Appendix A. +functions $B_l(r)$ are summarized in Appendix A. (N.B. there is one +important distinction between the two kernels, which is the behavior +of $\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$. The former +is zero everywhere except for a delta function evaluated at the +origin. The latter also has delta function behavior, but is non-zero +for $r \neq 0$. Thus the standard justification for using a traceless +quadrupole tensor fails for the damped case.) The main goal of this work is to smoothly cut off the interaction energy as well as forces and torques as $r\rightarrow r_c$. To @@ -418,7 +432,7 @@ U= (\text{prefactor}) (\text{derivatives}) f_n(r) In general, we can write % \begin{equation} -U= (\text{prefactor}) (\text{derivatives}) f_n(r) +U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) \label{generic} \end{equation} % @@ -439,6 +453,11 @@ presented in tables \ref{tab:tableenergy} and \ref{tab of the same index $n$. The algebra required to evaluate energies, forces and torques is somewhat tedious, so only the final forms are presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}. +One of the principal findings of our work is that the individual +orientational contributions to the various multipole-multipole +interactions must be treated with distinct radial functions, but each +of these contributions is independently force shifted at the cutoff +radius. \subsection{Gradient-shifted force (GSF) electrostatics} The second, and conceptually simpler approach to force-shifting @@ -446,19 +465,23 @@ U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \n expansion, and has a similar interaction energy for all multipole orders: \begin{equation} -U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big -\lvert _{r_c} . +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] \label{generic2} \end{equation} -Here the gradient for force shifting is evaluated for an image -multipole projected onto the surface of the cutoff sphere (see fig -\ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$ -appear. The primary difference between the TSF and GSF methods is the -stage at which the Taylor Series is applied; in the Taylor-shifted -approach, it is applied to the kernel itself. In the Gradient-shifted -approach, it is applied to individual radial interactions terms in the -multipole expansion. Energies from this method thus have the general -form: +where the sum describes a separate force-shifting that is applied to +each orientational contribution to the energy. Both the potential and +the gradient for force shifting are evaluated for an image multipole +projected onto the surface of the cutoff sphere (see fig +\ref{fig:shiftedMultipoles}). The image multipole retains the +orientation ($\hat{\mathbf{b}}$) of the interacting multipole. No +higher order terms $(r-r_c)^n$ appear. The primary difference between +the TSF and GSF methods is the stage at which the Taylor Series is +applied; in the Taylor-shifted approach, it is applied to the kernel +itself. In the Gradient-shifted approach, it is applied to individual +radial interactions terms in the multipole expansion. Energies from +this method thus have the general form: \begin{equation} U= \sum (\text{angular factor}) (\text{radial factor}). \label{generic3} @@ -471,32 +494,20 @@ the radial factors differ between the two methods. the radial factors differ between the two methods. \subsection{\label{sec:level2}Body and space axes} +Although objects $\bf a$ and $\bf b$ rotate during a molecular +dynamics (MD) simulation, their multipole tensors remain fixed in +body-frame coordinates. While deriving force and torque expressions, +it is therefore convenient to write the energies, forces, and torques +in intermediate forms involving the vectors of the rotation matrices. +We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors +$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$. +In a typical simulation , the initial axes are obtained by +diagonalizing the moment of inertia tensors for the objects. (N.B., +the body axes are generally {\it not} the same as those for which the +quadrupole moment is diagonal.) The rotation matrices are then +propagated during the simulation. -[XXX Do we need this section in the main paper? or should it go in the -extra materials?] - -So far, all energies and forces have been written in terms of fixed -space coordinates. Interaction energies are computed from the generic -formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine -orientational prefactors with radial functions. Because objects $\bf -a$ and $\bf b$ both translate and rotate during a molecular dynamics -(MD) simulation, it is desirable to contract all $r$-dependent terms -with dipole and quadrupole moments expressed in terms of their body -axes. To do so, we have followed the methodology of Allen and -Germano,\cite{Allen:2006fk} which was itself based on earlier work by -Price {\em et al.}\cite{Price:1984fk} - -We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors -$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$ -referring to a convenient set of inertial body axes. (N.B., these -body axes are generally not the same as those for which the quadrupole -moment is diagonal.) Then, -% -\begin{eqnarray} -\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z} \\ -\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z} . -\end{eqnarray} -Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be +The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be expressed using these unit vectors: \begin{eqnarray} \hat{\mathbf {a}} = @@ -504,35 +515,82 @@ expressed using these unit vectors: \hat{a}_1 \\ \hat{a}_2 \\ \hat{a}_3 -\end{pmatrix} -= -\begin{pmatrix} -a_{1x} \quad a_{1y} \quad a_{1z} \\ -a_{2x} \quad a_{2y} \quad a_{2z} \\ -a_{3x} \quad a_{3y} \quad a_{3z} -\end{pmatrix}\\ +\end{pmatrix}, \qquad \hat{\mathbf {b}} = \begin{pmatrix} \hat{b}_1 \\ \hat{b}_2 \\ \hat{b}_3 \end{pmatrix} -= -\begin{pmatrix} -b_{1x} \quad b_{1y} \quad b_{1z} \\ -b_{2x} \quad b_{2y} \quad b_{2z} \\ -b_{3x} \quad b_{3y} \quad b_{3z} -\end{pmatrix} . -\end{eqnarray} +\end{eqnarray} % These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$ -coordinates. All contractions of prefactors with derivatives of -functions can be written in terms of these matrices. It proves to be -equally convenient to just write any contraction in terms of unit -vectors $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$. In the torque -expressions, it is useful to have the angular-dependent terms -available in three different fashions, e.g. for the dipole-dipole -contraction: +coordinates. + +Allen and Germano,\cite{Allen:2006fk} following earlier work by Price +{\em et al.},\cite{Price:1984fk} showed that if the interaction +energies are written explicitly in terms of $\hat{r}$ and the body +axes ($\hat{a}_m$, $\hat{b}_n$) : +% +\begin{equation} +U(r, \{\hat{a}_m \cdot \hat{r} \}, +\{\hat{b}_n\cdot \hat{r} \}, +\{\hat{a}_m \cdot \hat{b}_n \}) . +\label{ugeneral} +\end{equation} +% +the forces come out relatively cleanly, +% +\begin{equation} +\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} += \frac{\partial U}{\partial r} \hat{r} + + \sum_m \left[ +\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} +\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} ++ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} +\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} +\right] \label{forceequation}. +\end{equation} + +The torques can also be found in a relatively similar +manner, +% +\begin{eqnarray} +\mathbf{\tau}_{\bf a} = + \sum_m +\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} +( \hat{r} \times \hat{a}_m ) +-\sum_{mn} +\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} +(\hat{a}_m \times \hat{b}_n) \\ +% +\mathbf{\tau}_{\bf b} = + \sum_m +\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} +( \hat{r} \times \hat{b}_m) ++\sum_{mn} +\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} +(\hat{a}_m \times \hat{b}_n) . +\end{eqnarray} + +Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ +is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} +We also made use of the identities, +% +\begin{align} +\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} +=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} +\right) \\ +\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} +=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} +\right) . +\end{align} + +Many of the multipole contractions required can be written in one of +three equivalent forms using the unit vectors $\hat{r}$, $\hat{a}_m$, +and $\hat{b}_n$. In the torque expressions, it is useful to have the +angular-dependent terms available in all three fashions, e.g. for the +dipole-dipole contraction: % \begin{equation} \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} @@ -546,12 +604,17 @@ contractions using space indices. explicit sums over body indices and the dot products now indicate contractions using space indices. +In computing our force and torque expressions, we carried out most of +the work in body coordinates, and have transformed the expressions +back to space-frame coordinates, which are reported below. Interested +readers may consult the supplemental information for this paper for +the intermediate body-frame expressions. \subsection{The Self-Interaction \label{sec:selfTerm}} In addition to cutoff-sphere neutralization, the Wolf summation~\cite{Wolf99} and the damped shifted force (DSF) -extension~\cite{Fennell:2006zl} also included self-interactions that +extension~\cite{Fennell:2006zl} also include self-interactions that are handled separately from the pairwise interactions between sites. The self-term is normally calculated via a single loop over all sites in the system, and is relatively cheap to evaluate. The @@ -598,7 +661,7 @@ site which posesses a charge, dipole, and multipole, b reciprocal-space portion is identical to half of the self-term obtained by Smith and Aguado and Madden for the application of the Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given -site which posesses a charge, dipole, and multipole, both types of +site which posesses a charge, dipole, and quadrupole, both types of contribution are given in table \ref{tab:tableSelf}. \begin{table*} @@ -611,7 +674,7 @@ Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \tex Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\ Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) + \frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\ -Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & +Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & $- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ & $-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\ Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left( @@ -741,8 +804,8 @@ +2 \text{Tr} \left( U_{Q_{\bf a}Q_{\bf b}}(r)=& \Bigl[ \text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}} -+2 \text{Tr} \left( -\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r) ++2 +\mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r) \\ % 2 &+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}} @@ -981,49 +1044,8 @@ F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_ F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r} \quad \text{and} \quad F_{\bf b \alpha} = - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r} . -\end{equation} -% -Obtaining the force from the interaction energy expressions is the -same for higher-order multipole interactions -- the trick is to make -sure that all $r$-dependent derivatives are considered. This is -straighforward if the interaction energies are written explicitly in -terms of $\hat{r}$ and the body axes ($\hat{a}_m$, -$\hat{b}_n$) : -% -\begin{equation} -U(r,\{\hat{a}_m \cdot \hat{r} \}, -\{\hat{b}_n\cdot \hat{r} \}, -\{\hat{a}_m \cdot \hat{b}_n \}) . -\label{ugeneral} \end{equation} % -Allen and Germano,\cite{Allen:2006fk} showed that if the energy is -written in this form, the forces come out relatively cleanly, -% -\begin{equation} -\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} -= \frac{\partial U}{\partial r} \hat{r} - + \sum_m \left[ -\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} -\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} -+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} -\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} -\right] \label{forceequation}. -\end{equation} -% -Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ -is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} -In simplifying the algebra, we have also used: -% -\begin{align} -\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} -=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} -\right) \\ -\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} -=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} -\right) . -\end{align} -% We list below the force equations written in terms of lab-frame coordinates. The radial functions used in the two methods are listed in Table \ref{tab:tableFORCE} @@ -1134,8 +1156,8 @@ w_i(r) \begin{split} \mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& \Bigl[ -\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r} -+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\ +\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} ++ 2 \mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] w_k(r) \hat{r} \\ % 2 &+ \Bigl[ 2\text{Tr}\mathbf{Q}_{\mathbf{b}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} ) @@ -1171,29 +1193,8 @@ When energies are written in the form of Eq.~({\ref{ug % Torques SECTION ----------------------------------------------------------------------------------------- % \subsection{Torques} -When energies are written in the form of Eq.~({\ref{ugeneral}), then - torques can be found in a relatively straightforward - manner,\cite{Allen:2006fk} + % -\begin{eqnarray} -\mathbf{\tau}_{\bf a} = - \sum_m -\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} -( \hat{r} \times \hat{a}_m ) --\sum_{mn} -\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} -(\hat{a}_m \times \hat{b}_n) \\ -% -\mathbf{\tau}_{\bf b} = - \sum_m -\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} -( \hat{r} \times \hat{b}_m) -+\sum_{mn} -\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} -(\hat{a}_m \times \hat{b}_n) . -\end{eqnarray} -% -% The torques for both the Taylor-Shifted as well as Gradient-Shifted methods are given in space-frame coordinates: % @@ -1355,18 +1356,40 @@ above. $\mathbf{b}$ can be obtained by swapping indices in the expressions above. +\section{Related real-space methods} +One can also formulate an extension of the Wolf approach for point +multipoles by simply projecting the image multipole onto the surface +of the cutoff sphere, and including the 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(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] +\label{eq:SP} +\end{equation} +where the sum describes separate potential shifting that is applied to +each orientational contribution to the energy. + +The energies and torques for the shifted potential (SP) can be easily +obtained by zeroing out the $(r-r_c)$ terms in the final column of +table \ref{tab:tableenergy}. Forces for the SP method retain +discontinuities at the cutoff sphere, and can be obtained by +eliminating all functions that depend on $r_c$ in the last column of +table \ref{tab:tableFORCE}. The self-energy contributions for the SP +potential are identical to both the GSF and TSF methods. + \section{Comparison to known multipolar energies} To understand how these new real-space multipole methods behave in computer simulations, it is vital to test against established methods for computing electrostatic interactions in periodic systems, and to evaluate the size and sources of any errors that arise from the -real-space cutoffs. In this paper we test Taylor-shifted and -Gradient-shifted electrostatics against analytical methods for -computing the energies of ordered multipolar arrays. In the following -paper, we test the new methods against the multipolar Ewald sum for -computing the energies, forces and torques for a wide range of typical -condensed-phase (disordered) systems. +real-space cutoffs. In this paper we test both TSF and GSF +electrostatics against analytical methods for computing the energies +of ordered multipolar arrays. In the following paper, we test the new +methods against the multipolar Ewald sum for computing the energies, +forces and torques for a wide range of typical condensed-phase +(disordered) systems. Because long-range electrostatic effects can be significant in crystalline materials, ordered multipolar arrays present one of the @@ -1376,76 +1399,61 @@ and other periodic structures. We have repeated the L magnetization and obtained a number of these constants.\cite{Sauer} This theory was developed more completely by Luttinger and Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays -and other periodic structures. We have repeated the Luttinger \& -Tisza series summations to much higher order and obtained the energy -constants (converged to one part in $10^9$) in table \ref{tab:LT}. +and other periodic structures. -\begin{table*}[h] -\centering{ - \caption{Luttinger \& Tisza arrays and their associated - energy constants. Type "A" arrays have nearest neighbor strings of - antiparallel dipoles. Type "B" arrays have nearest neighbor - strings of antiparallel dipoles if the dipoles are contained in a - plane perpendicular to the dipole direction that passes through - the dipole.} -} -\label{tab:LT} -\begin{ruledtabular} -\begin{tabular}{cccc} -Array Type & Lattice & Dipole Direction & Energy constants \\ \hline - A & SC & 001 & -2.676788684 \\ - A & BCC & 001 & 0 \\ - A & BCC & 111 & -1.770078733 \\ - A & FCC & 001 & 2.166932835 \\ - A & FCC & 011 & -1.083466417 \\ - B & SC & 001 & -2.676788684 \\ - B & BCC & 001 & -1.338394342 \\ - B & BCC & 111 & -1.770078733 \\ - B & FCC & 001 & -1.083466417 \\ - B & FCC & 011 & -1.807573634 \\ - -- & BCC & minimum & -1.985920929 \\ -\end{tabular} -\end{ruledtabular} -\end{table*} - -In addition to the A and B arrays, there is an additional minimum +To test the new electrostatic methods, we have constructed very large, +$N=$ 16,000~(bcc) arrays of dipoles in the orientations described in +Ref. \onlinecite{LT}. These structures include ``A'' lattices with +nearest neighbor chains of antiparallel dipoles, as well as ``B'' +lattices with nearest neighbor strings of antiparallel dipoles if the +dipoles are contained in a plane perpendicular to the dipole direction +that passes through the dipole. We have also studied the minimum energy structure for the BCC lattice that was found by Luttinger \& Tisza. The total electrostatic energy for any of the arrays is given by: \begin{equation} E = C N^2 \mu^2 \end{equation} -where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is -the number of dipoles per unit volume, and $\mu$ is the strength of -the dipole. +where $C$ is the energy constant (equivalent to the Madelung +constant), $N$ is the number of dipoles per unit volume, and $\mu$ is +the strength of the dipole. Energy constants (converged to 1 part in +$10^9$) are given in the supplemental information. -To test the new electrostatic methods, we have constructed very large, -$N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in -the orientations described in table \ref{tab:LT}. For the purposes of -testing the energy expressions and the self-neutralization schemes, -the primary quantity of interest is the analytic energy constant for -the perfect arrays. Convergence to these constants are shown as a -function of both the cutoff radius, $r_c$, and the damping parameter, -$\alpha$ in Figs. \ref{fig:energyConstVsCutoff} and XXX. We have -simultaneously tested a hard cutoff (where the kernel is simply -truncated at the cutoff radius), as well as a shifted potential (SP) -form which includes a potential-shifting and self-interaction term, -but does not shift the forces and torques smoothly at the cutoff -radius. - \begin{figure} -\includegraphics[width=4.5in]{energyConstVsCutoff} -\caption{Convergence to the analytic energy constants as a function of - cutoff radius (normalized by the lattice constant) for the different - real-space methods. The two crystals shown here are the ``B'' array - for bcc crystals with the dipoles along the 001 direction (upper), - as well as the minimum energy bcc lattice (lower). The analytic - energy constants are shown as a grey dashed line. The left panel - shows results for the undamped kernel ($1/r$), while the damped - error function kernel, $B_0(r)$ was used in the right panel. } -\label{fig:energyConstVsCutoff} +\includegraphics[width=\linewidth]{Dipoles_rCutNew.pdf} +\caption{Convergence of the lattice energy constants as a function of + cutoff radius (normalized by the lattice constant, $a$) for the new + real-space methods. Three dipolar crystal structures were sampled, + and the analytic energy constants for the three lattices are + indicated with grey dashed lines. The left panel shows results for + the undamped kernel ($1/r$), while the damped error function kernel, + $B_0(r)$ was used in the right panel.} +\label{fig:Dipoles_rCut} +\end{figure} + +\begin{figure} +\includegraphics[width=\linewidth]{Dipoles_alphaNew.pdf} +\caption{Convergence to the lattice energy constants as a function of + the reduced damping parameter ($\alpha^* = \alpha a$) for the + different real-space methods in the same three dipolar crystals in + Figure \ref{fig:Dipoles_rCut}. The left panel shows results for a + relatively small cutoff radius ($r_c = 4.5 a$) while a larger cutoff + radius ($r_c = 6 a$) was used in the right panel. } +\label{fig:Dipoles_alpha} \end{figure} +For the purposes of testing the energy expressions and the +self-neutralization schemes, the primary quantity of interest is the +analytic energy constant for the perfect arrays. Convergence to these +constants are shown as a function of both the cutoff radius, $r_c$, +and the damping parameter, $\alpha$ in Figs. \ref{fig:Dipoles_rCut} +and \ref{fig:Dipoles_alpha}. We have simultaneously tested a hard +cutoff (where the kernel is simply truncated at the cutoff radius), as +well as a shifted potential (SP) form which includes a +potential-shifting and self-interaction term, but does not shift the +forces and torques smoothly at the cutoff radius. The SP method is +essentially an extension of the original Wolf method for multipoles. + The Hard cutoff exhibits oscillations around the analytic energy constants, and converges to incorrect energies when the complementary error function damping kernel is used. The shifted potential (SP) and @@ -1455,36 +1463,51 @@ cutoff region to provide accurate measures of the ener for obtaining accurate energies. The Taylor-shifted force (TSF) approximation appears to perturb the potential too much inside the cutoff region to provide accurate measures of the energy constants. - - {\it Quadrupolar} analogues to the Madelung constants were first worked out by Nagai and Nakamura who computed the energies of selected quadrupole arrays based on extensions to the Luttinger and Tisza approach.\cite{Nagai01081960,Nagai01091963} We have compared the energy constants for the lowest energy configurations for linear -quadrupoles shown in table \ref{tab:NNQ} +quadrupoles. -\begin{table*} -\centering{ - \caption{Nagai and Nakamura Quadurpolar arrays}} -\label{tab:NNQ} -\begin{ruledtabular} -\begin{tabular}{ccc} - Lattice & Quadrupole Direction & Energy constants \\ \hline - SC & 111 & -8.3 \\ - BCC & 011 & -21.7 \\ - FCC & 111 & -80.5 -\end{tabular} -\end{ruledtabular} -\end{table*} - In analogy to the dipolar arrays, the total electrostatic energy for the quadrupolar arrays is: \begin{equation} - E = C \frac{3}{4} N^2 Q^2 + E = C N \frac{3\bar{Q}^2}{4a^5} \end{equation} -where $Q$ is the quadrupole moment. +where $a$ is the lattice parameter, and $\bar{Q}$ is the effective +quadrupole moment, +\begin{equation} +\bar{Q}^2 = 2 \left(3 Q : Q - (\text{Tr} Q)^2 \right) +\end{equation} +for the primitive quadrupole as defined in Eq. \ref{eq:quadrupole}. +(For the traceless quadrupole tensor, $\Theta = 3 Q - \text{Tr} Q$, +the effective moment, $\bar{Q}^2 = \frac{2}{3} \Theta : \Theta$.) +\begin{figure} +\includegraphics[width=\linewidth]{Quadrupoles_rcutConvergence.pdf} +\caption{Convergence of the lattice energy constants as a function of + cutoff radius (normalized by the lattice constant, $a$) for the new + real-space methods. Three quadrupolar crystal structures were + sampled, and the analytic energy constants for the three lattices + are indicated with grey dashed lines. The left panel shows results + for the undamped kernel ($1/r$), while the damped error function + kernel, $B_0(r)$ was used in the right panel.} +\label{fig:QuadrupolesrcutCovergence} +\end{figure} + + +\begin{figure}[!htbp] +\includegraphics[width=3.5in]{Quadrupoles_alphaConvergence-crop.pdf} +\caption{Convergence to the analytic energy constants as a function of + cutoff damping alpha for the different + real-space methods for (a) dipolar and (b) quadrupolar crystals.The energy constants for hard, SP, GSF, TSF, and analytic methods are represented by the black sold-circle, red solid-square, green solid-diamond, and grey dashed line respectively. + The left panel shows results for the undamped kernel ($1/r$), while the damped + error function kernel, $B_0(r)$ was used in the right panel. } +\label{fig:Quadrupoles_alphaCovergence-crop.pdf} +\end{figure} + + \section{Conclusion} We have presented two efficient real-space methods for computing the interactions between point multipoles. These methods have the benefit @@ -1513,7 +1536,7 @@ for a wide range of chemical environments follows in t \begin{acknowledgments} JDG acknowledges helpful discussions with Christopher Fennell. Support for this project was provided by the National - Science Foundation under grant CHE-0848243. Computational time was + Science Foundation under grant CHE-1362211. Computational time was provided by the Center for Research Computing (CRC) at the University of Notre Dame. \end{acknowledgments} @@ -1609,33 +1632,33 @@ u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) . \begin{equation} u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) . \end{equation} - +% The functions +% needed are listed schematically below: +% % +% \begin{eqnarray} +% f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ +% g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ +% h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ +% s_2 \quad s_3 \quad &s_4 \nonumber \\ +% t_3 \quad &t_4 \nonumber \\ +% &u_4 \nonumber . +% \end{eqnarray} The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and -stored on a grid for values of $r$ from $0$ to $r_c$. The functions -needed are listed schematically below: +stored on a grid for values of $r$ from $0$ to $r_c$. Using these +functions, we find % -\begin{eqnarray} -f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ -g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ -h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ -s_2 \quad s_3 \quad &s_4 \nonumber \\ -t_3 \quad &t_4 \nonumber \\ -&u_4 \nonumber . -\end{eqnarray} - -Using these functions, we find -% \begin{align} \frac{\partial f_n}{\partial r_\alpha} =&r_\alpha \frac {g_n}{r} \label{eq:b9}\\ \frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =&\delta_{\alpha \beta}\frac {g_n}{r} +r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) \\ -\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} =& +\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial r_\gamma} =& \left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + \delta_{ \beta \gamma} r_\alpha \right) -\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) -+ r_\alpha r_\beta r_\gamma +\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\ +& + r_\alpha r_\beta r_\gamma \left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\ -\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} =& +\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial + r_\gamma \partial r_\delta} =& \left( \delta_{\alpha \beta} \delta_{\gamma \delta} + \delta_{\alpha \gamma} \delta_{\beta \delta} +\delta_{ \beta \gamma} \delta_{\alpha \delta} \right) @@ -1648,7 +1671,8 @@ + \frac{t_n}{r^4} \right)\\ \left( -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5} + \frac{t_n}{r^4} \right)\\ \frac{\partial^5 f_n} -{\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} =& +{\partial r_\alpha \partial r_\beta \partial r_\gamma \partial + r_\delta \partial r_\epsilon} =& \left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon + \text{14 permutations} \right) \left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\ @@ -1670,8 +1694,8 @@ we generalize the notation of the previous appendix. rather the individual terms in the multipole interaction energies. For damped charges , this still brings into the algebra multiple derivatives of the Smith's $B_0(r)$ function. To denote these terms, -we generalize the notation of the previous appendix. For $f(r)=1/r$ -(bare Coulomb) or $f(r)=B_0(r)$ (smeared charge) +we generalize the notation of the previous appendix. For either +$f(r)=1/r$ (undamped) or $f(r)=B_0(r)$ (damped), % \begin{align} g(r)=& \frac{df}{d r}\\ @@ -1681,10 +1705,9 @@ For undamped charges, $f(r)=1/r$, Table I lists these u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} . \end{align} % -For undamped charges, $f(r)=1/r$, Table I lists these derivatives -under the column ``Bare Coulomb.'' Equations \ref{eq:b9} to -\ref{eq:b13} are still correct for GSF electrostatics if the subscript -$n$ is eliminated. +For undamped charges Table I lists these derivatives under the column +``Bare Coulomb.'' Equations \ref{eq:b9} to \ref{eq:b13} are still +correct for GSF electrostatics if the subscript $n$ is eliminated. \newpage