#LyX 2.0 created this file. For more info see http://www.lyx.org/
\lyxformat 413
\begin_document
\begin_header
\textclass memoir
\use_default_options false
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman palatino
\font_sans lmss
\font_typewriter luximono
\font_default_family rmdefault
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100
\font_tt_scale 100

\graphics default
\default_output_format default
\output_sync 1
\bibtex_command default
\index_command default
\paperfontsize default
\spacing single
\use_hyperref false
\papersize default
\use_geometry false
\use_amsmath 1
\use_esint 1
\use_mhchem 1
\use_mathdots 1
\cite_engine natbib_numerical
\use_bibtopic true
\use_indices false
\paperorientation portrait
\suppress_date false
\use_refstyle 1
\index Index
\shortcut idx
\color #008000
\end_index
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\quotes_language english
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header

\begin_body

\begin_layout Chapter
Ewald summation
\begin_inset CommandInset label
LatexCommand label
name "chap:ewald-summation"

\end_inset


\end_layout

\begin_layout Standard

\emph on
Ewald summation
\emph default

\begin_inset CommandInset citation
LatexCommand citep
key "Ewald1923-Ewald-summation"

\end_inset

 is a method of rapidly evaluating fields at positions in an infinite lattice.
 It involves performing two sums, the first in real space over a lattice
 which has had an artificial screening added to reduce the importance of
 long-range interactions, and the second in reciprocal space.
 The screening in real space reduces the number of poles which need to be
 summed over and, even after the additional computational effort required
 to perform the reciprocal space sum is taken into account, this significantly
 improves the efficiency of the calculation.
\end_layout

\begin_layout Standard
First, Madelung constants are described which, though they are not directly
 relevant to the ultimate goal of magnetic dipole field simulation, provide
 a more comprehensible system in which to explain the process of Ewald summation.
 Then, the Madelung constant of sodium chloride is investigated and used
 to illustrate the problem of conditional convergence.
 Ewald summation is presented as a solution to this issue, and its generalisatio
n to other quantities, specifically dipole fields, is then briefly explored.
\end_layout

\begin_layout Section
Madelung constants
\end_layout

\begin_layout Standard
The 
\emph on
Madelung energy
\emph default

\begin_inset CommandInset citation
LatexCommand citep
key "Madelung1918"

\end_inset

 of a crystal is the electrostatic potential energy of an ion within it.
 This is dependent on the charges on and configuration of the ions which
 surround it.
 Evaluating it for an ion 
\begin_inset Formula $i$
\end_inset

 requires performing the sum
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
E_{i}=\sum_{j\neq i}\frac{q_{i}q_{j}}{4\pi\epsilon_{0}r_{j}}\,,
\end{equation}

\end_inset

where 
\begin_inset Formula $q_{i,j}$
\end_inset

 are the charges on the 
\begin_inset Formula $i$
\end_inset


\begin_inset script superscript

\begin_layout Plain Layout
th
\end_layout

\end_inset

 and 
\begin_inset Formula $j$
\end_inset


\begin_inset script superscript

\begin_layout Plain Layout
th
\end_layout

\end_inset

 ions, 
\begin_inset Formula $\epsilon_{0}$
\end_inset

 is the permittivity of free space, and the index 
\begin_inset Formula $j$
\end_inset

 implies summation over the infinite array of ions in the crystal.
 In many simple systems where the charges on ions have the same magnitude
 but may take either sign, the sum can be simplified to
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
E=\frac{q^{2}}{4\pi\epsilon_{0}}\sum_{j\neq i}\frac{(\pm)_{j}}{r_{j}},
\end{equation}

\end_inset

where 
\begin_inset Formula $(\pm)_{j}$
\end_inset

 is 
\begin_inset Note Note
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
num[explicit-sign]{+1}
\end_layout

\end_inset


\end_layout

\end_inset


\begin_inset Formula $+1$
\end_inset

 for a positively-charged ion 
\begin_inset Formula $j$
\end_inset

 and 
\begin_inset ERT
status collapsed

