44 |
|
%\preprint{AIP/123-QED} |
45 |
|
|
46 |
|
\title{Real space alternatives to the Ewald |
47 |
< |
Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles} |
47 |
> |
Sum. I. Shifted electrostatics for multipoles} |
48 |
|
|
49 |
|
\author{Madan Lamichhane} |
50 |
|
\affiliation{Department of Physics, University |
152 |
|
An efficient real-space electrostatic method involves the use of a |
153 |
|
pair-wise functional form, |
154 |
|
\begin{equation} |
155 |
< |
V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) + |
156 |
< |
\sum_i V_i^\mathrm{correction} |
155 |
> |
U = \sum_i \sum_{j>i} U_\mathrm{pair}(\mathbf{r}_{ij}, \Omega_i, \Omega_j) + |
156 |
> |
\sum_i U_i^\mathrm{self} |
157 |
|
\end{equation} |
158 |
|
that is short-ranged and easily truncated at a cutoff radius, |
159 |
|
\begin{equation} |
160 |
< |
V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) = \left\{ |
160 |
> |
U_\mathrm{pair}(\mathbf{r}_{ij},\Omega_i, \Omega_j) = \left\{ |
161 |
|
\begin{array}{ll} |
162 |
< |
V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\ |
163 |
< |
0 & \quad r > r_c , |
162 |
> |
U_\mathrm{approx} (\mathbf{r}_{ij}, \Omega_i, \Omega_j) & \quad \left| \mathbf{r}_{ij} \right| \le r_c \\ |
163 |
> |
0 & \quad \left| \mathbf{r}_{ij} \right| > r_c , |
164 |
|
\end{array} |
165 |
|
\right. |
166 |
|
\end{equation} |
167 |
< |
along with an easily computed correction term ($\sum_i |
168 |
< |
V_i^\mathrm{correction}$) which has linear-scaling with the number of |
167 |
> |
along with an easily computed self-interaction term ($\sum_i |
168 |
> |
U_i^\mathrm{self}$) which scales linearly with the number of |
169 |
|
particles. Here $\Omega_i$ and $\Omega_j$ represent orientational |
170 |
< |
coordinates of the two sites. The computational efficiency, energy |
170 |
> |
coordinates of the two sites, and $\mathbf{r}_{ij}$ is the vector |
171 |
> |
between the two sites. The computational efficiency, energy |
172 |
|
conservation, and even some physical properties of a simulation can |
173 |
< |
depend dramatically on how the $V_\mathrm{approx}$ function behaves at |
173 |
> |
depend dramatically on how the $U_\mathrm{approx}$ function behaves at |
174 |
|
the cutoff radius. The goal of any approximation method should be to |
175 |
|
mimic the real behavior of the electrostatic interactions as closely |
176 |
|
as possible without sacrificing the near-linear scaling of a cutoff |
181 |
|
contained within the cutoff sphere surrounding each particle. This is |
182 |
|
accomplished by shifting the potential functions to generate image |
183 |
|
charges on the surface of the cutoff sphere for each pair interaction |
184 |
< |
computed within $r_c$. Damping using a complementary error |
185 |
< |
function is applied to the potential to accelerate convergence. The |
186 |
< |
potential for the DSF method, where $\alpha$ is the adjustable damping |
186 |
< |
parameter, is given by |
184 |
> |
computed within $r_c$. Damping using a complementary error function is |
185 |
> |
applied to the potential to accelerate convergence. The interaction |
186 |
> |
for a pair of charges ($C_i$ and $C_j$) in the DSF method, |
187 |
|
\begin{equation*} |
188 |
< |
V_\mathrm{DSF}(r) = C_i C_j \Biggr{[} |
188 |
> |
U_\mathrm{DSF}(r_{ij}) = C_i C_j \Biggr{[} |
189 |
|
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} |
190 |
|
- \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2} |
191 |
|
+ \frac{2\alpha}{\pi^{1/2}} |
193 |
|
\right)\left(r_{ij}-r_c\right)\ \Biggr{]} |
194 |
|
\label{eq:DSFPot} |
195 |
|
\end{equation*} |
196 |
< |
Note that in this potential and in all electrostatic quantities that |
197 |
< |
follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for |
198 |
< |
clarity. |
196 |
> |
where $\alpha$ is the adjustable damping parameter. Note that in this |
197 |
> |
potential and in all electrostatic quantities that follow, the |
198 |
> |
standard $1/4 \pi \epsilon_{0}$ has been omitted for clarity. |
199 |
|
|
200 |
|
To insure net charge neutrality within each cutoff sphere, an |
201 |
|
additional ``self'' term is added to the potential. This term is |
244 |
|
a$. Then the electrostatic potential of object $\bf a$ at |
245 |
|
$\mathbf{r}$ is given by |
246 |
|
\begin{equation} |
247 |
< |
V_a(\mathbf r) = |
247 |
> |
\phi_a(\mathbf r) = |
248 |
|
\sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}. |
249 |
|
\end{equation} |
250 |
|
The Taylor expansion in $r$ can be written using an implied summation |
251 |
|
notation. Here Greek indices are used to indicate space coordinates |
252 |
|
($x$, $y$, $z$) and the subscripts $k$ and $j$ are reserved for |
253 |
< |
labelling specific charges in $\bf a$ and $\bf b$ respectively. The |
253 |
> |
labeling specific charges in $\bf a$ and $\bf b$ respectively. The |
254 |
|
Taylor expansion, |
255 |
|
\begin{equation} |
256 |
|
\frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = |
263 |
|
can then be used to express the electrostatic potential on $\bf a$ in |
264 |
|
terms of multipole operators, |
265 |
|
\begin{equation} |
266 |
< |
V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} |
266 |
> |
\phi_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} |
267 |
|
\end{equation} |
268 |
|
where |
269 |
|
\begin{equation} |
276 |
|
a}\alpha\beta}$, respectively. These are the primitive multipoles |
277 |
|
which can be expressed as a distribution of charges, |
278 |
|
\begin{align} |
279 |
< |
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\ |
280 |
< |
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha} ,\\ |
281 |
< |
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k r_{k\alpha} r_{k\beta} . |
279 |
> |
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \label{eq:charge} \\ |
280 |
> |
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha}, \label{eq:dipole}\\ |
281 |
> |
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k |
282 |
> |
r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} |
283 |
|
\end{align} |
284 |
|
Note that the definition of the primitive quadrupole here differs from |
285 |
|
the standard traceless form, and contains an additional Taylor-series |
286 |
< |
based factor of $1/2$. |
286 |
> |
based factor of $1/2$. We are essentially treating the mass |
287 |
> |
distribution with higher priority; the moment of inertia tensor, |
288 |
> |
$\overleftrightarrow{\mathsf I}$, is diagonalized to obtain body-fixed |
289 |
> |
axes, and the charge distribution may result in a quadrupole tensor |
290 |
> |
that is not necessarily diagonal in the body frame. Additional |
291 |
> |
reasons for utilizing the primitive quadrupole are discussed in |
292 |
> |
section \ref{sec:damped}. |
293 |
|
|
294 |
|
It is convenient to locate charges $q_j$ relative to the center of mass of $\bf b$. Then with $\bf{r}$ pointing from |
295 |
|
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by |
313 |
|
of $\bf a$ interacting with the same multipoles on $\bf b$. |
314 |
|
|
315 |
|
\subsection{Damped Coulomb interactions} |
316 |
+ |
\label{sec:damped} |
317 |
|
In the standard multipole expansion, one typically uses the bare |
318 |
|
Coulomb potential, with radial dependence $1/r$, as shown in |
319 |
|
Eq.~(\ref{kernel}). It is also quite common to use a damped Coulomb |
329 |
|
either $1/r$ or $B_0(r)$, and all of the techniques can be applied to |
330 |
|
bare or damped Coulomb kernels (or any other function) as long as |
331 |
|
derivatives of these functions are known. Smith's convenient |
332 |
< |
functions $B_l(r)$ are summarized in Appendix A. |
332 |
> |
functions $B_l(r)$ are summarized in Appendix A. (N.B. there is one |
333 |
> |
important distinction between the two kernels, which is the behavior |
334 |
> |
of $\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$. The former |
335 |
> |
is zero everywhere except for a delta function evaluated at the |
336 |
> |
origin. The latter also has delta function behavior, but is non-zero |
337 |
> |
for $r \neq 0$. Thus the standard justification for using a traceless |
338 |
> |
quadrupole tensor fails for the damped case.) |
339 |
|
|
340 |
|
The main goal of this work is to smoothly cut off the interaction |
341 |
|
energy as well as forces and torques as $r\rightarrow r_c$. To |
432 |
|
In general, we can write |
433 |
|
% |
434 |
|
\begin{equation} |
435 |
< |
U= (\text{prefactor}) (\text{derivatives}) f_n(r) |
435 |
> |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
436 |
|
\label{generic} |
437 |
|
\end{equation} |
438 |
|
% |
453 |
|
of the same index $n$. The algebra required to evaluate energies, |
454 |
|
forces and torques is somewhat tedious, so only the final forms are |
455 |
|
presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}. |
456 |
+ |
One of the principal findings of our work is that the individual |
457 |
+ |
orientational contributions to the various multipole-multipole |
458 |
+ |
interactions must be treated with distinct radial functions, but each |
459 |
+ |
of these contributions is independently force shifted at the cutoff |
460 |
+ |
radius. |
461 |
|
|
462 |
|
\subsection{Gradient-shifted force (GSF) electrostatics} |
463 |
|
The second, and conceptually simpler approach to force-shifting |
465 |
|
expansion, and has a similar interaction energy for all multipole |
466 |
|
orders: |
467 |
|
\begin{equation} |
468 |
< |
U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big |
469 |
< |
\lvert _{r_c} . |
468 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
469 |
> |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
470 |
> |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
471 |
|
\label{generic2} |
472 |
|
\end{equation} |
473 |
< |
Here the gradient for force shifting is evaluated for an image |
474 |
< |
multipole projected onto the surface of the cutoff sphere (see fig |
475 |
< |
\ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$ |
476 |
< |
appear. The primary difference between the TSF and GSF methods is the |
477 |
< |
stage at which the Taylor Series is applied; in the Taylor-shifted |
478 |
< |
approach, it is applied to the kernel itself. In the Gradient-shifted |
479 |
< |
approach, it is applied to individual radial interactions terms in the |
480 |
< |
multipole expansion. Energies from this method thus have the general |
481 |
< |
form: |
473 |
> |
where the sum describes a separate force-shifting that is applied to |
474 |
> |
each orientational contribution to the energy. Both the potential and |
475 |
> |
the gradient for force shifting are evaluated for an image multipole |
476 |
> |
projected onto the surface of the cutoff sphere (see fig |
477 |
> |
\ref{fig:shiftedMultipoles}). The image multipole retains the |
478 |
> |
orientation ($\hat{\mathbf{b}}$) of the interacting multipole. No |
479 |
> |
higher order terms $(r-r_c)^n$ appear. The primary difference between |
480 |
> |
the TSF and GSF methods is the stage at which the Taylor Series is |
481 |
> |
applied; in the Taylor-shifted approach, it is applied to the kernel |
482 |
> |
itself. In the Gradient-shifted approach, it is applied to individual |
483 |
> |
radial interactions terms in the multipole expansion. Energies from |
484 |
> |
this method thus have the general form: |
485 |
|
\begin{equation} |
486 |
|
U= \sum (\text{angular factor}) (\text{radial factor}). |
487 |
|
\label{generic3} |
494 |
|
the radial factors differ between the two methods. |
495 |
|
|
496 |
|
\subsection{\label{sec:level2}Body and space axes} |
497 |
+ |
Although objects $\bf a$ and $\bf b$ rotate during a molecular |
498 |
+ |
dynamics (MD) simulation, their multipole tensors remain fixed in |
499 |
+ |
body-frame coordinates. While deriving force and torque expressions, |
500 |
+ |
it is therefore convenient to write the energies, forces, and torques |
501 |
+ |
in intermediate forms involving the vectors of the rotation matrices. |
502 |
+ |
We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors |
503 |
+ |
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$. |
504 |
+ |
In a typical simulation , the initial axes are obtained by |
505 |
+ |
diagonalizing the moment of inertia tensors for the objects. (N.B., |
506 |
+ |
the body axes are generally {\it not} the same as those for which the |
507 |
+ |
quadrupole moment is diagonal.) The rotation matrices are then |
508 |
+ |
propagated during the simulation. |
509 |
|
|
510 |
< |
[XXX Do we need this section in the main paper? or should it go in the |
476 |
< |
extra materials?] |
477 |
< |
|
478 |
< |
So far, all energies and forces have been written in terms of fixed |
479 |
< |
space coordinates. Interaction energies are computed from the generic |
480 |
< |
formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine |
481 |
< |
orientational prefactors with radial functions. Because objects $\bf |
482 |
< |
a$ and $\bf b$ both translate and rotate during a molecular dynamics |
483 |
< |
(MD) simulation, it is desirable to contract all $r$-dependent terms |
484 |
< |
with dipole and quadrupole moments expressed in terms of their body |
485 |
< |
axes. To do so, we have followed the methodology of Allen and |
486 |
< |
Germano,\cite{Allen:2006fk} which was itself based on earlier work by |
487 |
< |
Price {\em et al.}\cite{Price:1984fk} |
488 |
< |
|
489 |
< |
We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors |
490 |
< |
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$ |
491 |
< |
referring to a convenient set of inertial body axes. (N.B., these |
492 |
< |
body axes are generally not the same as those for which the quadrupole |
493 |
< |
moment is diagonal.) Then, |
494 |
< |
% |
495 |
< |
\begin{eqnarray} |
496 |
< |
\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z} \\ |
497 |
< |
\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z} . |
498 |
< |
\end{eqnarray} |
499 |
< |
Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be |
510 |
> |
The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be |
511 |
|
expressed using these unit vectors: |
512 |
|
\begin{eqnarray} |
513 |
|
\hat{\mathbf {a}} = |
515 |
|
\hat{a}_1 \\ |
516 |
|
\hat{a}_2 \\ |
517 |
|
\hat{a}_3 |
518 |
< |
\end{pmatrix} |
508 |
< |
= |
509 |
< |
\begin{pmatrix} |
510 |
< |
a_{1x} \quad a_{1y} \quad a_{1z} \\ |
511 |
< |
a_{2x} \quad a_{2y} \quad a_{2z} \\ |
512 |
< |
a_{3x} \quad a_{3y} \quad a_{3z} |
513 |
< |
\end{pmatrix}\\ |
518 |
> |
\end{pmatrix}, \qquad |
519 |
|
\hat{\mathbf {b}} = |
520 |
|
\begin{pmatrix} |
521 |
|
\hat{b}_1 \\ |
522 |
|
\hat{b}_2 \\ |
523 |
|
\hat{b}_3 |
524 |
|
\end{pmatrix} |
520 |
– |
= |
521 |
– |
\begin{pmatrix} |
522 |
– |
b_{1x} \quad b_{1y} \quad b_{1z} \\ |
523 |
– |
b_{2x} \quad b_{2y} \quad b_{2z} \\ |
524 |
– |
b_{3x} \quad b_{3y} \quad b_{3z} |
525 |
– |
\end{pmatrix} . |
525 |
|
\end{eqnarray} |
526 |
|
% |
527 |
|
These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$ |
528 |
< |
coordinates. All contractions of prefactors with derivatives of |
529 |
< |
functions can be written in terms of these matrices. It proves to be |
530 |
< |
equally convenient to just write any contraction in terms of unit |
531 |
< |
vectors $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$. In the torque |
532 |
< |
expressions, it is useful to have the angular-dependent terms |
533 |
< |
available in three different fashions, e.g. for the dipole-dipole |
535 |
< |
contraction: |
528 |
> |
coordinates. |
529 |
> |
|
530 |
> |
Allen and Germano,\cite{Allen:2006fk} following earlier work by Price |
531 |
> |
{\em et al.},\cite{Price:1984fk} showed that if the interaction |
532 |
> |
energies are written explicitly in terms of $\hat{r}$ and the body |
533 |
> |
axes ($\hat{a}_m$, $\hat{b}_n$) : |
534 |
|
% |
535 |
|
\begin{equation} |
536 |
+ |
U(r, \{\hat{a}_m \cdot \hat{r} \}, |
537 |
+ |
\{\hat{b}_n\cdot \hat{r} \}, |
538 |
+ |
\{\hat{a}_m \cdot \hat{b}_n \}) . |
539 |
+ |
\label{ugeneral} |
540 |
+ |
\end{equation} |
541 |
+ |
% |
542 |
+ |
the forces come out relatively cleanly, |
543 |
+ |
% |
544 |
+ |
\begin{equation} |
545 |
+ |
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} |
546 |
+ |
= \frac{\partial U}{\partial r} \hat{r} |
547 |
+ |
+ \sum_m \left[ |
548 |
+ |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
549 |
+ |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
550 |
+ |
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
551 |
+ |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
552 |
+ |
\right] \label{forceequation}. |
553 |
+ |
\end{equation} |
554 |
+ |
|
555 |
+ |
The torques can also be found in a relatively similar |
556 |
+ |
manner, |
557 |
+ |
% |
558 |
+ |
\begin{eqnarray} |
559 |
+ |
\mathbf{\tau}_{\bf a} = |
560 |
+ |
\sum_m |
561 |
+ |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
562 |
+ |
( \hat{r} \times \hat{a}_m ) |
563 |
+ |
-\sum_{mn} |
564 |
+ |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
565 |
+ |
(\hat{a}_m \times \hat{b}_n) \\ |
566 |
+ |
% |
567 |
+ |
\mathbf{\tau}_{\bf b} = |
568 |
+ |
\sum_m |
569 |
+ |
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
570 |
+ |
( \hat{r} \times \hat{b}_m) |
571 |
+ |
+\sum_{mn} |
572 |
+ |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
573 |
+ |
(\hat{a}_m \times \hat{b}_n) . |
574 |
+ |
\end{eqnarray} |
575 |
+ |
|
576 |
+ |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ |
577 |
+ |
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} |
578 |
+ |
We also made use of the identities, |
579 |
+ |
% |
580 |
+ |
\begin{align} |
581 |
+ |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
582 |
+ |
=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} |
583 |
+ |
\right) \\ |
584 |
+ |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
585 |
+ |
=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} |
586 |
+ |
\right) . |
587 |
+ |
\end{align} |
588 |
+ |
|
589 |
+ |
Many of the multipole contractions required can be written in one of |
590 |
+ |
three equivalent forms using the unit vectors $\hat{r}$, $\hat{a}_m$, |
591 |
+ |
and $\hat{b}_n$. In the torque expressions, it is useful to have the |
592 |
+ |
angular-dependent terms available in all three fashions, e.g. for the |
593 |
+ |
dipole-dipole contraction: |
594 |
+ |
% |
595 |
+ |
\begin{equation} |
596 |
|
\mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} |
597 |
|
= D_{\bf {a}\alpha} D_{\bf {b}\alpha} = |
598 |
|
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} |
604 |
|
explicit sums over body indices and the dot products now indicate |
605 |
|
contractions using space indices. |
606 |
|
|
607 |
+ |
In computing our force and torque expressions, we carried out most of |
608 |
+ |
the work in body coordinates, and have transformed the expressions |
609 |
+ |
back to space-frame coordinates, which are reported below. Interested |
610 |
+ |
readers may consult the supplemental information for this paper for |
611 |
+ |
the intermediate body-frame expressions. |
612 |
|
|
613 |
|
\subsection{The Self-Interaction \label{sec:selfTerm}} |
614 |
|
|
615 |
|
In addition to cutoff-sphere neutralization, the Wolf |
616 |
|
summation~\cite{Wolf99} and the damped shifted force (DSF) |
617 |
< |
extension~\cite{Fennell:2006zl} also included self-interactions that |
617 |
> |
extension~\cite{Fennell:2006zl} also include self-interactions that |
618 |
|
are handled separately from the pairwise interactions between |
619 |
|
sites. The self-term is normally calculated via a single loop over all |
620 |
|
sites in the system, and is relatively cheap to evaluate. The |
626 |
|
the cutoff sphere. For a system of undamped charges, the total |
627 |
|
self-term is |
628 |
|
\begin{equation} |
629 |
< |
V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2} |
629 |
> |
U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2} |
630 |
|
\label{eq:selfTerm} |
631 |
|
\end{equation} |
632 |
|
|
641 |
|
complexity to the Ewald sum). For a system containing only damped |
642 |
|
charges, the complete self-interaction can be written as |
643 |
|
\begin{equation} |
644 |
< |
V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
644 |
> |
U_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
645 |
|
\frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N |
646 |
|
C_{\bf a}^{2}. |
647 |
|
\label{eq:dampSelfTerm} |
661 |
|
reciprocal-space portion is identical to half of the self-term |
662 |
|
obtained by Smith and Aguado and Madden for the application of the |
663 |
|
Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given |
664 |
< |
site which posesses a charge, dipole, and multipole, both types of |
664 |
> |
site which posesses a charge, dipole, and quadrupole, both types of |
665 |
|
contribution are given in table \ref{tab:tableSelf}. |
666 |
|
|
667 |
|
\begin{table*} |
674 |
|
Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\ |
675 |
|
Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) + |
676 |
|
\frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\ |
677 |
< |
Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
677 |
> |
Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
678 |
|
$- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ & |
679 |
|
$-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\ |
680 |
|
Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left( |
804 |
|
U_{Q_{\bf a}Q_{\bf b}}(r)=& |
805 |
|
\Bigl[ |
806 |
|
\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}} |
807 |
< |
+2 \text{Tr} \left( |
808 |
< |
\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r) |
807 |
> |
+2 |
808 |
> |
\mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r) |
809 |
|
\\ |
810 |
|
% 2 |
811 |
|
&+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}} |
838 |
|
\caption{\label{tab:tableenergy}Radial functions used in the energy |
839 |
|
and torque equations. The $f, g, h, s, t, \mathrm{and} u$ |
840 |
|
functions used in this table are defined in Appendices B and C.} |
841 |
< |
\begin{tabular}{|c|c|l|l|} \hline |
842 |
< |
Generic&Bare Coulomb&Taylor-Shifted&Gradient-Shifted |
841 |
> |
\begin{tabular}{|c|c|l|l|l|} \hline |
842 |
> |
Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF) |
843 |
|
\\ \hline |
844 |
|
% |
845 |
|
% |
848 |
|
$v_{01}(r)$ & |
849 |
|
$\frac{1}{r}$ & |
850 |
|
$f_0(r)$ & |
851 |
< |
$f(r)-f(r_c)-(r-r_c)g(r_c)$ |
851 |
> |
$f(r)-f(r_c)$ & |
852 |
> |
SP $-(r-r_c)g(r_c)$ |
853 |
|
\\ |
854 |
|
% |
855 |
|
% |
858 |
|
$v_{11}(r)$ & |
859 |
|
$-\frac{1}{r^2}$ & |
860 |
|
$g_1(r)$ & |
861 |
< |
$g(r)-g(r_c)-(r-r_c)h(r_c)$ \\ |
861 |
> |
$g(r)-g(r_c)$ & |
862 |
> |
SP $-(r-r_c)h(r_c)$ \\ |
863 |
|
% |
864 |
|
% |
865 |
|
% |
867 |
|
$v_{21}(r)$ & |
868 |
|
$-\frac{1}{r^3} $ & |
869 |
|
$\frac{g_2(r)}{r} $ & |
870 |
< |
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c} -(r-r_c) |
870 |
> |
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c}$ & |
871 |
> |
SP $-(r-r_c) |
872 |
|
\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
873 |
|
$v_{22}(r)$ & |
874 |
|
$\frac{3}{r^3} $ & |
875 |
|
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ & |
876 |
< |
$\left(-\frac{g(r)}{r}+h(r) \right) |
877 |
< |
-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ \\ |
812 |
< |
&&& $ ~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$ |
813 |
< |
\\ |
876 |
> |
$\left(-\frac{g(r)}{r}+h(r) \right)$ & SP \\ |
877 |
> |
&&& $~~~-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ & $~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\ |
878 |
|
% |
879 |
|
% |
880 |
|
% |
882 |
|
$v_{31}(r)$ & |
883 |
|
$\frac{3}{r^4} $ & |
884 |
|
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ & |
885 |
< |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right) |
886 |
< |
-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\ |
823 |
< |
&&&$ ~~~ -(r-r_c) \left(\frac{2g(r_c)}{r_c^3}-\frac{2h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ |
824 |
< |
\\ |
885 |
> |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)$ & SP \\ |
886 |
> |
&&& $-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $ & $~~~-(r-r_c) \left(\frac{2g(r_c)}{r_c^3}-\frac{2h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\ |
887 |
|
% |
888 |
|
$v_{32}(r)$ & |
889 |
|
$-\frac{15}{r^4} $ & |
890 |
|
$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ & |
891 |
< |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right) |
892 |
< |
- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\ |
893 |
< |
&&&$ ~~~ -(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}-\frac{3s(r_c)}{r_c}+t(r_c) \right)$ |
832 |
< |
\\ |
891 |
> |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$ & SP \\ |
892 |
> |
&&& $~~~- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) |
893 |
> |
\right)$ & $~~~-(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}-\frac{3s(r_c)}{r_c}+t(r_c) \right)$ \\ |
894 |
|
% |
895 |
|
% |
896 |
|
% |
898 |
|
$v_{41}(r)$ & |
899 |
|
$\frac{3}{r^5} $ & |
900 |
|
$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $ & |
901 |
< |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right) |
902 |
< |
- \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\ |
842 |
< |
&&&$ ~~~ -(r-r_c) \left( \frac{3g(r_c)}{r_c^4}-\frac{3h(r_c)}{r_c^3}+\frac{s(r_c)}{r_c^2} \right)$ |
901 |
> |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)$ & SP \\ |
902 |
> |
&&& $~~~- \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$& $~~~-(r-r_c) \left( \frac{3g(r_c)}{r_c^4}-\frac{3h(r_c)}{r_c^3}+\frac{s(r_c)}{r_c^2} \right)$ |
903 |
|
\\ |
904 |
|
% 2 |
905 |
|
$v_{42}(r)$ & |
906 |
|
$- \frac{15}{r^5} $ & |
907 |
|
$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ & |
908 |
< |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right) |
909 |
< |
-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\ |
910 |
< |
&&&$ ~~~ -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3} |
911 |
< |
-\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$ |
912 |
< |
\\ |
908 |
> |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} |
909 |
> |
\right)$ & SP \\ |
910 |
> |
&&& $~~~-\left( \frac{3g(r_c)}{r_c^3} - |
911 |
> |
\frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ & $~~~-(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3} |
912 |
> |
-\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$ \\ |
913 |
|
% 3 |
914 |
|
$v_{43}(r)$ & |
915 |
|
$ \frac{105}{r^5} $ & |
916 |
|
$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ & |
917 |
< |
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\ |
918 |
< |
&&&$~~~ -\left(-\frac{15g(r_c)}{r_c^3}+\frac{15h(r_c)}{r_c^2}-\frac{6s(r_c)}{r_c} + t(r_c)\right)$ \\ |
919 |
< |
&&&$~~~ -(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}+\frac{21s(r_c)}{r_c^2} |
920 |
< |
-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline |
917 |
> |
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + |
918 |
> |
t(r)\right)$ & SP $-(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}\right.$\\ |
919 |
> |
&&& |
920 |
> |
$~~~-\left(-\frac{15g(r_c)}{r_c^3}+\frac{15h(r_c)}{r_c^2}-\frac{6s(r_c)}{r_c} |
921 |
> |
+ t(r_c)\right)$ & $~~~~~~~\left.+\frac{21s(r_c)}{r_c^2} |
922 |
> |
-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ |
923 |
> |
\hline |
924 |
|
\end{tabular} |
925 |
|
\end{sidewaystable} |
926 |
|
% |
1044 |
|
F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r} |
1045 |
|
\quad \text{and} \quad F_{\bf b \alpha} |
1046 |
|
= - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r} . |
984 |
– |
\end{equation} |
985 |
– |
% |
986 |
– |
Obtaining the force from the interaction energy expressions is the |
987 |
– |
same for higher-order multipole interactions -- the trick is to make |
988 |
– |
sure that all $r$-dependent derivatives are considered. This is |
989 |
– |
straighforward if the interaction energies are written explicitly in |
990 |
– |
terms of $\hat{r}$ and the body axes ($\hat{a}_m$, |
991 |
– |
$\hat{b}_n$) : |
992 |
– |
% |
993 |
– |
\begin{equation} |
994 |
– |
U(r,\{\hat{a}_m \cdot \hat{r} \}, |
995 |
– |
\{\hat{b}_n\cdot \hat{r} \}, |
996 |
– |
\{\hat{a}_m \cdot \hat{b}_n \}) . |
997 |
– |
\label{ugeneral} |
998 |
– |
\end{equation} |
999 |
– |
% |
1000 |
– |
Allen and Germano,\cite{Allen:2006fk} showed that if the energy is |
1001 |
– |
written in this form, the forces come out relatively cleanly, |
1002 |
– |
% |
1003 |
– |
\begin{equation} |
1004 |
– |
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} |
1005 |
– |
= \frac{\partial U}{\partial r} \hat{r} |
1006 |
– |
+ \sum_m \left[ |
1007 |
– |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
1008 |
– |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
1009 |
– |
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
1010 |
– |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
1011 |
– |
\right] \label{forceequation}. |
1047 |
|
\end{equation} |
1048 |
|
% |
1014 |
– |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ |
1015 |
– |
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} |
1016 |
– |
In simplifying the algebra, we have also used: |
1017 |
– |
% |
1018 |
– |
\begin{align} |
1019 |
– |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
1020 |
– |
=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} |
1021 |
– |
\right) \\ |
1022 |
– |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
1023 |
– |
=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} |
1024 |
– |
\right) . |
1025 |
– |
\end{align} |
1026 |
– |
% |
1049 |
|
We list below the force equations written in terms of lab-frame |
1050 |
|
coordinates. The radial functions used in the two methods are listed |
1051 |
|
in Table \ref{tab:tableFORCE} |
1156 |
|
\begin{split} |
1157 |
|
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& |
1158 |
|
\Bigl[ |
1159 |
< |
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r} |
1160 |
< |
+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\ |
1159 |
> |
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} |
1160 |
> |
+ 2 \mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] w_k(r) \hat{r} \\ |
1161 |
|
% 2 |
1162 |
|
&+ \Bigl[ |
1163 |
|
2\text{Tr}\mathbf{Q}_{\mathbf{b}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} ) |
1193 |
|
% Torques SECTION ----------------------------------------------------------------------------------------- |
1194 |
|
% |
1195 |
|
\subsection{Torques} |
1196 |
< |
When energies are written in the form of Eq.~({\ref{ugeneral}), then |
1175 |
< |
torques can be found in a relatively straightforward |
1176 |
< |
manner,\cite{Allen:2006fk} |
1177 |
< |
% |
1178 |
< |
\begin{eqnarray} |
1179 |
< |
\mathbf{\tau}_{\bf a} = |
1180 |
< |
\sum_m |
1181 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
1182 |
< |
( \hat{r} \times \hat{a}_m ) |
1183 |
< |
-\sum_{mn} |
1184 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
1185 |
< |
(\hat{a}_m \times \hat{b}_n) \\ |
1186 |
< |
% |
1187 |
< |
\mathbf{\tau}_{\bf b} = |
1188 |
< |
\sum_m |
1189 |
< |
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
1190 |
< |
( \hat{r} \times \hat{b}_m) |
1191 |
< |
+\sum_{mn} |
1192 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
1193 |
< |
(\hat{a}_m \times \hat{b}_n) . |
1194 |
< |
\end{eqnarray} |
1196 |
> |
|
1197 |
|
% |
1196 |
– |
% |
1198 |
|
The torques for both the Taylor-Shifted as well as Gradient-Shifted |
1199 |
|
methods are given in space-frame coordinates: |
1200 |
|
% |
1356 |
|
$\mathbf{b}$ can be obtained by swapping indices in the expressions |
1357 |
|
above. |
1358 |
|
|
1359 |
+ |
\section{Related real-space methods} |
1360 |
+ |
One can also formulate an extension of the Wolf approach for point |
1361 |
+ |
multipoles by simply projecting the image multipole onto the surface |
1362 |
+ |
of the cutoff sphere, and including the interactions with the central |
1363 |
+ |
multipole and the image. This effectively shifts the total potential |
1364 |
+ |
to zero at the cutoff radius, |
1365 |
+ |
\begin{equation} |
1366 |
+ |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
1367 |
+ |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
1368 |
+ |
\label{eq:SP} |
1369 |
+ |
\end{equation} |
1370 |
+ |
where the sum describes separate potential shifting that is applied to |
1371 |
+ |
each orientational contribution to the energy. |
1372 |
+ |
|
1373 |
+ |
The energies and torques for the shifted potential (SP) can be easily |
1374 |
+ |
obtained by zeroing out the $(r-r_c)$ terms in the final column of |
1375 |
+ |
table \ref{tab:tableenergy}. Forces for the SP method retain |
1376 |
+ |
discontinuities at the cutoff sphere, and can be obtained by |
1377 |
+ |
eliminating all functions that depend on $r_c$ in the last column of |
1378 |
+ |
table \ref{tab:tableFORCE}. The self-energy contributions for the SP |
1379 |
+ |
potential are identical to both the GSF and TSF methods. |
1380 |
+ |
|
1381 |
|
\section{Comparison to known multipolar energies} |
1382 |
|
|
1383 |
|
To understand how these new real-space multipole methods behave in |
1384 |
|
computer simulations, it is vital to test against established methods |
1385 |
|
for computing electrostatic interactions in periodic systems, and to |
1386 |
|
evaluate the size and sources of any errors that arise from the |
1387 |
< |
real-space cutoffs. In this paper we test Taylor-shifted and |
1388 |
< |
Gradient-shifted electrostatics against analytical methods for |
1389 |
< |
computing the energies of ordered multipolar arrays. In the following |
1390 |
< |
paper, we test the new methods against the multipolar Ewald sum for |
1391 |
< |
computing the energies, forces and torques for a wide range of typical |
1392 |
< |
condensed-phase (disordered) systems. |
1387 |
> |
real-space cutoffs. In this paper we test both TSF and GSF |
1388 |
> |
electrostatics against analytical methods for computing the energies |
1389 |
> |
of ordered multipolar arrays. In the following paper, we test the new |
1390 |
> |
methods against the multipolar Ewald sum for computing the energies, |
1391 |
> |
forces and torques for a wide range of typical condensed-phase |
1392 |
> |
(disordered) systems. |
1393 |
|
|
1394 |
|
Because long-range electrostatic effects can be significant in |
1395 |
|
crystalline materials, ordered multipolar arrays present one of the |
1399 |
|
magnetization and obtained a number of these constants.\cite{Sauer} |
1400 |
|
This theory was developed more completely by Luttinger and |
1401 |
|
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays |
1402 |
< |
and other periodic structures. We have repeated the Luttinger \& |
1380 |
< |
Tisza series summations to much higher order and obtained the energy |
1381 |
< |
constants (converged to one part in $10^9$) in table \ref{tab:LT}. |
1402 |
> |
and other periodic structures. |
1403 |
|
|
1404 |
< |
\begin{table*}[h] |
1405 |
< |
\centering{ |
1406 |
< |
\caption{Luttinger \& Tisza arrays and their associated |
1407 |
< |
energy constants. Type "A" arrays have nearest neighbor strings of |
1408 |
< |
antiparallel dipoles. Type "B" arrays have nearest neighbor |
1409 |
< |
strings of antiparallel dipoles if the dipoles are contained in a |
1410 |
< |
plane perpendicular to the dipole direction that passes through |
1390 |
< |
the dipole.} |
1391 |
< |
} |
1392 |
< |
\label{tab:LT} |
1393 |
< |
\begin{ruledtabular} |
1394 |
< |
\begin{tabular}{cccc} |
1395 |
< |
Array Type & Lattice & Dipole Direction & Energy constants \\ \hline |
1396 |
< |
A & SC & 001 & -2.676788684 \\ |
1397 |
< |
A & BCC & 001 & 0 \\ |
1398 |
< |
A & BCC & 111 & -1.770078733 \\ |
1399 |
< |
A & FCC & 001 & 2.166932835 \\ |
1400 |
< |
A & FCC & 011 & -1.083466417 \\ |
1401 |
< |
B & SC & 001 & -2.676788684 \\ |
1402 |
< |
B & BCC & 001 & -1.338394342 \\ |
1403 |
< |
B & BCC & 111 & -1.770078733 \\ |
1404 |
< |
B & FCC & 001 & -1.083466417 \\ |
1405 |
< |
B & FCC & 011 & -1.807573634 \\ |
1406 |
< |
-- & BCC & minimum & -1.985920929 \\ |
1407 |
< |
\end{tabular} |
1408 |
< |
\end{ruledtabular} |
1409 |
< |
\end{table*} |
1410 |
< |
|
1411 |
< |
In addition to the A and B arrays, there is an additional minimum |
1404 |
> |
To test the new electrostatic methods, we have constructed very large, |
1405 |
> |
$N=$ 16,000~(bcc) arrays of dipoles in the orientations described in |
1406 |
> |
Ref. \onlinecite{LT}. These structures include ``A'' lattices with |
1407 |
> |
nearest neighbor chains of antiparallel dipoles, as well as ``B'' |
1408 |
> |
lattices with nearest neighbor strings of antiparallel dipoles if the |
1409 |
> |
dipoles are contained in a plane perpendicular to the dipole direction |
1410 |
> |
that passes through the dipole. We have also studied the minimum |
1411 |
|
energy structure for the BCC lattice that was found by Luttinger \& |
1412 |
|
Tisza. The total electrostatic energy for any of the arrays is given |
1413 |
|
by: |
1414 |
|
\begin{equation} |
1415 |
|
E = C N^2 \mu^2 |
1416 |
|
\end{equation} |
1417 |
< |
where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is |
1418 |
< |
the number of dipoles per unit volume, and $\mu$ is the strength of |
1419 |
< |
the dipole. |
1417 |
> |
where $C$ is the energy constant (equivalent to the Madelung |
1418 |
> |
constant), $N$ is the number of dipoles per unit volume, and $\mu$ is |
1419 |
> |
the strength of the dipole. Energy constants (converged to 1 part in |
1420 |
> |
$10^9$) are given in the supplemental information. |
1421 |
|
|
1422 |
< |
To test the new electrostatic methods, we have constructed very large, |
1423 |
< |
$N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in |
1424 |
< |
the orientations described in table \ref{tab:LT}. For the purposes of |
1425 |
< |
testing the energy expressions and the self-neutralization schemes, |
1426 |
< |
the primary quantity of interest is the analytic energy constant for |
1427 |
< |
the perfect arrays. Convergence to these constants are shown as a |
1428 |
< |
function of both the cutoff radius, $r_c$, and the damping parameter, |
1429 |
< |
$\alpha$ in Figs. \ref{fig:energyConstVsCutoff} and XXX. We have |
1430 |
< |
simultaneously tested a hard cutoff (where the kernel is simply |
1431 |
< |
truncated at the cutoff radius), as well as a shifted potential (SP) |
1432 |
< |
form which includes a potential-shifting and self-interaction term, |
1433 |
< |
but does not shift the forces and torques smoothly at the cutoff |
1434 |
< |
radius. |
1422 |
> |
\begin{figure} |
1423 |
> |
\includegraphics[width=\linewidth]{Dipoles_rCutNew.pdf} |
1424 |
> |
\caption{Convergence of the lattice energy constants as a function of |
1425 |
> |
cutoff radius (normalized by the lattice constant, $a$) for the new |
1426 |
> |
real-space methods. Three dipolar crystal structures were sampled, |
1427 |
> |
and the analytic energy constants for the three lattices are |
1428 |
> |
indicated with grey dashed lines. The left panel shows results for |
1429 |
> |
the undamped kernel ($1/r$), while the damped error function kernel, |
1430 |
> |
$B_0(r)$ was used in the right panel.} |
1431 |
> |
\label{fig:Dipoles_rCut} |
1432 |
> |
\end{figure} |
1433 |
|
|
1434 |
|
\begin{figure} |
1435 |
< |
\includegraphics[width=4.5in]{energyConstVsCutoff} |
1436 |
< |
\caption{Convergence to the analytic energy constants as a function of |
1437 |
< |
cutoff radius (normalized by the lattice constant) for the different |
1438 |
< |
real-space methods. The two crystals shown here are the ``B'' array |
1439 |
< |
for bcc crystals with the dipoles along the 001 direction (upper), |
1440 |
< |
as well as the minimum energy bcc lattice (lower). The analytic |
1441 |
< |
energy constants are shown as a grey dashed line. The left panel |
1442 |
< |
shows results for the undamped kernel ($1/r$), while the damped |
1445 |
< |
error function kernel, $B_0(r)$ was used in the right panel. } |
1446 |
< |
\label{fig:energyConstVsCutoff} |
1435 |
> |
\includegraphics[width=\linewidth]{Dipoles_alphaNew.pdf} |
1436 |
> |
\caption{Convergence to the lattice energy constants as a function of |
1437 |
> |
the reduced damping parameter ($\alpha^* = \alpha a$) for the |
1438 |
> |
different real-space methods in the same three dipolar crystals in |
1439 |
> |
Figure \ref{fig:Dipoles_rCut}. The left panel shows results for a |
1440 |
> |
relatively small cutoff radius ($r_c = 4.5 a$) while a larger cutoff |
1441 |
> |
radius ($r_c = 6 a$) was used in the right panel. } |
1442 |
> |
\label{fig:Dipoles_alpha} |
1443 |
|
\end{figure} |
1444 |
|
|
1445 |
< |
The Hard cutoff exhibits oscillations around the analytic energy |
1445 |
> |
For the purposes of testing the energy expressions and the |
1446 |
> |
self-neutralization schemes, the primary quantity of interest is the |
1447 |
> |
analytic energy constant for the perfect arrays. Convergence to these |
1448 |
> |
constants are shown as a function of both the cutoff radius, $r_c$, |
1449 |
> |
and the damping parameter, $\alpha$ in Figs.\ref{fig:Dipoles_rCut} |
1450 |
> |
and \ref{fig:Dipoles_alpha}. We have simultaneously tested a hard |
1451 |
> |
cutoff (where the kernel is simply truncated at the cutoff radius), as |
1452 |
> |
well as a shifted potential (SP) form which includes a |
1453 |
> |
potential-shifting and self-interaction term, but does not shift the |
1454 |
> |
forces and torques smoothly at the cutoff radius. The SP method is |
1455 |
> |
essentially an extension of the original Wolf method for multipoles. |
1456 |
> |
|
1457 |
> |
The hard cutoff exhibits oscillations around the analytic energy |
1458 |
|
constants, and converges to incorrect energies when the complementary |
1459 |
|
error function damping kernel is used. The shifted potential (SP) and |
1460 |
|
gradient-shifted force (GSF) approximations converge to the correct |
1461 |
< |
energy smoothly by $r_c / 6 a$ even for the undamped case. This |
1461 |
> |
energy smoothly by $r_c = 6 a$ even for the undamped case. This |
1462 |
|
indicates that the correction provided by the self term is required |
1463 |
|
for obtaining accurate energies. The Taylor-shifted force (TSF) |
1464 |
|
approximation appears to perturb the potential too much inside the |
1465 |
|
cutoff region to provide accurate measures of the energy constants. |
1466 |
|
|
1459 |
– |
|
1467 |
|
{\it Quadrupolar} analogues to the Madelung constants were first |
1468 |
|
worked out by Nagai and Nakamura who computed the energies of selected |
1469 |
|
quadrupole arrays based on extensions to the Luttinger and Tisza |
1470 |
< |
approach.\cite{Nagai01081960,Nagai01091963} We have compared the |
1464 |
< |
energy constants for the lowest energy configurations for linear |
1465 |
< |
quadrupoles shown in table \ref{tab:NNQ} |
1470 |
> |
approach.\cite{Nagai01081960,Nagai01091963} |
1471 |
|
|
1467 |
– |
\begin{table*} |
1468 |
– |
\centering{ |
1469 |
– |
\caption{Nagai and Nakamura Quadurpolar arrays}} |
1470 |
– |
\label{tab:NNQ} |
1471 |
– |
\begin{ruledtabular} |
1472 |
– |
\begin{tabular}{ccc} |
1473 |
– |
Lattice & Quadrupole Direction & Energy constants \\ \hline |
1474 |
– |
SC & 111 & -8.3 \\ |
1475 |
– |
BCC & 011 & -21.7 \\ |
1476 |
– |
FCC & 111 & -80.5 |
1477 |
– |
\end{tabular} |
1478 |
– |
\end{ruledtabular} |
1479 |
– |
\end{table*} |
1480 |
– |
|
1472 |
|
In analogy to the dipolar arrays, the total electrostatic energy for |
1473 |
|
the quadrupolar arrays is: |
1474 |
|
\begin{equation} |
1475 |
< |
E = C \frac{3}{4} N^2 Q^2 |
1475 |
> |
E = C N \frac{3\bar{Q}^2}{4a^5} |
1476 |
|
\end{equation} |
1477 |
< |
where $Q$ is the quadrupole moment. |
1477 |
> |
where $a$ is the lattice parameter, and $\bar{Q}$ is the effective |
1478 |
> |
quadrupole moment, |
1479 |
> |
\begin{equation} |
1480 |
> |
\bar{Q}^2 = 2 \left(3 Q : Q - (\text{Tr} Q)^2 \right) |
1481 |
> |
\end{equation} |
1482 |
> |
for the primitive quadrupole as defined in Eq. \ref{eq:quadrupole}. |
1483 |
> |
(For the traceless quadrupole tensor, $\Theta = 3 Q - \text{Tr} Q$, |
1484 |
> |
the effective moment, $\bar{Q}^2 = \frac{2}{3} \Theta : \Theta$.) |
1485 |
> |
|
1486 |
> |
To test the new electrostatic methods for quadrupoles, we have |
1487 |
> |
constructed very large, $N=$ 8,000~(sc), 16,000~(bcc), and |
1488 |
> |
32,000~(fcc) arrays of linear quadrupoles in the orientations |
1489 |
> |
described in Ref. \onlinecite{Nagai01081960}. We have compared the |
1490 |
> |
energy constants for the lowest energy configurations for these linear |
1491 |
> |
quadrupoles. Convergence to these constants are shown as a function |
1492 |
> |
of both the cutoff radius, $r_c$, and the damping parameter, $\alpha$ |
1493 |
> |
in Figs.~\ref{fig:Quadrupoles_rCut} and \ref{fig:Quadrupoles_alpha}. |
1494 |
> |
|
1495 |
> |
\begin{figure} |
1496 |
> |
\includegraphics[width=\linewidth]{Quadrupoles_rcutConvergence.pdf} |
1497 |
> |
\caption{Convergence of the lattice energy constants as a function of |
1498 |
> |
cutoff radius (normalized by the lattice constant, $a$) for the new |
1499 |
> |
real-space methods. Three quadrupolar crystal structures were |
1500 |
> |
sampled, and the analytic energy constants for the three lattices |
1501 |
> |
are indicated with grey dashed lines. The left panel shows results |
1502 |
> |
for the undamped kernel ($1/r$), while the damped error function |
1503 |
> |
kernel, $B_0(r)$ was used in the right panel.} |
1504 |
> |
\label{fig:Quadrupoles_rCut} |
1505 |
> |
\end{figure} |
1506 |
> |
|
1507 |
> |
|
1508 |
> |
\begin{figure}[!htbp] |
1509 |
> |
\includegraphics[width=\linewidth]{Quadrupoles_newAlpha.pdf} |
1510 |
> |
\caption{Convergence to the lattice energy constants as a function of |
1511 |
> |
the reduced damping parameter ($\alpha^* = \alpha a$) for the |
1512 |
> |
different real-space methods in the same three quadrupolar crystals |
1513 |
> |
in Figure \ref{fig:Quadrupoles_rCut}. The left panel shows |
1514 |
> |
results for a relatively small cutoff radius ($r_c = 4.5 a$) while a |
1515 |
> |
larger cutoff radius ($r_c = 6 a$) was used in the right panel. } |
1516 |
> |
\label{fig:Quadrupoles_alpha} |
1517 |
> |
\end{figure} |
1518 |
|
|
1519 |
+ |
Again, we find that the hard cutoff exhibits oscillations around the |
1520 |
+ |
analytic energy constants. The shifted potential (SP) and |
1521 |
+ |
gradient-shifted force (GSF) approximations converge to the correct |
1522 |
+ |
energy smoothly by $r_c = 4 a$ even for the undamped case. The |
1523 |
+ |
Taylor-shifted force (TSF) approximation again appears to perturb the |
1524 |
+ |
potential too much inside the cutoff region to provide accurate |
1525 |
+ |
measures of the energy constants. |
1526 |
+ |
|
1527 |
+ |
|
1528 |
|
\section{Conclusion} |
1529 |
< |
We have presented two efficient real-space methods for computing the |
1529 |
> |
We have presented three efficient real-space methods for computing the |
1530 |
|
interactions between point multipoles. These methods have the benefit |
1531 |
|
of smoothly truncating the energies, forces, and torques at the cutoff |
1532 |
|
radius, making them attractive for both molecular dynamics (MD) and |
1534 |
|
(GSF) and the Shifted-Potential (SP) methods converge rapidly to the |
1535 |
|
correct lattice energies for ordered dipolar and quadrupolar arrays, |
1536 |
|
while the Taylor-Shifted Force (TSF) is too severe an approximation to |
1537 |
< |
provide accurate convergence to lattice energies. |
1537 |
> |
provide accurate convergence to lattice energies. |
1538 |
|
|
1539 |
|
In most cases, GSF can obtain nearly quantitative agreement with the |
1540 |
|
lattice energy constants with reasonably small cutoff radii. The only |
1553 |
|
\begin{acknowledgments} |
1554 |
|
JDG acknowledges helpful discussions with Christopher |
1555 |
|
Fennell. Support for this project was provided by the National |
1556 |
< |
Science Foundation under grant CHE-0848243. Computational time was |
1556 |
> |
Science Foundation under grant CHE-1362211. Computational time was |
1557 |
|
provided by the Center for Research Computing (CRC) at the |
1558 |
|
University of Notre Dame. |
1559 |
|
\end{acknowledgments} |
1649 |
|
\begin{equation} |
1650 |
|
u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) . |
1651 |
|
\end{equation} |
1652 |
< |
|
1652 |
> |
% The functions |
1653 |
> |
% needed are listed schematically below: |
1654 |
> |
% % |
1655 |
> |
% \begin{eqnarray} |
1656 |
> |
% f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ |
1657 |
> |
% g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ |
1658 |
> |
% h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ |
1659 |
> |
% s_2 \quad s_3 \quad &s_4 \nonumber \\ |
1660 |
> |
% t_3 \quad &t_4 \nonumber \\ |
1661 |
> |
% &u_4 \nonumber . |
1662 |
> |
% \end{eqnarray} |
1663 |
|
The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and |
1664 |
< |
stored on a grid for values of $r$ from $0$ to $r_c$. The functions |
1665 |
< |
needed are listed schematically below: |
1664 |
> |
stored on a grid for values of $r$ from $0$ to $r_c$. Using these |
1665 |
> |
functions, we find |
1666 |
|
% |
1617 |
– |
\begin{eqnarray} |
1618 |
– |
f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ |
1619 |
– |
g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ |
1620 |
– |
h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ |
1621 |
– |
s_2 \quad s_3 \quad &s_4 \nonumber \\ |
1622 |
– |
t_3 \quad &t_4 \nonumber \\ |
1623 |
– |
&u_4 \nonumber . |
1624 |
– |
\end{eqnarray} |
1625 |
– |
|
1626 |
– |
Using these functions, we find |
1627 |
– |
% |
1667 |
|
\begin{align} |
1668 |
|
\frac{\partial f_n}{\partial r_\alpha} =&r_\alpha \frac {g_n}{r} \label{eq:b9}\\ |
1669 |
|
\frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =&\delta_{\alpha \beta}\frac {g_n}{r} |
1670 |
|
+r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) \\ |
1671 |
< |
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} =& |
1671 |
> |
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial r_\gamma} =& |
1672 |
|
\left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + |
1673 |
|
\delta_{ \beta \gamma} r_\alpha \right) |
1674 |
< |
\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) |
1675 |
< |
+ r_\alpha r_\beta r_\gamma |
1674 |
> |
\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\ |
1675 |
> |
& + r_\alpha r_\beta r_\gamma |
1676 |
|
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\ |
1677 |
< |
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} =& |
1677 |
> |
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial |
1678 |
> |
r_\gamma \partial r_\delta} =& |
1679 |
|
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} |
1680 |
|
+ \delta_{\alpha \gamma} \delta_{\beta \delta} |
1681 |
|
+\delta_{ \beta \gamma} \delta_{\alpha \delta} \right) |
1688 |
|
\left( -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5} |
1689 |
|
+ \frac{t_n}{r^4} \right)\\ |
1690 |
|
\frac{\partial^5 f_n} |
1691 |
< |
{\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} =& |
1691 |
> |
{\partial r_\alpha \partial r_\beta \partial r_\gamma \partial |
1692 |
> |
r_\delta \partial r_\epsilon} =& |
1693 |
|
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon |
1694 |
|
+ \text{14 permutations} \right) |
1695 |
|
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\ |
1711 |
|
rather the individual terms in the multipole interaction energies. |
1712 |
|
For damped charges , this still brings into the algebra multiple |
1713 |
|
derivatives of the Smith's $B_0(r)$ function. To denote these terms, |
1714 |
< |
we generalize the notation of the previous appendix. For $f(r)=1/r$ |
1715 |
< |
(bare Coulomb) or $f(r)=B_0(r)$ (smeared charge) |
1714 |
> |
we generalize the notation of the previous appendix. For either |
1715 |
> |
$f(r)=1/r$ (undamped) or $f(r)=B_0(r)$ (damped), |
1716 |
|
% |
1717 |
|
\begin{align} |
1718 |
|
g(r)=& \frac{df}{d r}\\ |
1722 |
|
u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} . |
1723 |
|
\end{align} |
1724 |
|
% |
1725 |
< |
For undamped charges, $f(r)=1/r$, Table I lists these derivatives |
1726 |
< |
under the column ``Bare Coulomb.'' Equations \ref{eq:b9} to |
1727 |
< |
\ref{eq:b13} are still correct for GSF electrostatics if the subscript |
1687 |
< |
$n$ is eliminated. |
1725 |
> |
For undamped charges Table I lists these derivatives under the column |
1726 |
> |
``Bare Coulomb.'' Equations \ref{eq:b9} to \ref{eq:b13} are still |
1727 |
> |
correct for GSF electrostatics if the subscript $n$ is eliminated. |
1728 |
|
|
1729 |
|
\newpage |
1730 |
|
|