| 446 |
|
expansion, and has a similar interaction energy for all multipole |
| 447 |
|
orders: |
| 448 |
|
\begin{equation} |
| 449 |
< |
U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big |
| 450 |
< |
\lvert _{r_c} . |
| 449 |
> |
U^{\text{GSF}} = |
| 450 |
> |
U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 451 |
> |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
| 452 |
> |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} . |
| 453 |
|
\label{generic2} |
| 454 |
|
\end{equation} |
| 455 |
< |
Here the gradient for force shifting is evaluated for an image |
| 456 |
< |
multipole projected onto the surface of the cutoff sphere (see fig |
| 457 |
< |
\ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$ |
| 458 |
< |
appear. The primary difference between the TSF and GSF methods is the |
| 459 |
< |
stage at which the Taylor Series is applied; in the Taylor-shifted |
| 460 |
< |
approach, it is applied to the kernel itself. In the Gradient-shifted |
| 461 |
< |
approach, it is applied to individual radial interactions terms in the |
| 462 |
< |
multipole expansion. Energies from this method thus have the general |
| 463 |
< |
form: |
| 455 |
> |
Both the potential and the gradient for force shifting are evaluated |
| 456 |
> |
for an image multipole projected onto the surface of the cutoff sphere |
| 457 |
> |
(see fig \ref{fig:shiftedMultipoles}). The image multipole retains |
| 458 |
> |
the orientation ($\hat{\mathbf{b}}$) of the interacting multipole. No |
| 459 |
> |
higher order terms $(r-r_c)^n$ appear. The primary difference between |
| 460 |
> |
the TSF and GSF methods is the stage at which the Taylor Series is |
| 461 |
> |
applied; in the Taylor-shifted approach, it is applied to the kernel |
| 462 |
> |
itself. In the Gradient-shifted approach, it is applied to individual |
| 463 |
> |
radial interactions terms in the multipole expansion. Energies from |
| 464 |
> |
this method thus have the general form: |
| 465 |
|
\begin{equation} |
| 466 |
|
U= \sum (\text{angular factor}) (\text{radial factor}). |
| 467 |
|
\label{generic3} |
| 1336 |
|
$\mathbf{b}$ can be obtained by swapping indices in the expressions |
| 1337 |
|
above. |
| 1338 |
|
|
| 1339 |
+ |
\section{Related real-space methods} |
| 1340 |
+ |
One can also formulate a shifted potential, |
| 1341 |
+ |
\begin{equation} |
| 1342 |
+ |
U^{\text{SP}} = U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 1343 |
+ |
U(\mathbf{r}_c, \hat{\mathbf{a}}, \hat{\mathbf{b}}), |
| 1344 |
+ |
\label{eq:SP} |
| 1345 |
+ |
\end{equation} |
| 1346 |
+ |
obtained by projecting the image multipole onto the surface of the |
| 1347 |
+ |
cutoff sphere. The shifted potential (SP) can be thought of as a |
| 1348 |
+ |
simple extension to the original Wolf method. The energies and |
| 1349 |
+ |
torques for the SP can be easily obtained by zeroing out the $(r-r_c)$ |
| 1350 |
+ |
terms in the final column of table \ref{tab:tableenergy}. SP forces |
| 1351 |
+ |
(which retain discontinuities at the cutoff sphere) can be obtained by |
| 1352 |
+ |
eliminating all functions that depend on $r_c$ in the last column of |
| 1353 |
+ |
table \ref{tab:tableFORCE}. The self-energy contributions to the SP |
| 1354 |
+ |
potential are identical to both the GSF and TSF methods. |
| 1355 |
+ |
|
| 1356 |
|
\section{Comparison to known multipolar energies} |
| 1357 |
|
|
| 1358 |
|
To understand how these new real-space multipole methods behave in |
| 1359 |
|
computer simulations, it is vital to test against established methods |
| 1360 |
|
for computing electrostatic interactions in periodic systems, and to |
| 1361 |
|
evaluate the size and sources of any errors that arise from the |
| 1362 |
< |
real-space cutoffs. In this paper we test Taylor-shifted and |
| 1363 |
< |
Gradient-shifted electrostatics against analytical methods for |
| 1364 |
< |
computing the energies of ordered multipolar arrays. In the following |
| 1365 |
< |
paper, we test the new methods against the multipolar Ewald sum for |
| 1366 |
< |
computing the energies, forces and torques for a wide range of typical |
| 1367 |
< |
condensed-phase (disordered) systems. |
| 1362 |
> |
real-space cutoffs. In this paper we test both TSF and GSF |
| 1363 |
> |
electrostatics against analytical methods for computing the energies |
| 1364 |
> |
of ordered multipolar arrays. In the following paper, we test the new |
| 1365 |
> |
methods against the multipolar Ewald sum for computing the energies, |
| 1366 |
> |
forces and torques for a wide range of typical condensed-phase |
| 1367 |
> |
(disordered) systems. |
| 1368 |
|
|
| 1369 |
|
Because long-range electrostatic effects can be significant in |
| 1370 |
|
crystalline materials, ordered multipolar arrays present one of the |
| 1374 |
|
magnetization and obtained a number of these constants.\cite{Sauer} |
| 1375 |
|
This theory was developed more completely by Luttinger and |
| 1376 |
|
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays |
| 1377 |
< |
and other periodic structures. We have repeated the Luttinger \& |
| 1358 |
< |
Tisza series summations to much higher order and obtained the energy |
| 1359 |
< |
constants (converged to one part in $10^9$) in table \ref{tab:LT}. |
| 1377 |
> |
and other periodic structures. |
| 1378 |
|
|
| 1379 |
< |
\begin{table*}[h] |
| 1380 |
< |
\centering{ |
| 1381 |
< |
\caption{Luttinger \& Tisza arrays and their associated |
| 1382 |
< |
energy constants. Type ``A'' arrays have nearest neighbor strings of |
| 1383 |
< |
antiparallel dipoles. Type ``B'' arrays have nearest neighbor |
| 1384 |
< |
strings of antiparallel dipoles if the dipoles are contained in a |
| 1385 |
< |
plane perpendicular to the dipole direction that passes through |
| 1368 |
< |
the dipole.} |
| 1369 |
< |
} |
| 1370 |
< |
\label{tab:LT} |
| 1371 |
< |
\begin{ruledtabular} |
| 1372 |
< |
\begin{tabular}{cccc} |
| 1373 |
< |
Array Type & Lattice & Dipole Direction & Energy constants \\ \hline |
| 1374 |
< |
A & SC & 001 & -2.676788684 \\ |
| 1375 |
< |
A & BCC & 001 & 0 \\ |
| 1376 |
< |
A & BCC & 111 & -1.770078733 \\ |
| 1377 |
< |
A & FCC & 001 & 2.166932835 \\ |
| 1378 |
< |
A & FCC & 011 & -1.083466417 \\ |
| 1379 |
< |
B & SC & 001 & -2.676788684 \\ |
| 1380 |
< |
B & BCC & 001 & -1.338394342 \\ |
| 1381 |
< |
B & BCC & 111 & -1.770078733 \\ |
| 1382 |
< |
B & FCC & 001 & -1.083466417 \\ |
| 1383 |
< |
B & FCC & 011 & -1.807573634 \\ |
| 1384 |
< |
-- & BCC & minimum & -1.985920929 \\ |
| 1385 |
< |
\end{tabular} |
| 1386 |
< |
\end{ruledtabular} |
| 1387 |
< |
\end{table*} |
| 1388 |
< |
|
| 1389 |
< |
In addition to the A and B arrays, there is an additional minimum |
| 1379 |
> |
To test the new electrostatic methods, we have constructed very large, |
| 1380 |
> |
$N=$ 16,000~(bcc) arrays of dipoles in the orientations described in |
| 1381 |
> |
Ref. \onlinecite{LT}. These structures include ``A'' lattices with |
| 1382 |
> |
nearest neighbor chains of antiparallel dipoles, as well as ``B'' |
| 1383 |
> |
lattices with nearest neighbor strings of antiparallel dipoles if the |
| 1384 |
> |
dipoles are contained in a plane perpendicular to the dipole direction |
| 1385 |
> |
that passes through the dipole. We have also studied the minimum |
| 1386 |
|
energy structure for the BCC lattice that was found by Luttinger \& |
| 1387 |
|
Tisza. The total electrostatic energy for any of the arrays is given |
| 1388 |
|
by: |
| 1389 |
|
\begin{equation} |
| 1390 |
|
E = C N^2 \mu^2 |
| 1391 |
|
\end{equation} |
| 1392 |
< |
where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is |
| 1393 |
< |
the number of dipoles per unit volume, and $\mu$ is the strength of |
| 1394 |
< |
the dipole. |
| 1392 |
> |
where $C$ is the energy constant (equivalent to the Madelung |
| 1393 |
> |
constant), $N$ is the number of dipoles per unit volume, and $\mu$ is |
| 1394 |
> |
the strength of the dipole. Energy constants (converged to 1 part in |
| 1395 |
> |
$10^9$) are given in the supplemental information. |
| 1396 |
|
|
| 1397 |
< |
To test the new electrostatic methods, we have constructed very large, |
| 1398 |
< |
$N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in |
| 1399 |
< |
the orientations described in table \ref{tab:LT}. For the purposes of |
| 1400 |
< |
testing the energy expressions and the self-neutralization schemes, |
| 1401 |
< |
the primary quantity of interest is the analytic energy constant for |
| 1402 |
< |
the perfect arrays. Convergence to these constants are shown as a |
| 1403 |
< |
function of both the cutoff radius, $r_c$, and the damping parameter, |
| 1404 |
< |
$\alpha$ in Figs. \ref{fig:energyConstVsCutoff} and XXX. We have |
| 1405 |
< |
simultaneously tested a hard cutoff (where the kernel is simply |
| 1406 |
< |
truncated at the cutoff radius), as well as a shifted potential (SP) |
| 1407 |
< |
form which includes a potential-shifting and self-interaction term, |
| 1411 |
< |
but does not shift the forces and torques smoothly at the cutoff |
| 1412 |
< |
radius. The SP method is essentially an extension of the original |
| 1413 |
< |
Wolf method for multipoles. |
| 1397 |
> |
For the purposes of testing the energy expressions and the |
| 1398 |
> |
self-neutralization schemes, the primary quantity of interest is the |
| 1399 |
> |
analytic energy constant for the perfect arrays. Convergence to these |
| 1400 |
> |
constants are shown as a function of both the cutoff radius, $r_c$, |
| 1401 |
> |
and the damping parameter, $\alpha$ in Figs. |
| 1402 |
> |
\ref{fig:energyConstVsCutoff} and XXX. We have simultaneously tested a |
| 1403 |
> |
hard cutoff (where the kernel is simply truncated at the cutoff |
| 1404 |
> |
radius), as well as a shifted potential (SP) form which includes a |
| 1405 |
> |
potential-shifting and self-interaction term, but does not shift the |
| 1406 |
> |
forces and torques smoothly at the cutoff radius. The SP method is |
| 1407 |
> |
essentially an extension of the original Wolf method for multipoles. |
| 1408 |
|
|
| 1409 |
|
\begin{figure}[!htbp] |
| 1410 |
|
\includegraphics[width=4.5in]{energyConstVsCutoff} |
| 1429 |
|
approximation appears to perturb the potential too much inside the |
| 1430 |
|
cutoff region to provide accurate measures of the energy constants. |
| 1431 |
|
|
| 1438 |
– |
|
| 1432 |
|
{\it Quadrupolar} analogues to the Madelung constants were first |
| 1433 |
|
worked out by Nagai and Nakamura who computed the energies of selected |
| 1434 |
|
quadrupole arrays based on extensions to the Luttinger and Tisza |
| 1435 |
|
approach.\cite{Nagai01081960,Nagai01091963} We have compared the |
| 1436 |
|
energy constants for the lowest energy configurations for linear |
| 1437 |
< |
quadrupoles shown in table \ref{tab:NNQ} |
| 1437 |
> |
quadrupoles. |
| 1438 |
|
|
| 1446 |
– |
\begin{table*} |
| 1447 |
– |
\centering{ |
| 1448 |
– |
\caption{Nagai and Nakamura Quadurpolar arrays}} |
| 1449 |
– |
\label{tab:NNQ} |
| 1450 |
– |
\begin{ruledtabular} |
| 1451 |
– |
\begin{tabular}{ccc} |
| 1452 |
– |
Lattice & Quadrupole Direction & Energy constants \\ \hline |
| 1453 |
– |
SC & 111 & -8.3 \\ |
| 1454 |
– |
BCC & 011 & -21.7 \\ |
| 1455 |
– |
FCC & 111 & -80.5 |
| 1456 |
– |
\end{tabular} |
| 1457 |
– |
\end{ruledtabular} |
| 1458 |
– |
\end{table*} |
| 1459 |
– |
|
| 1439 |
|
In analogy to the dipolar arrays, the total electrostatic energy for |
| 1440 |
|
the quadrupolar arrays is: |
| 1441 |
|
\begin{equation} |
| 1442 |
|
E = C \frac{3}{4} N^2 Q^2 |
| 1443 |
|
\end{equation} |
| 1444 |
< |
where $Q$ is the quadrupole moment. |
| 1444 |
> |
where $Q$ is the quadrupole moment. The lowest energy |
| 1445 |
|
|
| 1446 |
|
\section{Conclusion} |
| 1447 |
|
We have presented two efficient real-space methods for computing the |