\begin_layout Plain Layout


\backslash
num[explicit-sign]{-1}
\end_layout

\end_inset

 if the ion carries negative charge.
 The 
\emph on
Madelung constant
\emph default
 is then defined as
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
M=\sum_{j\neq i}(\pm)_{j}\frac{r_{0}}{r_{j}},
\end{equation}

\end_inset

where 
\begin_inset Formula $r_{0}$
\end_inset

 is the nearest-neighbour distance in the lattice.
 It is thus independent of both the crystal dimensions and the units of
 charge, meaning that one Madelung constant applies to many systems but
 for a prefactor.
 Thus defined, the constant must be positive in order for a lattice configuratio
n to be energetically stable.
\end_layout

\begin_layout Standard
To take a very simple example, the Madelung constant of a 1D chain of alternatin
g charges can be expressed analytically:
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{eqnarray}
M & = & \sum_{j\neq i}\left(\pm\right){}_{j}\frac{r_{0}}{r_{j}}\nonumber \\
 & = & 2\left(1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+...\right)\nonumber \\
 & = & 2\ln2\nonumber \\
 & = & \num{1.386294}...,\label{eq:madelung-1d-plus-minus}
\end{eqnarray}

\end_inset

where the factor of 2 arises from the fact that the chain of ions extends
 to infinity in both directions from the central ion, and using the identity
 
\begin_inset Formula $\ln\left(1+x\right)=x-\frac{x^{2}}{2}+\frac{x^{3}}{3}-\frac{x^{4}}{4}...$
\end_inset

.
 
\end_layout

\begin_layout Section
NaCl and conditional convergence
\end_layout

\begin_layout Standard
Evaluating Madelung constants for 3D crystals is somewhat more complex.
 Let us take as an example sodium chloride, and consider evaluating its
 Madelung constant
\begin_inset CommandInset citation
LatexCommand citep
key "Grosso2001-cubic-Madelung"

\end_inset

.
 (Madelung constants for several common crystal structures are shown in
 
\begin_inset CommandInset ref
LatexCommand formatted
reference "tab:madelung-constants"

\end_inset

.)
\begin_inset Float table
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Tabular
<lyxtabular version="3" rows="4" columns="2">
<features booktabs="true" tabularvalignment="middle">
<column alignment="center" valignment="top" width="0pt">
<column alignment="center" valignment="top" width="0">
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $M$
\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
sodium chloride (NaCl)
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
num{1.747564595}
\end_layout

\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
cesium chloride (CsCl)
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
num{1.762674773}
\end_layout

\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" bottomline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
zincblende (ZnS)
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" bottomline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
num{1.638055062}
\end_layout

\end_inset


\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Selected Madelung constants, referred to the nearest-neighbour distance
\begin_inset CommandInset citation
LatexCommand citep
key "Grosso2001-cubic-Madelung"

\end_inset

.
\begin_inset CommandInset label
LatexCommand label
name "tab:madelung-constants"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset

The NaCl structure comprises two interpenetrating face-centred cubic lattices,
 one of Na
\begin_inset script superscript

\begin_layout Plain Layout
+
\end_layout

\end_inset

 and the other of Cl
\begin_inset script superscript

\begin_layout Plain Layout
\begin_inset Formula $-$
\end_inset


\end_layout

\end_inset

 ions.
 This can be imagined as a simple cubic lattice, with alternating ions on
 each lattice site, like a three-dimensional chessboard of charge (see 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:nacl-structure"

\end_inset

).
 
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename ../chapters/figures/dipole/ewald/nacl-structure.eps

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The structure of sodium chloride.
 
\begin_inset CommandInset label
LatexCommand label
name "fig:nacl-structure"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset

The Madelung constant must be extracted numerically.
 However, the algorithm chosen has a curious effect on whether or not it
 is possible to obtain the answer.
 An intuitive method for evaluating this sum would be to take spheres of
 increasing radius about the origin and add up the electrostatic contribution
 from each new shell of ions.
 Since ions at larger radii contribute less to the potential energy, it
 seems logical that the contribution from each successive shell should fall
 off and, with sufficiently large radius, an accurate answer might be obtained.
 However, this procedure does not result in a converging series---it has,
 in fact, been shown mathematically that this series 
