| 65 |
|
|
| 66 |
|
\begin{abstract} |
| 67 |
|
We have extended the original damped-shifted force (DSF) |
| 68 |
< |
electrostatic kernel and have been able to derive two new |
| 68 |
> |
electrostatic kernel and have been able to derive three new |
| 69 |
|
electrostatic potentials for higher-order multipoles that are based |
| 70 |
< |
on truncated Taylor expansions around the cutoff radius. For |
| 71 |
< |
multipole-multipole interactions, we find that each of the distinct |
| 72 |
< |
orientational contributions has a separate radial function to ensure |
| 73 |
< |
that the overall forces and torques vanish at the cutoff radius. In |
| 74 |
< |
this paper, we present energy, force, and torque expressions for the |
| 75 |
< |
new models, and compare these real-space interaction models to exact |
| 76 |
< |
results for ordered arrays of multipoles. |
| 70 |
> |
on truncated Taylor expansions around the cutoff radius. These |
| 71 |
> |
include a shifted potential (SP) that generalizes the Wolf method |
| 72 |
> |
for point multipoles, and Taylor-shifted force (TSF) and |
| 73 |
> |
gradient-shifted force (GSF) potentials that are both |
| 74 |
> |
generalizations of DSF electrostatics for multipoles. We find that |
| 75 |
> |
each of the distinct orientational contributions requires a separate |
| 76 |
> |
radial function to ensure that pairwise energies, forces and torques |
| 77 |
> |
all vanish at the cutoff radius. In this paper, we present energy, |
| 78 |
> |
force, and torque expressions for the new models, and compare these |
| 79 |
> |
real-space interaction models to exact results for ordered arrays of |
| 80 |
> |
multipoles. We find that the GSF and SP methods converge rapidly to |
| 81 |
> |
the correct lattice energies for ordered dipolar and quadrupolar |
| 82 |
> |
arrays, while the Taylor-Shifted Force (TSF) is too severe an |
| 83 |
> |
approximation to provide accurate convergence to lattice energies. |
| 84 |
> |
Because real-space methods can be made to scale linearly with system |
| 85 |
> |
size, the SP and GSF are attractive options for large Monte Carlo |
| 86 |
> |
and molecular dynamics simulations. |
| 87 |
|
\end{abstract} |
| 88 |
|
|
| 89 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 302 |
|
section \ref{sec:damped}. |
| 303 |
|
|
| 304 |
|
It is convenient to locate charges $q_j$ relative to the center of mass of $\bf b$. Then with $\bf{r}$ pointing from |
| 305 |
< |
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by |
| 305 |
> |
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $), the interaction energy is given by |
| 306 |
|
\begin{equation} |
| 307 |
|
U_{\bf{ab}}(r) |
| 308 |
|
= \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} . |
| 320 |
|
\end{equation} |
| 321 |
|
This form has the benefit of separating out the energies of |
| 322 |
|
interaction into contributions from the charge, dipole, and quadrupole |
| 323 |
< |
of $\bf a$ interacting with the same multipoles on $\bf b$. |
| 323 |
> |
of $\bf a$ interacting with the same multipoles in $\bf b$. |
| 324 |
|
|
| 325 |
|
\subsection{Damped Coulomb interactions} |
| 326 |
|
\label{sec:damped} |
| 339 |
|
either $1/r$ or $B_0(r)$, and all of the techniques can be applied to |
| 340 |
|
bare or damped Coulomb kernels (or any other function) as long as |
| 341 |
|
derivatives of these functions are known. Smith's convenient |
| 342 |
< |
functions $B_l(r)$ are summarized in Appendix A. (N.B. there is one |
| 343 |
< |
important distinction between the two kernels, which is the behavior |
| 344 |
< |
of $\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$. The former |
| 345 |
< |
is zero everywhere except for a delta function evaluated at the |
| 346 |
< |
origin. The latter also has delta function behavior, but is non-zero |
| 347 |
< |
for $r \neq 0$. Thus the standard justification for using a traceless |
| 342 |
> |
functions $B_l(r)$, which are used for derivatives of the damped |
| 343 |
> |
kernel, are summarized in Appendix A. (N.B. there is one important |
| 344 |
> |
distinction between the two kernels, which is the behavior of |
| 345 |
> |
$\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$. The former is |
| 346 |
> |
zero everywhere except for a delta function evaluated at the origin. |
| 347 |
> |
The latter also has delta function behavior, but is non-zero for $r |
| 348 |
> |
\neq 0$. Thus the standard justification for using a traceless |
| 349 |
|
quadrupole tensor fails for the damped case.) |
| 350 |
|
|
| 351 |
|
The main goal of this work is to smoothly cut off the interaction |
| 504 |
|
approaches below -- for each separate orientational contribution, only |
| 505 |
|
the radial factors differ between the two methods. |
| 506 |
|
|
| 507 |
+ |
\subsection{Generalization of the Wolf shifted potential (SP)} |
| 508 |
+ |
It is also possible to formulate an extension of the Wolf approach for |
| 509 |
+ |
multipoles by simply projecting the image multipole onto the surface |
| 510 |
+ |
of the cutoff sphere, and including the interactions with the central |
| 511 |
+ |
multipole and the image. This effectively shifts the pair potential |
| 512 |
+ |
to zero at the cutoff radius, |
| 513 |
+ |
\begin{equation} |
| 514 |
+ |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 515 |
+ |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 516 |
+ |
\label{eq:SP} |
| 517 |
+ |
\end{equation} |
| 518 |
+ |
independent of the orientations of the two multipoles. The sum again |
| 519 |
+ |
describes separate potential shifting that is applied to each |
| 520 |
+ |
orientational contribution to the energy. |
| 521 |
+ |
|
| 522 |
+ |
The shifted potential (SP) method is a simple truncation of the GSF |
| 523 |
+ |
method for each orientational contribution, leaving out the $(r-r_c)$ |
| 524 |
+ |
terms that multiply the gradient. Functional forms for the |
| 525 |
+ |
shifted-potential (SP) method can also be summarized using the form of |
| 526 |
+ |
Eq.~(\ref{eq:generic}). The energy, force, and torque expressions are |
| 527 |
+ |
tabulated below for all three methods. As in the GSF and TSF methods, |
| 528 |
+ |
for each separate orientational contribution, only the radial factors |
| 529 |
+ |
differ between the SP, GSF, and TSF methods. |
| 530 |
+ |
|
| 531 |
+ |
|
| 532 |
|
\subsection{\label{sec:level2}Body and space axes} |
| 533 |
|
Although objects $\bf a$ and $\bf b$ rotate during a molecular |
| 534 |
|
dynamics (MD) simulation, their multipole tensors remain fixed in |
| 537 |
|
in intermediate forms involving the vectors of the rotation matrices. |
| 538 |
|
We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors |
| 539 |
|
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$. |
| 540 |
< |
In a typical simulation , the initial axes are obtained by |
| 540 |
> |
In a typical simulation, the initial axes are obtained by |
| 541 |
|
diagonalizing the moment of inertia tensors for the objects. (N.B., |
| 542 |
|
the body axes are generally {\it not} the same as those for which the |
| 543 |
|
quadrupole moment is diagonal.) The rotation matrices are then |
| 609 |
|
(\hat{a}_m \times \hat{b}_n) . |
| 610 |
|
\end{eqnarray} |
| 611 |
|
|
| 612 |
< |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ |
| 612 |
> |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $ |
| 613 |
|
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} |
| 614 |
|
We also made use of the identities, |
| 615 |
|
% |
| 619 |
|
\right) \\ |
| 620 |
|
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
| 621 |
|
=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} |
| 622 |
< |
\right) . |
| 622 |
> |
\right). |
| 623 |
|
\end{align} |
| 624 |
|
|
| 625 |
|
Many of the multipole contractions required can be written in one of |
| 631 |
|
\begin{equation} |
| 632 |
|
\mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} |
| 633 |
|
= D_{\bf {a}\alpha} D_{\bf {b}\alpha} = |
| 634 |
< |
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} |
| 634 |
> |
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}. |
| 635 |
|
\end{equation} |
| 636 |
|
% |
| 637 |
|
The first two forms are written using space coordinates. The first |
| 662 |
|
the cutoff sphere. For a system of undamped charges, the total |
| 663 |
|
self-term is |
| 664 |
|
\begin{equation} |
| 665 |
< |
U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2} |
| 665 |
> |
U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}. |
| 666 |
|
\label{eq:selfTerm} |
| 667 |
|
\end{equation} |
| 668 |
|
|
| 695 |
|
cutoff sphere, and averaging over the surface of the cutoff sphere. |
| 696 |
|
Because the self term is carried out as a single sum over sites, the |
| 697 |
|
reciprocal-space portion is identical to half of the self-term |
| 698 |
< |
obtained by Smith and Aguado and Madden for the application of the |
| 699 |
< |
Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given |
| 700 |
< |
site which posesses a charge, dipole, and quadrupole, both types of |
| 701 |
< |
contribution are given in table \ref{tab:tableSelf}. |
| 698 |
> |
obtained by Smith, and also by Aguado and Madden for the application |
| 699 |
> |
of the Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a |
| 700 |
> |
given site which posesses a charge, dipole, and quadrupole, both types |
| 701 |
> |
of contribution are given in table \ref{tab:tableSelf}. |
| 702 |
|
|
| 703 |
|
\begin{table*} |
| 704 |
|
\caption{\label{tab:tableSelf} Self-interaction contributions for |
| 873 |
|
\begin{sidewaystable} |
| 874 |
|
\caption{\label{tab:tableenergy}Radial functions used in the energy |
| 875 |
|
and torque equations. The $f, g, h, s, t, \mathrm{and} u$ |
| 876 |
< |
functions used in this table are defined in Appendices B and C.} |
| 876 |
> |
functions used in this table are defined in Appendices B and C. |
| 877 |
> |
Gradient shifted (GSF) functions are constructed using the shifted |
| 878 |
> |
potential (SP) functions.} |
| 879 |
|
\begin{tabular}{|c|c|l|l|l|} \hline |
| 880 |
|
Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF) |
| 881 |
|
\\ \hline |
| 906 |
|
$-\frac{1}{r^3} $ & |
| 907 |
|
$\frac{g_2(r)}{r} $ & |
| 908 |
|
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c}$ & |
| 909 |
< |
SP $-(r-r_c) |
| 910 |
< |
\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
| 909 |
> |
SP $-(r-r_c) \left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
| 910 |
> |
% |
| 911 |
> |
% |
| 912 |
> |
% |
| 913 |
|
$v_{22}(r)$ & |
| 914 |
|
$\frac{3}{r^3} $ & |
| 915 |
|
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ & |
| 916 |
< |
$\left(-\frac{g(r)}{r}+h(r) \right)$ & SP \\ |
| 917 |
< |
&&& $~~~-\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)$\\ |
| 916 |
> |
$\left(-\frac{g(r)}{r}+h(r) \right) -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ |
| 917 |
> |
& SP $-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\ |
| 918 |
|
% |
| 919 |
|
% |
| 920 |
|
% |
| 922 |
|
$v_{31}(r)$ & |
| 923 |
|
$\frac{3}{r^4} $ & |
| 924 |
|
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ & |
| 925 |
< |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)$ & SP \\ |
| 926 |
< |
&&& $-\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)$ \\ |
| 925 |
> |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r}\right)-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right)$ |
| 926 |
> |
& SP $-(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)$ \\ |
| 927 |
|
% |
| 928 |
+ |
% |
| 929 |
+ |
% |
| 930 |
|
$v_{32}(r)$ & |
| 931 |
|
$-\frac{15}{r^4} $ & |
| 932 |
|
$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ & |
| 933 |
< |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$ & SP \\ |
| 934 |
< |
&&& $~~~- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) |
| 935 |
< |
\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)$ \\ |
| 933 |
> |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$& |
| 934 |
> |
SP $-(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}\right.$ \\ |
| 935 |
> |
&&& $~~~-\left(\frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c)\right)$ & |
| 936 |
> |
$\phantom{SP-(r-r_c)}\left.-\frac{3s(r_c)}{r_c}+t(r_c) \right)$\\ |
| 937 |
|
% |
| 938 |
|
% |
| 939 |
|
% |
| 941 |
|
$v_{41}(r)$ & |
| 942 |
|
$\frac{3}{r^5} $ & |
| 943 |
|
$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $ & |
| 944 |
< |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)$ & SP \\ |
| 945 |
< |
&&& $~~~- \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)$ |
| 944 |
> |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)- \left(-\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ & |
| 945 |
> |
SP $-(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)$ |
| 946 |
|
\\ |
| 947 |
|
% 2 |
| 948 |
|
$v_{42}(r)$ & |
| 949 |
|
$- \frac{15}{r^5} $ & |
| 950 |
|
$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ & |
| 951 |
< |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} |
| 952 |
< |
\right)$ & SP \\ |
| 953 |
< |
&&& $~~~-\left( \frac{3g(r_c)}{r_c^3} - |
| 954 |
< |
\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)$ \\ |
| 951 |
> |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)$ & |
| 952 |
> |
SP$-(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}\right.$ \\ |
| 953 |
> |
&&& $~~~-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ & |
| 954 |
> |
$\phantom{SP-(r-r_c)}\left. -\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$\\ |
| 955 |
|
% 3 |
| 956 |
+ |
% |
| 957 |
+ |
% |
| 958 |
|
$v_{43}(r)$ & |
| 959 |
|
$ \frac{105}{r^5} $ & |
| 960 |
|
$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ & |
| 961 |
< |
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + |
| 962 |
< |
t(r)\right)$ & SP $-(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}\right.$\\ |
| 963 |
< |
&&& |
| 964 |
< |
$~~~-\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)$ \\ |
| 961 |
> |
$ \left(-\frac{15g(r)}{r^3} +\frac{15h(r)}{r^2}-\frac{6s(r)}{r}+t(r)\right) $ & |
| 962 |
> |
SP $-(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}\right.$\\ |
| 963 |
> |
&&& $~~~-\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)$ & |
| 964 |
> |
$\phantom{SP-(r-r_c)}\left.+\frac{21s(r_c)}{r_c^2}-\frac{6t(r_c)}{r_c}+u(r_c) \right)$\\ |
| 965 |
|
\hline |
| 966 |
|
\end{tabular} |
| 967 |
|
\end{sidewaystable} |
| 971 |
|
% |
| 972 |
|
|
| 973 |
|
\begin{sidewaystable} |
| 974 |
< |
\caption{\label{tab:tableFORCE}Radial functions used in the force equations.} |
| 975 |
< |
\begin{tabular}{|c|c|l|l|} \hline |
| 976 |
< |
Function&Definition&Taylor-Shifted&Gradient-Shifted |
| 974 |
> |
\caption{\label{tab:tableFORCE}Radial functions used in the force |
| 975 |
> |
equations. Gradient shifted (GSF) functions are constructed using the shifted |
| 976 |
> |
potential (SP) functions. Some of these functions are simple |
| 977 |
> |
modifications of the functions found in table \ref{tab:tableenergy}} |
| 978 |
> |
\begin{tabular}{|c|c|l|l|l|} \hline |
| 979 |
> |
Function&Definition&Taylor-Shifted (TSF)& Shifted Potential (SP) |
| 980 |
> |
&Gradient-Shifted (GSF) |
| 981 |
|
\\ \hline |
| 982 |
|
% |
| 983 |
|
% |
| 985 |
|
$w_a(r)$& |
| 986 |
|
$\frac{d v_{01}}{dr}$& |
| 987 |
|
$g_0(r)$& |
| 988 |
< |
$g(r)-g(r_c)$ \\ |
| 988 |
> |
$g(r)$& |
| 989 |
> |
SP $-g(r_c)$ \\ |
| 990 |
|
% |
| 991 |
|
% |
| 992 |
|
$w_b(r)$ & |
| 993 |
|
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $& |
| 994 |
|
$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ & |
| 995 |
< |
$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\ |
| 995 |
> |
$h(r) - \frac{v_{11}(r)}{r} $ & |
| 996 |
> |
SP $- h(r_c)$ \\ |
| 997 |
|
% |
| 998 |
|
$w_c(r)$ & |
| 999 |
|
$\frac{v_{11}(r)}{r}$ & |
| 1000 |
|
$\frac{g_1(r)}{r} $ & |
| 1001 |
+ |
$\frac{v_{11}(r)}{r}$& |
| 1002 |
|
$\frac{v_{11}(r)}{r}$\\ |
| 1003 |
|
% |
| 1004 |
|
% |
| 1005 |
|
$w_d(r)$& |
| 1006 |
|
$\frac{d v_{21}}{dr}$& |
| 1007 |
|
$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ & |
| 1008 |
< |
$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right) |
| 1009 |
< |
-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\ |
| 1008 |
> |
$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)$ & |
| 1009 |
> |
SP $-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\ |
| 1010 |
|
% |
| 1011 |
|
$w_e(r)$ & |
| 1012 |
|
$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ & |
| 1013 |
|
$\frac{v_{22}(r)}{r}$ & |
| 1014 |
+ |
$\frac{v_{22}(r)}{r}$ & |
| 1015 |
|
$\frac{v_{22}(r)}{r}$ \\ |
| 1016 |
|
% |
| 1017 |
|
% |
| 1018 |
|
$w_f(r)$& |
| 1019 |
|
$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$& |
| 1020 |
|
$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ & |
| 1021 |
< |
$ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) $ \\ |
| 1022 |
< |
&&& $ ~~~- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) |
| 973 |
< |
\right)-\frac{2v_{22}(r)}{r}$\\ |
| 1021 |
> |
$ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) -\frac{2v_{22}(r)}{r}$& |
| 1022 |
> |
SP $- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\ |
| 1023 |
|
% |
| 1024 |
|
$w_g(r)$& |
| 1025 |
|
$\frac{v_{31}(r)}{r}$& |
| 1026 |
|
$ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& |
| 1027 |
+ |
$\frac{v_{31}(r)}{r}$& |
| 1028 |
|
$\frac{v_{31}(r)}{r}$\\ |
| 1029 |
|
% |
| 1030 |
|
$w_h(r)$ & |
| 1031 |
|
$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$& |
| 1032 |
|
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $ & |
| 1033 |
< |
$ \left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - \left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\ |
| 1034 |
< |
&&& $ ~~~ -\frac{v_{31}(r)}{r}$ \\ |
| 1033 |
> |
$ \left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) -\frac{v_{31}(r)}{r}$ & |
| 1034 |
> |
SP $ - \left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\ |
| 1035 |
|
% 2 |
| 1036 |
|
$w_i(r)$ & |
| 1037 |
|
$\frac{v_{32}(r)}{r}$ & |
| 1038 |
|
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $ & |
| 1039 |
+ |
$\frac{v_{32}(r)}{r}$& |
| 1040 |
|
$\frac{v_{32}(r)}{r}$\\ |
| 1041 |
|
% |
| 1042 |
|
$w_j(r)$ & |
| 1043 |
|
$\frac{d v_{32}}{dr} - \frac{3v_{32}}{r}$& |
| 1044 |
|
$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right) $ & |
| 1045 |
< |
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right)$ \\ |
| 1046 |
< |
&&& $~~~-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2} |
| 1047 |
< |
-\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\ |
| 1045 |
> |
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) -\frac{3v_{32}}{r}$ & |
| 1046 |
> |
SP $-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2} |
| 1047 |
> |
-\frac{3s(r_c)}{r_c} +t(r_c) \right)$ \\ |
| 1048 |
|
% |
| 1049 |
|
$w_k(r)$ & |
| 1050 |
|
$\frac{d v_{41}}{dr} $ & |
| 1051 |
|
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2} \right)$ & |
| 1052 |
< |
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2} \right) |
| 1053 |
< |
-\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2} \right)$ \\ |
| 1052 |
> |
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2} |
| 1053 |
> |
\right)$ & |
| 1054 |
> |
SP $-\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2} \right)$ \\ |
| 1055 |
|
% |
| 1056 |
|
$w_l(r)$ & |
| 1057 |
|
$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ & |
| 1058 |
|
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ & |
| 1059 |
< |
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\ |
| 1060 |
< |
&&& $~~~ -\left(-\frac{9g(r_c)}{r_c^4} +\frac{9h(r_c)}{r_c^3} -\frac{4s(r_c)}{r_c^2} +\frac{t(r_c)}{r_c} \right) |
| 1061 |
< |
-\frac{2v_{42}(r)}{r}$\\ |
| 1059 |
> |
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} |
| 1060 |
> |
+\frac{t(r)}{r} \right) -\frac{2v_{42}(r)}{r}$& |
| 1061 |
> |
SP$-\left(-\frac{9g(r_c)}{r_c^4} +\frac{9h(r_c)}{r_c^3} -\frac{4s(r_c)}{r_c^2} +\frac{t(r_c)}{r_c} \right)$\\ |
| 1062 |
|
% |
| 1063 |
|
$w_m(r)$ & |
| 1064 |
|
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$& |
| 1065 |
< |
$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} + \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$ & |
| 1066 |
< |
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$\\ |
| 1067 |
< |
&&& $~~~- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3} |
| 1068 |
< |
+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\ |
| 1069 |
< |
&&& $~~~-\frac{4v_{43}(r)}{r}$ \\ |
| 1065 |
> |
$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} \right.$ & |
| 1066 |
> |
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2}\right.$ & |
| 1067 |
> |
SP $- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}\right.$ \\ |
| 1068 |
> |
&& $~~~\left.+ \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$ |
| 1069 |
> |
& $~~~\left. -\frac{6t(r)}{r} +u(r) \right) -\frac{4v_{43}(r)}{r}$ & |
| 1070 |
> |
$\phantom{SP-} \left.+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\ |
| 1071 |
|
% |
| 1072 |
|
$w_n(r)$ & |
| 1073 |
|
$\frac{v_{42}(r)}{r}$ & |
| 1074 |
|
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2} \right)$ & |
| 1075 |
+ |
$\frac{v_{42}(r)}{r}$& |
| 1076 |
|
$\frac{v_{42}(r)}{r}$\\ |
| 1077 |
|
% |
| 1078 |
|
$w_o(r)$ & |
| 1079 |
|
$\frac{v_{43}(r)}{r}$& |
| 1080 |
|
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ & |
| 1081 |
+ |
$\frac{v_{43}(r)}{r}$& |
| 1082 |
|
$\frac{v_{43}(r)}{r}$ \\ \hline |
| 1083 |
|
% |
| 1084 |
|
|
| 1399 |
|
\left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta} |
| 1400 |
|
-\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta} |
| 1401 |
|
\right] |
| 1402 |
+ |
\label{eq:matrixCross} |
| 1403 |
|
\end{equation} |
| 1404 |
|
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic |
| 1405 |
|
permuations of the matrix indices. |
| 1411 |
|
$\mathbf{a}$ interacting with those of lower order on site |
| 1412 |
|
$\mathbf{b}$ can be obtained by swapping indices in the expressions |
| 1413 |
|
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. |
| 1414 |
|
|
| 1415 |
|
\section{Comparison to known multipolar energies} |
| 1416 |
|
|
| 1443 |
|
dipoles are contained in a plane perpendicular to the dipole direction |
| 1444 |
|
that passes through the dipole. We have also studied the minimum |
| 1445 |
|
energy structure for the BCC lattice that was found by Luttinger \& |
| 1446 |
< |
Tisza. The total electrostatic energy for any of the arrays is given |
| 1447 |
< |
by: |
| 1446 |
> |
Tisza. The total electrostatic energy density for any of the arrays |
| 1447 |
> |
is given by: |
| 1448 |
|
\begin{equation} |
| 1449 |
|
E = C N^2 \mu^2 |
| 1450 |
|
\end{equation} |
| 1490 |
|
|
| 1491 |
|
The hard cutoff exhibits oscillations around the analytic energy |
| 1492 |
|
constants, and converges to incorrect energies when the complementary |
| 1493 |
< |
error function damping kernel is used. The shifted potential (SP) and |
| 1494 |
< |
gradient-shifted force (GSF) approximations converge to the correct |
| 1495 |
< |
energy smoothly by $r_c = 6 a$ even for the undamped case. This |
| 1496 |
< |
indicates that the correction provided by the self term is required |
| 1497 |
< |
for obtaining accurate energies. The Taylor-shifted force (TSF) |
| 1498 |
< |
approximation appears to perturb the potential too much inside the |
| 1499 |
< |
cutoff region to provide accurate measures of the energy constants. |
| 1493 |
> |
error function damping kernel is used. The shifted potential (SP) |
| 1494 |
> |
converges to the correct energy smoothly by $r_c = 4.5 a$ even for the |
| 1495 |
> |
undamped case. This indicates that the shifting and the correction |
| 1496 |
> |
provided by the self term are required for obtaining accurate energies. |
| 1497 |
> |
The Taylor-shifted force (TSF) approximation appears to perturb the |
| 1498 |
> |
potential too much inside the cutoff region to provide accurate |
| 1499 |
> |
measures of the energy constants. GSF is a compromise, converging to |
| 1500 |
> |
the correct energies within $r_c = 6 a$. |
| 1501 |
|
|
| 1502 |
|
{\it Quadrupolar} analogues to the Madelung constants were first |
| 1503 |
|
worked out by Nagai and Nakamura who computed the energies of selected |
| 1522 |
|
constructed very large, $N=$ 8,000~(sc), 16,000~(bcc), and |
| 1523 |
|
32,000~(fcc) arrays of linear quadrupoles in the orientations |
| 1524 |
|
described in Ref. \onlinecite{Nagai01081960}. We have compared the |
| 1525 |
< |
energy constants for the lowest energy configurations for these linear |
| 1526 |
< |
quadrupoles. Convergence to these constants are shown as a function |
| 1527 |
< |
of both the cutoff radius, $r_c$, and the damping parameter, $\alpha$ |
| 1528 |
< |
in Figs.~\ref{fig:Quadrupoles_rCut} and \ref{fig:Quadrupoles_alpha}. |
| 1525 |
> |
energy constants for these low-energy configurations for linear |
| 1526 |
> |
quadrupoles. Convergence to these constants are shown as a function of |
| 1527 |
> |
both the cutoff radius, $r_c$, and the damping parameter, $\alpha$ in |
| 1528 |
> |
Figs.~\ref{fig:Quadrupoles_rCut} and \ref{fig:Quadrupoles_alpha}. |
| 1529 |
|
|
| 1530 |
|
\begin{figure} |
| 1531 |
|
\includegraphics[width=\linewidth]{Quadrupoles_rcutConvergence.pdf} |
| 1552 |
|
\end{figure} |
| 1553 |
|
|
| 1554 |
|
Again, we find that the hard cutoff exhibits oscillations around the |
| 1555 |
< |
analytic energy constants. The shifted potential (SP) and |
| 1556 |
< |
gradient-shifted force (GSF) approximations converge to the correct |
| 1557 |
< |
energy smoothly by $r_c = 4 a$ even for the undamped case. The |
| 1558 |
< |
Taylor-shifted force (TSF) approximation again appears to perturb the |
| 1559 |
< |
potential too much inside the cutoff region to provide accurate |
| 1560 |
< |
measures of the energy constants. |
| 1555 |
> |
analytic energy constants. The shifted potential (SP) approximation |
| 1556 |
> |
converges to the correct energy smoothly by $r_c = 3 a$ even for the |
| 1557 |
> |
undamped case. The Taylor-shifted force (TSF) approximation again |
| 1558 |
> |
appears to perturb the potential too much inside the cutoff region to |
| 1559 |
> |
provide accurate measures of the energy constants. GSF again provides |
| 1560 |
> |
a compromise between the two methods -- energies are converged by $r_c |
| 1561 |
> |
= 4.5 a$, and the approximation is not as perturbative at short range |
| 1562 |
> |
as TSF. |
| 1563 |
|
|
| 1564 |
+ |
It is also useful to understand the convergence to the lattice energy |
| 1565 |
+ |
constants as a function of the reduced damping parameter ($\alpha^* = |
| 1566 |
+ |
\alpha a$) for the different real-space methods. |
| 1567 |
+ |
Figures. \ref{fig:Dipoles_alpha} and \ref{fig:Quadrupoles_alpha} show |
| 1568 |
+ |
this comparison for the dipolar and quadrupolar lattices, |
| 1569 |
+ |
respectively. All of the methods (except for TSF) have excellent |
| 1570 |
+ |
behavior for the undamped or weakly-damped cases. All of the methods |
| 1571 |
+ |
can be forced to converge by increasing the value of $\alpha$, which |
| 1572 |
+ |
effectively shortens the overall range of the potential, but which |
| 1573 |
+ |
equalizing the truncation effects on the different orientational |
| 1574 |
+ |
contributions. In the second paper in the series, we discuss how |
| 1575 |
+ |
large values of $\alpha$ can perturb the force and torque vectors, but |
| 1576 |
+ |
both weakly-damped or over-damped electrostatics appear to generate |
| 1577 |
+ |
reasonable values for the total electrostatic energies under both the |
| 1578 |
+ |
SP and GSF approximations. |
| 1579 |
|
|
| 1580 |
|
\section{Conclusion} |
| 1581 |
|
We have presented three efficient real-space methods for computing the |
| 1598 |
|
crystals with net-zero moments, so this is not expected to be an issue |
| 1599 |
|
in most simulations. |
| 1600 |
|
|
| 1601 |
+ |
The techniques used here to derive the force, torque and energy |
| 1602 |
+ |
expressions can be extended to higher order multipoles, although some |
| 1603 |
+ |
of the objects (e.g. the matrix cross product in |
| 1604 |
+ |
Eq. \ref{eq:matrixCross}) will need to be generalized for higher-rank |
| 1605 |
+ |
tensors. We also note that the definitions of the multipoles used |
| 1606 |
+ |
here are in a primitive form, and these need some care when comparing |
| 1607 |
+ |
with experiment or other computational techniques. |
| 1608 |
+ |
|
| 1609 |
|
In large systems, these new methods can be made to scale approximately |
| 1610 |
|
linearly with system size, and detailed comparisons with the Ewald sum |
| 1611 |
|
for a wide range of chemical environments follows in the second paper. |
| 1769 |
|
|
| 1770 |
|
In Gradient-shifted force electrostatics, the kernel is not expanded, |
| 1771 |
|
rather the individual terms in the multipole interaction energies. |
| 1772 |
< |
For damped charges , this still brings into the algebra multiple |
| 1772 |
> |
For damped charges, this still brings into the algebra multiple |
| 1773 |
|
derivatives of the Smith's $B_0(r)$ function. To denote these terms, |
| 1774 |
|
we generalize the notation of the previous appendix. For either |
| 1775 |
|
$f(r)=1/r$ (undamped) or $f(r)=B_0(r)$ (damped), |
| 1776 |
|
% |
| 1777 |
|
\begin{align} |
| 1778 |
< |
g(r)=& \frac{df}{d r}\\ |
| 1779 |
< |
h(r)=& \frac{dg}{d r} = \frac{d^2f}{d r^2} \\ |
| 1780 |
< |
s(r)=& \frac{dh}{d r} = \frac{d^3f}{d r^3} \\ |
| 1781 |
< |
t(r)=& \frac{ds}{d r} = \frac{d^4f}{d r^4} \\ |
| 1782 |
< |
u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} . |
| 1778 |
> |
g(r) &= \frac{df}{d r} && &&=-\frac{1}{r^2} |
| 1779 |
> |
&&\mathrm{or~~~} -rB_1(r) \\ |
| 1780 |
> |
h(r) &= \frac{dg}{d r} &&= \frac{d^2f}{d r^2} &&= \frac{2}{r^3} &&\mathrm{or~~~}-B_1(r) + r^2 B_2(r) \\ |
| 1781 |
> |
s(r) &= \frac{dh}{d r} &&= \frac{d^3f}{d r^3} &&=-\frac{6}{r^4}&&\mathrm{or~~~}3rB_2(r) - r^3 B_3(r)\\ |
| 1782 |
> |
t(r) &= \frac{ds}{d r} &&= \frac{d^4f}{d r^4} &&= \frac{24}{r^5} &&\mathrm{or~~~} 3 |
| 1783 |
> |
B_2(r) - 6r^2 B_3(r) + r^4 B_4(r) \\ |
| 1784 |
> |
u(r) &= \frac{dt}{d r} &&= \frac{d^5f}{d r^5} &&=-\frac{120}{r^6} &&\mathrm{or~~~} -15 |
| 1785 |
> |
r B_3(r) + 10 r^3B_4(r) -r^5B_5(r). |
| 1786 |
|
\end{align} |
| 1787 |
|
% |
| 1788 |
< |
For undamped charges Table I lists these derivatives under the column |
| 1789 |
< |
``Bare Coulomb.'' Equations \ref{eq:b9} to \ref{eq:b13} are still |
| 1788 |
> |
For undamped charges, Table I lists these derivatives under the Bare |
| 1789 |
> |
Coulomb column. Equations \ref{eq:b9} to \ref{eq:b13} are still |
| 1790 |
|
correct for GSF electrostatics if the subscript $n$ is eliminated. |
| 1791 |
|
|
| 1792 |
|
\newpage |