\emph on
never
\emph default
 converges
\begin_inset CommandInset citation
LatexCommand citep
key "Borwein1985-Madelung-convergence"

\end_inset

.
 If you instead add up successive cubic shells centred on the origin, the
 series does converge on the correct answer
\begin_inset CommandInset citation
LatexCommand citep
key "Borwein1985-Madelung-convergence"

\end_inset

 (as shown in 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:madelung-conditional-convergence"

\end_inset

).
 This phenomenon, where convergence is contingent on the order in which
 you add the terms in a given series, is known as 
\emph on
conditional convergence
\emph default

\begin_inset CommandInset citation
LatexCommand citep
key "deLeeuw1980-lattice-sums"

\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename ../chapters/figures/dipole/madelung-conditional-convergence/madelung-conditional-convergence.eps

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Conditional convergence in the Madelung constant for sodium chloride.
 Summing over successive spherical shells fails to cause the series to converge,
 but the true value is arrived at relatively quickly if cubic sections are
 considered instead.
 The calculated value for the Madelung constant, 
\begin_inset Formula $M$
\end_inset

, is shown against the radius of the sphere, or half-edge of the cube, 
\begin_inset Formula $R$
\end_inset

 in terms of nearest-neighbour distances, equal to half the cubic lattice
 constant, 
\begin_inset Formula $a$
\end_inset

.
 The right 
\begin_inset Formula $y$
\end_inset

-axis label indicates the true value, 
\begin_inset Formula $M=\num{1.748}$
\end_inset

.
\begin_inset CommandInset label
LatexCommand label
name "fig:madelung-conditional-convergence"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
In the case of the sodium chloride lattice sum, the additional condition
 which needs to be satisfied in order to extract the correct lattice constant
 is charge neutrality, and charges in concentric spherical shells very rarely
 cancel out precisely.
 However, in cases where this is approached, so too is the correct value
 of the Madelung constant: the 
\begin_inset ERT
status collapsed

\begin_layout Plain Layout


\backslash
num{342648}
\end_layout

\end_inset


\begin_inset script superscript

\begin_layout Plain Layout
th
\end_layout

\end_inset

 term in the infinite sum over spherical shells has a net charge of 
\begin_inset ERT
status collapsed

\begin_layout Plain Layout


\backslash
num{2}
\end_layout

\end_inset

, and the calculated Madelung constant is 
\begin_inset ERT
status collapsed

\begin_layout Plain Layout


\backslash
num{1.748042464}
\end_layout

\end_inset

, relatively close to the correct value
\begin_inset CommandInset citation
LatexCommand citep
key "Grosso2001-cubic-Madelung"

\end_inset

.
\end_layout

\begin_layout Standard
As well as being a mathematical curiosity, conditional convergence is a
 potential hazard in performing any lattice sum: experimentation with a
 variety of summation methodologies would be required to have any assurance
 that, if convergence were achieved at all, the answer arrived at was correct.
 Worse, those sums which do converge often do so very slowly; with long-range
 interactions which fall off as 
\begin_inset Formula $1/r$
\end_inset

, many terms in the sum must be evaluated for every additional decimal place
 of accuracy.
 These issues can be overcome by introducing Ewald summation.
\end_layout

\begin_layout Section
Ewald summation
\end_layout

\begin_layout Standard
Ewald summation
\begin_inset CommandInset citation
LatexCommand citep
key "Ewald1923-Ewald-summation"

\end_inset

 splits a sum into two parts---a real space sum, which deals with the short-rang
e part of an interaction, and a reciprocal space sum, which accounts for
 long-range component.
 By sharing the load wisely between these two different sums, rapid convergence
 can be achieved and time spent computing lattice sums much reduced.
 Another advantage is that, by judicious choice of the 
\emph on
convergence factor
\emph default
 which makes the real-space interactions short-range, the sum can be made
 
\emph on
absolutely
\emph default
 and 
\emph on
uniformly
\emph default
 convergent for lattice sums involving spherical shells
\begin_inset CommandInset citation
LatexCommand citep
key "deLeeuw1980-lattice-sums"

\end_inset

, meaning that choice of summation algorithm is no longer an issue: naïvely
 using spheres alway works.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename ../chapters/figures/dipole/ewald/ewald.eps
	scale 90

\end_inset


\begin_inset Caption

\begin_layout Plain Layout
The components of an Ewald summation to evaluate the potential of the 
\begin_inset Formula $i$
\end_inset


\begin_inset script superscript

\begin_layout Plain Layout
th
\end_layout

\end_inset

 ion at the centre of an infinite 1D chain of alternating charges.
 All graphs show the electric charge 
\begin_inset Formula $Q$
\end_inset

 as a function of distance 
\begin_inset Formula $x$
\end_inset

 along the chain.
 (a)
\begin_inset space ~
\end_inset

The real-space potential 
\begin_inset Formula $V_{\mathrm{real}}$
\end_inset

 arises from the delta functions of charge in conjunction with the screening
 Gaussian distributions of opposite sign.
 The central ion 
\begin_inset Formula $i$
\end_inset

 is excluded from this sum.
 (b)
\begin_inset space ~
\end_inset

The reciprocal space potential, 
\begin_inset Formula $V_{\mathrm{reciprocal}}$
\end_inset

, is summed over an infinite lattice of Gaussians, with opposite sign to
 those in real space.
 Being periodic in real space, it is impossible to exclude the central charge
 
\begin_inset Formula $i$
\end_inset

.
 (c)
\begin_inset space ~
\end_inset

The self-term, 
\begin_inset Formula $V_{\mathrm{self}}$
\end_inset

, cancels the extraneous Gaussian at 
\begin_inset Formula $i$
\end_inset

 in 
\begin_inset Formula $V_{\mathrm{reciprocal}}$
\end_inset

.
 (d)
\begin_inset space ~
\end_inset

The resulting potential 
\begin_inset Formula $V$
\end_inset

 is thus calculated from a charge distribution containing only the original
 
\begin_inset Formula $\delta$
\end_inset

-functions.
\begin_inset CommandInset label
LatexCommand label
name "fig:ewald-summation"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-summation"

\end_inset

 illustrates the problem for the case of a 1D lattice comprising alternating
 positive and negative charges.
 We wish to know the electrostatic potential energy of the central ion in
 an infinite line of charges.
 Though this scenario involves a 1D lattice extending to infinity in the
 
\begin_inset Formula $\pm x$
\end_inset

 directions, the mathematics which follows is presented in terms of the
 3D vector 
\begin_inset Formula $\bm{r}$
\end_inset

 for generality.
 We start by defining a charge distribution 
\begin_inset Formula $\rho_{i}$
\end_inset

 at every lattice site 
\begin_inset Formula $i$
\end_inset

; simply a 
\begin_inset Formula $\delta$
\end_inset

-function of the appropriate charge 
\begin_inset Formula $q_{i}$
\end_inset

 located at the site's position 
\begin_inset Formula $\bm{r}_{i}$
\end_inset

,
\begin_inset Formula 
\begin{equation}
\rho_{i}(\bm{r})=q_{i}\delta(\bm{r}-\bm{r}_{i}).
\end{equation}

\end_inset

The first stage of Ewald summation is then shown in 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-summation"

\end_inset


\begin_inset space ~
\end_inset

(a): we reduce the range of the interactions in real space by adding screening
 charge distributions over every ion.
 A Gaussian screening distribution is usually chosen, though the choice
 of function is essentially arbitrary.
 The integral of this Gaussian charge density is identical to that of the
 
\begin_inset Formula $\delta$
\end_inset

-function point charge it screens, but is of opposite sign.
 This means that the net electric charge seen at large distances is zero,
 thus eliminating long-range interactions.
 The distribution used is
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
f_{i}(\bm{r})=q_{i}(\alpha^{2}/\pi)^{\frac{3}{2}}\mathrm{e}^{-\alpha^{2}(r-r_{i})^{2}},
\end{equation}

\end_inset

which means that the Gaussian is defined in terms of a screening parameter
 
\begin_inset Formula $\alpha$
\end_inset

.
 Higher 
\begin_inset Formula $\alpha$
\end_inset

 implies a narrower Gaussian, more sharply peaked around the 
\begin_inset Formula $\delta$
\end_inset

-function it screens, and therefore shrouds it more effectively.
\end_layout

\begin_layout Standard
The modified charge distribution in real space is therefore
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
\rho_{i}^{\prime}(\bm{r})=q_{i}\delta(\bm{r}-\bm{r}_{i})-q_{i}(\alpha^{2}/\pi)^{\frac{3}{2}}\mathrm{e}^{-\alpha^{2}(r-r_{i})^{2}}.
\end{equation}

\end_inset

Evaluating the potential at a point 
\begin_inset Formula $\bm{r}$
\end_inset

 due to this new distribution now involves three terms: the Coulomb term
 from the point charge at 
\begin_inset Formula $\bm{r}_{i}$
\end_inset

; the charge from the Gaussian inside the sphere centred on 
\begin_inset Formula $\bm{r}_{i}$
\end_inset

 with radius 
\begin_inset Formula $|\bm{r}-\bm{r}_{i}|$
\end_inset

, which is equivalent to that charge collapsed onto another 
\begin_inset Formula $\delta$
\end_inset

-function at 
\begin_inset Formula $\bm{r}_{i}$
\end_inset

; and the remaining Gaussian charge, outside this sphere.
 This must then be summed over all ions.
 Thus,
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
V_{\mathrm{real}}=\sum_{j\neq i}q_{j}\left(\frac{1}{r_{j}}-\frac{1}{r_{j}}\int_{0}^{r_{j}}\rho(\bm{r})\mathrm{d}\bm{r}-\int_{r_{j}}^{\infty}\frac{\rho(\bm{r})}{r}\mathrm{d}\bm{r}\right),
\end{equation}

\end_inset


\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename ../chapters/figures/dipole/ewald/erfc.eps

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The complementary error function, 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none

\begin_inset Formula $\mathrm{erfc}(x)=(2/\sqrt{\pi})\int_{x}^{\infty}\mathrm{e}^{-s^{2}}\mathrm{d}s$
\end_inset

.
\begin_inset CommandInset label
LatexCommand label
name "fig:erfc"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset

which simplifies pleasingly to
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
V_{\mathrm{real}}=\sum_{j\neq i}\frac{q_{j}}{r_{j}}\mathrm{erfc}(\alpha r_{j}),
\end{equation}

\end_inset

where 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
the complementary error function 
\begin_inset Formula $\mathrm{erfc}(x)=(2/\sqrt{\pi})\int_{x}^{\infty}\mathrm{e}^{-s^{2}}\mathrm{d}s$
\end_inset

 (see 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:erfc"

\end_inset

).
 The ion 
\begin_inset Formula $i$
\end_inset

 does not interact with itself, and thus the sum is defined over 
\begin_inset Formula $j\neq i$
\end_inset

.
 Since 
\begin_inset Formula $0<r_{j}<\infty$
\end_inset

, the complementary error function takes values 
\begin_inset Formula $\mathrm{1<erfc}(\alpha r_{j})<0$
\end_inset

, and premultiplying by it will reduce the Coulombic potential.
 It is worth emphasising that the Gaussian screening charges and their width
 as parameterised by 
\begin_inset Formula $\alpha$
\end_inset

 have no physical significance whatsoever---the screening parameter is chosen
 purely for convenience, with an eye to making sure that the sums converge,
 and do so as rapidly as possible.
\end_layout

\begin_layout Standard
The next stage is to remove the Gaussian distributions, leaving us with
 the original distribution of point charges we were attempting to sum over.
 This is shown in 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-summation"

\end_inset


\begin_inset space ~
\end_inset

(b).
 The reason that this is not an exercise in pointless mathematical circularity
 is that the lattice of Gaussians can be summed in reciprocal space.
 Evaluating this sum takes several stages but, by taking Fourier transforms
 of potential and charge, solving the Poisson equation and then manipulating
 the resulting expression, a sum in reciprocal space results
\begin_inset CommandInset citation
LatexCommand citep
key "Kittel-9th-edition"

\end_inset

:
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
V_{\mathrm{reciprocal}}=\frac{4\pi}{\Delta}\sum_{\bm{G\neq0}}\frac{S(\bm{G})}{G^{2}}\mathrm{e}^{\mathrm{i}\bm{G\cdot r}-G^{2}/4\alpha^{2}},
\end{equation}

\end_inset

where 
\begin_inset Formula $\Delta=\bm{a}\cdot\bm{b}\times\bm{c}$
\end_inset

 is the volume of the primitive unit cell in real space, 
\begin_inset Formula $\bm{G}$
\end_inset

 are the reciprocal lattice vectors and 
\begin_inset Formula $S(\bm{G})$
\end_inset

 is the structure factor in units of charge [
\begin_inset Formula $S(\bm{G})=\sum_{p}q_{p}\mathrm{e}^{\mathrm{i}\bm{G}\cdot\bm{r}_{p}}$
\end_inset

, where the index 
\begin_inset Formula $p$
\end_inset

 implies summation over the atoms in the primitive unit cell].
\end_layout

\begin_layout Standard
Components of a sum in reciprocal space are necessarily periodic in real
 space, and consequently it is not possible to exclude the Gaussian on the
 
\begin_inset Formula $i$
\end_inset


\begin_inset script superscript

\begin_layout Plain Layout
th
\end_layout

\end_inset

 atom, which is compensating for a Gaussian which we excluded in the 
\begin_inset Formula $V_{\mathrm{real}}$
\end_inset

 term.
 So, finally, the 
\emph on
self term
\emph default
 must be included, which removes this contribution:
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
V_{\mathrm{self}}=\int_{0}^{\infty}(4\pi r^{2}\mathrm{d}r)(\rho_{i}(r)/r)=\frac{2q_{i}\alpha}{\sqrt{\pi}}.
\end{equation}

\end_inset

This is shown in 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-summation"

\end_inset


\begin_inset space ~
\end_inset

(c).
\end_layout

\begin_layout Standard
The resulting formula for the potential on the ion 
\begin_inset Formula $i$
\end_inset

 is, then
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{eqnarray}
V_{i} & = & V_{\mathrm{real}}+V_{\mathrm{reciprocal}}+V_{\mathrm{self}}\nonumber \\
 & = & \sum_{j\neq i}\frac{q_{j}}{r_{j}}\mathrm{erfc}(\alpha r_{j})+\frac{4\pi}{\Delta}\sum_{\bm{G\neq0}}\frac{S(\bm{G})}{G^{2}}\mathrm{e}^{\mathrm{i}\bm{G\cdot r}-G^{2}/4\alpha^{2}}-\frac{2q_{i}\alpha}{\sqrt{\pi}}.\label{eq:ewald-sum}
\end{eqnarray}

\end_inset

This can be transformed into a formula for the Madelung constant by referring
 all distances 
\begin_inset Formula $r$
\end_inset

 to the nearest-neighbour distance 
\begin_inset Formula $r_{0}$
\end_inset

.
 Evaluating the potential at a lattice site is then a question of evaluating
 these two sums, and the self term.
 The final distribution of charge being summed over is thus, as shown in
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-summation"

\end_inset


\begin_inset space ~
\end_inset

(d), the same as the original.
\end_layout

\begin_layout Standard
As was mentioned above, choice of 
\begin_inset Formula $\alpha$
\end_inset

 is entirely pragmatic---its value is chosen to optimise the computational
 efficiency of the Ewald sum.
 Obtaining a result as quickly as possible is then a three-way optimisation
 between sphere size in real space, sphere size in reciprocal space and
 convergence parameter 
\begin_inset Formula $\alpha$
\end_inset

.
 A large value of 
\begin_inset Formula $\alpha$
\end_inset

 screens real-space charges more effectively, reducing the size of real-space
 sphere required, but has the opposite effect in reciprocal space.
 It can be shown that convergence occurs at the same rate in real and reciprocal
 space when 
\begin_inset Formula $\alpha=\sqrt{\pi}/a$
\end_inset

, where 
\begin_inset Formula $a$
\end_inset

 is the cubic lattice parameter
\begin_inset CommandInset citation
LatexCommand citep
key "Nijboer1957-lattice-sums"

\end_inset

.
 However, typically the reciprocal-space sum is more computationally efficient,
 and thus a larger value of 
\begin_inset Formula $\alpha$
\end_inset

 is often chosen to bias the summation to reciprocal space.
 (As a rule of thumb, 
\begin_inset Formula $\alpha$
\end_inset

 is usually sufficiently large to screen all but the charges in one unit
 cell in real space.) Various algorithms, such as one in Ref.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand citealp
key "Smith1998-Ewald"

\end_inset

, exist to optimise the three parameters, but their exact values are not
 crucial because Ewald summations tend to converge in a relatively large
 region of 
\begin_inset Formula $(r_{\mathrm{real}},\, r_{\mathrm{reciprocal}},\,\alpha)$
\end_inset

-space.
 An example optimisation in 
\begin_inset Formula $\alpha$
\end_inset

 is shown in 
\begin_inset CommandInset ref
LatexCommand formatted
reference "fig:ewald-alpha-converge"

\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename ../chapters/figures/dipole/ewald/nacl-madelung-ewald.eps

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Optimising an Ewald summation for the 3D NaCl lattice in the screening parameter
 
\begin_inset Formula $\alpha$
\end_inset

.
 Real space and reciprocal space spheres were fixed at 
\begin_inset Formula $r_{\mathrm{real}}=5a$
\end_inset

 and 
\begin_inset Formula $r_{\mathrm{reciprocal}}=5a^{*}$
\end_inset

, respectively.
 The sum converges by these sphere sizes for a relatively wide range of
 
\begin_inset Formula $\alpha$
\end_inset

 values.
 The right 
\begin_inset Formula $y$
\end_inset

-axis label indicates the true value, 
\begin_inset Formula $M=\num{1.748}$
\end_inset

.
\begin_inset CommandInset label
LatexCommand label
name "fig:ewald-alpha-converge"

\end_inset


\end_layout

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Section
Ewald summation for dipole-field simulations
\end_layout

\begin_layout Standard
In order to make use of Ewald summation for dipole-field simulations, a
 threefold generalisation of the Ewald method is required.
 Firstly, the muon will not sit at an atomic site, but at an interatomic
 position somewhere in the crystal.
 Generalising Ewald's formula to work out the potential at any position
 is simple
\begin_inset CommandInset citation
LatexCommand citep
key "Smith1998-Ewald"

\end_inset

:
\end_layout

\begin_layout Enumerate
We now sum over all ions, the 
\begin_inset Formula $j\neq i$
\end_inset

 condition no longer being relevant since we are not at an ion site 
\begin_inset Formula $i$
\end_inset

.
\end_layout

\begin_layout Enumerate
We remove the self-term.
\end_layout

\begin_layout Standard
We also need to generalise the sum to work for dipoles rather than monopoles
 and, further, fields rather than potentials.
 Ewald summation can be extended to calculate fields, field gradients and
 so on, and to dipoles or higher-order multipoles.
 Though the process is conceptually similar to monopole potential Ewald
 summation, it is rather more mathematically involved.
 Details can be found in e.g.
\begin_inset space ~
\end_inset

Refs.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand citealp
key "Nymand2000-Ewald-summation,Smith1998-Ewald,Kittel-9th-edition"

\end_inset

.
\end_layout

\end_body
\end_document
