International
Tables for
Crystallography
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 3.5, pp. 460-471

Section 3.5.2. Lattice sums of point charges

T. A. Dardena*

aLaboratory of Structural Biology, National Institute of Environmental Health Sciences, 111 T. W. Alexander Drive, Research Triangle Park, NC 27709, USA
Correspondence e-mail: darden@gamera.niehs.nih.gov

3.5.2. Lattice sums of point charges

| top | pdf |

We begin by discussing an idealized infinite crystal C, made up of point charges.

3.5.2.1. Basic quantities

| top | pdf |

The periodicity of an idealized crystal C is specified by lattice basis vectors [{\bf a}_1], [{\bf a}_2] and [{\bf a}_3], which are linearly independent as vectors and thus form a basis for the usual vector space of points in three dimensions. The conjugate reciprocal vectors [{\bf a}_{\alpha}^*] are defined by the relations [{\bf a}_{\alpha}^{*} \cdot {\bf a}_{\beta} = \delta_{\alpha \beta} ] (the Kronecker delta) for [\alpha,\beta = 1,2,3]. Thus for example [{\bf a}_{1}^{*}] is given by [{\bf a}_{1}^{*}] = [({\bf a}_2 \times {\bf a}_3) / [{\bf a}_1 \cdot ({\bf a}_2 \times {\bf a}_3)]] . General lattice vectors n are given by integral combinations of the lattice basis vectors:[{\bf n} = n_1{ \bf a}_1 + n_2{ \bf a}_2 + n_3 {\bf a}_3,\eqno(3.5.2.1) ]where [n_1,n_2,n_3] are integers. General reciprocal-lattice vectors m are given by integral combinations of the reciprocal-lattice basis vectors:[{\bf m} = m_1 {\bf a}_1^* + m_2 {\bf a}_2^* + m_3 {\bf a}_3^*.\eqno(3.5.2.2) ]

An arbitrary point r in the crystal (r also denotes the vector from the origin of coordinates to the point) is given by[{\bf r} = s_1 {\bf a}_1 + s_2 {\bf a}_2 + s_3 {\bf a}_3,]where [s_1,s_2,s_3 ] are the fractional coordinates of the point. From the above equation and the defining relations for the reciprocal basis we see that [s_{\alpha} = {\bf r}\cdot{\bf a}_{\alpha}^{*} ] for [{\alpha} = 1,2,3]. The unit cell U is the set of points r with associated fractional coordinates [s_1,s_2,s_3] satisfying [-\textstyle{1\over2} \le s_{\alpha} \le {1\over2}], [{\alpha} = 1,2,3]. The idealized infinite crystal C is generated by the union of all periodic translations of the unit cell U, using the set of general lattice vectors n. A point charge q at position r in U has periodic `image' charges [q_{\bf n} ] at positions [{\bf r}_{\bf n} = {\bf r} + {\bf n}] for all lattice vector n. The periodic image of U, denoted [{\it U}_{\bf n} ], is given by the periodic images [q_{\bf n}] of all point charges q in U. Thus, the crystal is made up of the union of images [U_{\bf n}], including [{\it U} = {\it U}_{\bf 0}], or[C = \bigcup_{\bf n} U_{\bf n}.]

A point charge [q_i] in U at position [{\bf r}_i] interacts with other unit-cell charges [q_{j}, j \not= i] at positions [{\bf r}_j] as well as with all of their periodic images [q_{j{\bf n}} ] at positions [{\bf r}_{j{\bf n}} = {\bf r}_{j}+{\bf n}] for all lattice vectors n. It also interacts with its own periodic images [q_{i{\bf n}}] at [{\bf r}_{i{\bf n}}] for all nonzero lattice vectors n, that is all [{\bf n} = n_1 {\bf a}_1 + n_2 {\bf a}_2 + n_3 {\bf a}_3] with [(n_{1},n_{2},n_{3}) \ne (0,0,0)]. Suppose the unit cell U is made up of N point charges [q_1,\ldots,q_N] at positions [{\bf r}_1,\ldots,{\bf r}_N]. Then the electrostatic energy of U, given by the sum of the Coulombic interactions of charges in U with each other and with the rest of C, is written as[\eqalignno{E(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= E({\bf r}_{1}, \ldots ,{\bf r}_{N}) &\cr&= {\textstyle{{1}\over{2}}} \sum_{\bf n} {}^{'} \sum_{i=1}^N \sum_{j=1}^N {{q_{i}q_{j}}\over{ | {\bf r}_{i} - {\bf r}_{j} - {\bf n}|}},&(3.5.2.3)} ]where here and below [\underline{\bf r}^{\lbrace\bf N\rbrace}]denotes the set of points [\lbrace{\bf r}_{1}, \ldots ,{\bf r}_{N}\rbrace ] and where the outer sum on the right is over the general lattice vectors n, the prime indicating that terms with i = j and n = 0 are omitted. The treatment of this infinite sum requires some care. The energy diverges unless the unit cell is neutral (i.e. unless [\textstyle\sum_{i=1}^N q_i = 0]). When U is neutral, the sum, if carried out in a naive fashion, can still converge quite slowly (see Chapter 3.4[link] ). In addition, the convergence is conditional, that is the resulting energy can vary depending on the order in which the sum is carried out! To better understand these phenomena, it is helpful to re-express the above energy. Let E(U, U) denote the Coulomb interaction energy of U with itself and [E({\it U},{\it U}_{\bf n}), {\bf n} \ne {\bf 0}] the Coulomb interaction energy of U with [{\it U}_{\bf n}], that is[E({\it U},{\it U}) = \sum_{i \,\lt\, j} {{q_i q_j}\over{|{\bf r}_i - {\bf r}_j|}} ]and[E({\it U},{\it U}_{\bf n}) = \sum_{i=1}^N\sum_{j=1}^N {{q_i q_j}\over{|{\bf r}_i - {\bf r}_j - {\bf n}|}} = \sum_{q_i \in {\it U}} \sum_{q_{j\,{\bf n}} \in {\it U}_{\bf n}} {{q_i q_{j\,{\bf n}}}\over{ | {\bf r}_{i} - {\bf r}_{j\,{\bf n}}|}}. ]Then we can abbreviate equation (3.5.2.3)[link] by[E(\underline{\bf r}^{\lbrace \bf N \rbrace}) = E({\it U},{\it U}) + \textstyle{{1}\over{2}} \textstyle\sum\limits_{{\bf n} \ne {\bf 0}} E({\it U},{\it U}_{\bf n}).\eqno(3.5.2.4) ]Note that for large n (that is, n with large norm [|{\bf n}|]) [E({\it U},{\it U}_{\bf n})] can be approximated by a Taylor expansion. Let [\bar{\bf r}] denote the centroid of points in U, given by [\bar{\bf r} = (1/N) \textstyle\sum_{i=1}^N {\bf r}_i ], and let [\bar{\bf r}_{\bf n}] denote the centroid of [{\it U}_{\bf n} ]. Then, expanding in terms of [{\bf r}_i - \bar{\bf r} = \Delta{\bf r}_i ] and [{\bf r}_{j\,{\bf n}} - \bar{\bf r}_{\bf n} = {\bf r}_{j} - \bar{\bf r} = \Delta{\bf r}_j ] (which are small compared to [\bar{\bf r} - \bar{\bf r}_{\bf n} = -{\bf n} ]) we have[\eqalign{{{1}\over{| {\bf r}_{i} - {\bf r}_{j\,{\bf n}}|}} & = {{1}\over{|(\bar{\bf r} - \bar{\bf r}_{\bf n}) + ({\bf r}_{i}- \bar{\bf r}) - ({\bf r}_{j\,{\bf n}} - \bar{\bf r}_{\bf n})|}} \cr &= {{1}\over{|{\bf n} - (\Delta{\bf r}_{i} - \Delta{\bf r}_{j})|}} \cr &\simeq {{1}\over{|{\bf n}|}} + \sum_{\alpha=1}^3{{{\bf n}_{\alpha}}\over{|{\bf n}|^3}} (\Delta{\bf r}_{i} - \Delta{\bf r}_{j})_{\alpha}\cr &\quad + {\textstyle{{1}\over{2}}}\sum_{\alpha,\beta = 1}^3 {{(3{\bf n}_{\alpha}{\bf n}_{\beta} - \delta_{\alpha \beta}{\bf n}^2)}\over{|{\bf n}|^5}} \cr &\quad\times (\Delta{\bf r}_{i}-\Delta{\bf r}_{j})_{\alpha} (\Delta{\bf r}_{i}-\Delta{\bf r}_{j})_{\beta} \cr &\quad+ \ldots.} ]

Let [Q = \textstyle\sum_{i=1}^N q_i] and [{\bf D} = \textstyle\sum_{i=1}^N q_i \Delta{\bf r}_i ] denote the net charge and dipole moment of the unit cell U, and let [{\bf G}_{\alpha\beta}] = [\sum_{i=1}^N q_i \Delta{\bf r}_{i\alpha} \Delta{\bf r}_{i\beta}] denote the quadrupole moment of the unit cell. Using the expansion in the above equation, we can Taylor expand the interaction energy [E({\it U},{\it U}_{\bf n}) ] as[\eqalignno{E({\it U},{\it U}_{\bf n}) &= \sum_{i=1}^N \sum_{j=1}^N q_i q_j {{1}\over{|{\bf n} - (\Delta{\bf r}_{i} - \Delta{\bf r}_{j})|}} &\cr& \simeq {{1}\over{|{\bf n}|}} Q^2 &\cr &\quad+ \sum_{\alpha,\beta = 1}^3 {{(3{\bf n}_{\alpha}{\bf n}_{\beta} - \delta_{\alpha \beta}{\bf n}^2)}\over{|{\bf n}|^5}} (Q {\bf G}_{\alpha\beta} - {\bf D}_\alpha {\bf D}_\beta) &\cr &\quad+ \ldots.&(3.5.2.5)} ]Note that the first-order term in the expansion has dropped out, since the dipole moments D of U and [{\it U}_{\bf n}] are the same. If [Q \ne 0] then for large [|{\bf n}|] the leading term in the Taylor expansion of [E({\it U},{\it U}_{\bf n})] is [Q^2 / |{\bf n}|]. Since [\textstyle\int\int\int ({{1}/{|{\bf r}|}} )\,{\rm d}^3{\bf r} = 4 \pi \textstyle\int_0^\infty r \,{\rm d}r = \infty ], then using a standard integral comparison [\textstyle\sum_{{\bf n} \ne {\bf 0}} 1 / |{\bf n}| ] is infinite, and so the resulting electrostatic energy of U diverges. If [Q = 0] we see that the leading term in the expansion of [E({\it U},{\it U}_{\bf n})] is given by [-\textstyle\sum_{\alpha,\beta}[({{3{\bf n}_{\alpha}{\bf n}_{\beta} - \delta_{\alpha \beta}{\bf n}^2})/{|{\bf n}|^5}}]{\bf D}_{\alpha}{\bf D}_{\beta}], that is a quadratic expression in the dipole moment D times a prefactor that decays like [1 / |{\bf n}|^3]. Since [\textstyle\int\int\int {{1}/{|{\bf r}|^3}} \,{\rm d}^3{\bf r} = 4 \pi \textstyle\int_0^\infty {{1}/{r}}\,{\rm d}r = \infty ], [\textstyle\sum_{{\bf n} \ne {\bf 0}} 1 / |{\bf n}|^3] also diverges (logarithmically), and so the electrostatic energy of U cannot converge absolutely in this case either. Hence, if it converges at all, it does so conditionally, that is the energy of U depends on the order in which we sum the contributions from image cells [{\it U}_{\bf n} ]. If Q = 0 and D = 0, the leading term in the Taylor expansion of [E({\it U},{\it U}_{\bf n})] depends on the quadrupole moment of the unit cell and decays like [1 / |{\bf n}|^5] for large [|{\bf n}|], and thus since [\textstyle\sum_{{\bf n} \ne {\bf 0}} 1 / |{\bf n}|^5 ] does converge, the electrostatic energy of U converges absolutely. From these considerations we might hope to express the energy of a neutral unit cell U in terms of an absolutely convergent sum plus a term quadratic in the unit-cell dipole moment D which also depends on the method of summation. As we will see below this is precisely what the Ewald summation method gives us.

3.5.2.2. Ewald sum: first derivation

| top | pdf |

In this section we give the standard Ewald-sum derivation found in textbooks such as Kittel (1986[link]). To do this we need two preliminary results, namely the electrostatic potential due to an isolated Gaussian charge distribution having total charge of one, and then that due to a periodic array of such distributions.

First we derive the electrostatic potential [\Psi({\bf r})] due to an isolated Gaussian density [\rho({\bf r})] given by[\rho({\bf r}) = (\alpha/\pi)^{3/2} \exp(-\alpha {\bf r}^2).\eqno(3.5.2.6) ]Owing to the spherical symmetry of [\rho({\bf r})] we can define the one-dimensional function [\psi(r)] by [\psi(r) = \Psi({\bf r}) ] for any r with [|{\bf r}| = r]. We next invoke Gauss' law (Eyges, 1980[link]) under spherical symmetry to evaluate the electric field [{\bf E}({\bf r})] due to [\rho] at any distance [r = |{\bf r}|]. In the case of spherical symmetry Gauss' law states that the electric field at a point r, due to a spherically symmetric charge distribution centred at the origin, is the same as if all the charge of the distribution from the origin out to [|{\bf r}|] were collected into a point charge at the origin. More specifically, if S is the sphere of radius [r = |{\bf r}|] with surface [\partial {\bf S} ], Gauss' law states[({{1}/{4\pi}})\textstyle \int\limits_{\partial {\bf S}} {\bf E}\,{\rm d}{\bf s} = \textstyle\int\limits_{\bf S} \rho({\bf u}) \,{\rm d}^3{\bf u}, ]and thus, defining the scalar quantity [E(r) = |{\bf E}({\bf r})| ], we get[\eqalign{E(r) &= {1\over r^2} \int\limits_{\bf S} \rho({\bf u}) \,{\rm d}^3{\bf u}\cr & = {1\over r^2} \left(\,{\alpha\over\pi}\,\right)^{3/2} \int\limits_{u = 0}^{r} 4 \pi u^2 \exp(-\alpha u^2) \,{\rm d}u \cr & = {{\rm erf}(\alpha^{1/2} r)\over r^2} - {2(\alpha/\pi)^{1/2}\exp(-\alpha r^2)\over r}} ]and thus, since (again by symmetry) [E(r) = -({\rm d}/{\rm d}r)\psi(r) ], we can verify that[\psi(r) = {{\rm erf}(\alpha^{1/2} r)\over r} + C,\eqno(3.5.2.7) ]where C is an arbitrary constant of integration. We obtain a potential with the desired asymptotic behaviour, i.e. [\psi(r) \rightarrow 0 ] as [r \rightarrow \infty], by setting C = 0. Note then that [\psi(r) \simeq 1/r] for large r as expected.

Next we derive the electrostatic potential [\Phi^{\rm P}({\bf r})] due to a periodic array of Gaussian charge distributions, each of the form just studied. That is, for points r in the unit cell U, define the periodic density [\rho^{\rm P}({\bf r})] by[\rho^{\rm P}({\bf r}) = \textstyle\sum\limits_{\bf n} \rho({\bf r} + {\bf n}), ]where the density [\rho] is defined in equation (3.5.2.6)[link] and the sum is over all general lattice vectors n, defined in equation (3.5.2.1)[link]. Since [\rho^{\rm P} ] is periodic, i.e. [\rho^{\rm P}({\bf r}+{\bf n}_0) = \rho({\bf r}) ] for any general lattice vector [{\bf n}_0], we can expand it in a Fourier series over the unit cell:[\rho^{\rm P}({\bf r}) = \textstyle\sum\limits_{\bf m} C_{\bf m} \exp(2 \pi i {\bf m} \cdot {\bf r}),\eqno(3.5.2.8) ]where the sum is over all reciprocal-lattice vectors m, defined in equation (3.5.2.2)[link], and the Fourier coefficient [C_{\bf m}] is given by[\eqalignno{C_{\bf m} &= {1\over V} \int\limits_U \rho^{\rm P}({\bf r}) \exp(-2 \pi i {\bf m} \cdot {\bf r})\,{\rm d}^3{\bf r}&\cr &= {1\over V}\sum\limits_{\bf n} \int\limits_U \rho({\bf r} + {\bf n}) \exp[-2 \pi i {\bf m} \cdot ({\bf r}+{\bf n})] \,{\rm d}^3{\bf r}&\cr &= {1\over V} \int\limits_{\Re^3} \rho({\bf r}) \exp(-2 \pi i {\bf m} \cdot {\bf r})\,{\rm d}^3{\bf r}&\cr &= {1\over V} \left(\,{\alpha\over\pi}\,\right)^{3/2}\int\limits_{\Re^3} \exp[-\alpha({\bf r} + \pi i {\bf m})]^2 \exp( -\pi^2 {\bf m}^2 / \alpha)\,{\rm d}^3{\bf r} &\cr & = {1\over V}\exp( -\pi^2 {\bf m}^2 / \alpha),&(3.5.2.9)} ]where V is the volume of the unit cell U. The periodicity of [\rho^{\rm P}] implies that of [\Phi^{\rm P}], so we can also expand the latter in a Fourier series:[\Phi^{\rm P}({\bf r}) = \textstyle\sum\limits_{\bf m} B_{\bf m} \exp(2 \pi i {\bf m} \cdot {\bf r}).\eqno(3.5.2.10) ]Next we recall Poisson's equation (Eyges, 1980[link]) for [\Phi^{\rm P}] in terms of [\rho^{\rm P}]:[\nabla \cdot \nabla \Phi^{\rm P}({\bf r}) = -4 \pi \rho^{\rm P}({\bf r}).\eqno(3.5.2.11) ]Using equations (3.5.2.10)[link] and (3.5.2.8)[link] we can rewrite equation (3.5.2.11)[link] as[\textstyle\sum\limits_{\bf m}(-4 \pi^2 {\bf m}^2) B_{\bf m} \exp(2 \pi i {\bf m} \cdot {\bf r}) = -4 \pi \textstyle\sum\limits_{\bf m} C_{\bf m} \exp(2 \pi i {\bf m} \cdot {\bf r}), ]thus using the uniqueness of Fourier-series coefficients [\pi {\bf m}^2 B_{\bf m} = C_{\bf m} ]. If [{\bf m} \ne {\bf 0}] this latter equation can be solved for [B_{\bf m}]. However, for m = 0 it requires that [C_{\bf m} = 0], whereas from equation (3.5.2.9)[link] [C_{\bf m} = 1 / V]. To avoid this contradiction we modify [\rho^{\rm P} ] by subtracting 1/V (uniform neutralizing plasma). That is, [\rho^{\rm P}] is now defined by[\rho^{\rm P}({\bf r}) = \textstyle\sum\limits_{\bf n} \rho({\bf r} + {\bf n}) - (1/V).\eqno(3.5.2.12) ]Note that [\rho^{\rm P}] is still periodic. Also note that for [{\bf m} \ne {\bf 0}] [C_{\bf m}] is still given by equation (3.5.2.9)[link], whereas now [C_{\bf 0} = 0]. Although we have avoided the contradiction, the value of [B_{\bf 0}] is still not determined. Some authors define it by considering m as a real variable and taking the limit as [{\bf m} \to {\bf 0}] (this depends on which direction the origin is approached from), but most treatments take it to be 0 (however, see Section 3.5.2.3[link]). In the latter case the electrostatic potential [\Phi^{\rm P}] due to [\rho^{\rm P}] is given by[\Phi^{\rm P}({\bf r}) = {1\over \pi V} \sum\limits_{{\bf m} \ne {\bf 0}} {\exp( -\pi^2 {\bf m}^2 / \alpha)\over {\bf m}^2 }\exp(2 \pi i {\bf m} \cdot {\bf r}).\eqno(3.5.2.13) ]

Now as before let [q_1,q_2,\ldots,q_N] be point charges located at positions [{\bf r}_1,{\bf r}_2,\ldots,{\bf r}_N] within the unit cell U. From the discussion in Section 3.5.2.1[link] we assume the unit cell is neutral, i.e. [q_1+\ldots+q_N = 0 ]. Let r be a point in U with [{\bf r} \ne {\bf r}_i,i = 1,\ldots,N ]. We want to obtain the electrostatic potential [\Phi] due to the point charges in U together with their periodic images in the remainder of the crystal C. That is,[\Phi_{\rm lattice}({\bf r}) = \sum\limits_{i=1}^N \sum\limits_{\bf n} {q_i\over |{\bf r}_i + {\bf n} - {\bf r}|}. ]We next supplement the potential due to the point charge [q_i] at [{\bf r}_i+{\bf n}] by the potential due to an equally charged Gaussian counterion, that is [(-q_i)] times a density of the form [\rho] in equation (3.5.2.6)[link], centred at [{\bf r}_i+{\bf n} ]. Using equation (3.5.2.7)[link] the resulting direct sum potential is given by[\Phi_{\rm dir}({\bf r}) = \sum\limits_{i=1}^N \sum\limits_{\bf n}{ q_i\, {\rm erfc}(\alpha^{1/2} |{\bf r}_i + {\bf n} - {\bf r}|)\over |{\bf r}_i + {\bf n} - {\bf r}|}, ]where [{\rm erfc}(x) = 1 - {\rm erf}(x)]. Note that since [{\rm erfc}(x)] decays exponentially fast to zero as [x \rightarrow \infty ], the sum in the above equation converges rapidly and absolutely.

However, having added the potential due to counterions we now need to cancel it. That is, we must consider the reciprocal electrostatic potential due to [(+q_i)] times a density of the form [\rho], centred at [{\bf r}_i + {\bf n}], for all direct-lattice vectors n. Since the unit cell is neutral, the sum of these reciprocal potentials over [i = 1,\ldots,N] and all n is the same as the sum over [i = 1,\ldots,N] of [q_i] times the potential due to [\rho^{\rm P} ] centred at [{\bf r}_i] where [\rho^{\rm P}] is given by equation (3.5.2.12)[link]. Note that the extra term involving 1/V drops out by neutrality. [In the case of non-neutrality, such as discussion of the free energy of charging an ion (Hummer et al., 1996[link]), the potential due to the uniform compensating plasma must be included (see below). Hummer has derived the potential due to a uniformly charged unit cell in the special case of a cubic unit cell (Hummer, 1996[link]).] Thus, using equation (3.5.2.13)[link] we arrive at the reciprocal sum potential:[\eqalign{\Phi_{\rm rec}({\bf r}) &= \sum\limits_{i=1}^N q_i \Phi^{\rm P}({\bf r} - {\bf r}_i) \cr& = {1\over\pi V} \sum_{{\bf m} \ne {\bf 0}} {\exp( -\pi^2 {\bf m}^2 / \alpha)\over {\bf m}^2 }\cr &\quad \times \sum\limits_{i=1}^N q_i \exp[2 \pi i {\bf m} \cdot ({\bf r} - {\bf r}_i)] \cr &= {1\over\pi V} \sum_{{\bf m} \ne {\bf 0}} {\exp( -\pi^2 {\bf m}^2 / \alpha)\over{\bf m}^2 }\exp(2 \pi i {\bf m} \cdot {\bf r}) S(-{\bf m}),} ]where the structure factor [S({\bf m})] is given by[S({\bf m}) = \textstyle\sum\limits_i^N q_i \exp(2 \pi i {\bf m} \cdot {\bf r}_i).\eqno(3.5.2.14) ]

We have now expressed [\Phi_{\rm lattice}({\bf r})] as the sum of [\Phi_{\rm dir}({\bf r})] plus [\Phi_{\rm rec}({\bf r})] (with the caveat that in the process some infinite series have been casually rearranged and that the contribution of the reciprocal vector m = 0 has been ignored). To get a valid electrostatic potential [\Phi_{\rm lattice}({\bf r}) ] at the position of a charge [q_i] in U we must remove the Coulomb potential due to [q_i] and then take the limit as [{\bf r} \rightarrow {\bf r}_i]. Noting that [1 = {\rm erf}(x) + {\rm erfc}(x) ] and that [\lim_{x \to 0} {\rm erf}(x) / x = 2 / \sqrt{\pi}], we have[\eqalign{\Phi_{\rm lattice}({\bf r}_i) &= \lim_{{\bf r} \to {\bf r}_i} \left[\Phi_{\rm lattice}({\bf r}) - {{q_i}\over{|{\bf r} - {\bf r}_i|}} \right]\cr &= \sum\limits_{j \ne i} {{q_j \ {\rm erfc}(\alpha^{1/2}|{\bf r}_j - {\bf r}_i|)}\over{|{\bf r}_j - {\bf r}_i|}} \cr &\quad+ \sum\limits_{j = 1}^N \sum_{{\bf n}\ne {\bf 0}} {{q_j \ {\rm erfc}(\alpha^{1/2}|{\bf r}_j + {\bf n} - {\bf r}_i|)}\over{|{\bf r}_j + {\bf n} - {\bf r}_i|}} \cr &\quad+ {{1}\over{\pi V}} \sum\limits_{{\bf m} \ne {\bf 0}} {{\exp( -\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}} \exp(2 \pi i {\bf m} \cdot {\bf r}_i) S(-{\bf m}) \cr &\quad- 2 q_i \left(\,{{\alpha}\over{\pi}} \,\right)^{1/2}.} ]Thus we can write the electrostatic energy of U from equation (3.5.2.3)[link] as[\eqalignno{E_{\rm Ewald}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum\limits_{i=1}^N q_i \, \Phi_{\rm lattice}({\bf r}_i) &\cr &= {\textstyle{{1}\over{2}} }\sum\limits_{\bf n} {}^{'} \sum\limits_{i=1}^N \sum\limits_{j=1}^N {{ q_i q_j \ {\rm erfc}(\alpha^{1/2}|{\bf r}_j + {\bf n} - {\bf r}_i|)}\over{|{\bf r}_j + {\bf n} - {\bf r}_i|}} &\cr &\quad+ {{1}\over{2 \pi V}} \sum\limits_{{\bf m} \ne {\bf 0}} {{\exp( -\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}} S({\bf m}) S(-{\bf m}) &\cr &\quad - \left(\,{{\alpha}\over{\pi}} \,\right)^{1/2} \sum\limits_{i=1}^N q_i^2,&(3.5.2.15)} ]where as above the prime indicates that terms with i = j and n = 0 are omitted.

3.5.2.3. Ewald sum: more complete derivation

| top | pdf |

The derivation in Section 3.5.2.2[link] is incomplete. In addition to the caveats concerning the expression for [\Phi_{\rm lattice}({\bf r}) ] in terms of [\Phi_{\rm dir}({\bf r})] plus [\Phi_{\rm rec}({\bf r}) ], there is no mention of a quadratic term involving the unit-cell dipole moment M. In this section we provide a more complete derivation of the Ewald sum, which will in addition explain the conditional convergence discussed previously. We follow the developments in Smith (1981[link], 1994[link]).

The derivation depends on the following identity:[\eqalignno{{{1}\over{|{\bf r}|}} &= {{{\rm erfc}(\alpha^{1/2} |{\bf r}|)}\over{|{\bf r}|}} + {{1}\over{\pi}} \int\limits_{\Re^3} {{\exp(-\pi^2 {\bf u}^2 / \alpha)}\over{{\bf u}^2}} \exp(-2 \pi i {\bf u}\cdot {\bf r}) \,{\rm d}^3 {\bf u} &\cr &= {{{\rm erfc}(\alpha^{1/2} |{\bf r}|)}\over{|{\bf r}|}} + {{1}\over{\pi}} \sum\limits_{\bf m} \int\limits_{U^{*}} {{\exp(-\pi^2|{\bf v}+{\bf m}|^2 / \alpha)}\over{|{\bf v}+{\bf m}|^2}} &\cr&\quad \times \exp[-2 \pi i ({\bf v}+{\bf m})\cdot {\bf r}]\,{\rm d}^3{\bf v},&(3.5.2.16)} ]where [U^{*}] is the reciprocal unit cell, i.e. the set of points [{\bf v} = t_1 a_1^{*}+t_2 a_2^{*} + t_3 a_3^{*}] with [-\textstyle{1\over2} \leq t_{\alpha} \leq \textstyle{1\over2}], [\alpha=1,2,3 ]. Actually we derive (following Essmann et al., 1995[link]) a generalization of identity (3.5.2.16)[link], valid for general inverse powers [1/|{\bf r}|^p], p > 0. For this first note that for u > 0[\eqalignno{\Gamma(u) &= \textstyle\int\limits_0^\infty t^{\,u-1} \exp(-t)\,{\rm d}t&\cr &= \lambda^u \textstyle\int\limits_0^\infty t^{\,u-1} \exp(-\lambda t)\,{\rm d}t,&(3.5.2.17)} ]where [\Gamma(u)] is the Euler gamma function. Next note that for a > 0[\eqalignno{ \exp(-a {\bf r}^2) &= ( {{\pi}/{a}})^{3/2}\textstyle \int\limits_{\Re^3} \exp(-\pi^2 {\bf u}^2 / a)&\cr &\quad\times \exp(-2 \pi i {\bf u}\cdot {\bf r})\,{\rm d}^3{\bf u}, &(3.5.2.18)} ]which is the Fourier transform of the Gaussian, as derived above within equation (3.5.2.9)[link]. In equation (3.5.2.17)[link] substitute [\lambda = |{\bf r}|^2] and u = p/2 to get, for α > 0,[\eqalign{{{\Gamma(p/2)}\over{|{\bf r}|^p}} &= \int\limits_0^\alpha t^{\,p/2 - 1} \exp(-{\bf r}^2 t) \,{\rm d}t + \int\limits_\alpha^\infty t^{\,p/2 - 1} \exp(-{\bf r}^2 t) \,{\rm d}t \cr&= {\rm I}_p + {\rm II}_p.} ]In [{\rm II}_p] we substitute [{\bf r}^2 t = s^2], whereas in [{\rm I}_p] we use equation (3.5.2.18)[link], change the order of integration, and substitute [\pi^2 {\bf u}^2 / t = s^2] to obtain[\eqalignno{{{1}\over{|{\bf r}|^p}}&= {{1}\over{\Gamma(p/2)}} {{1}\over{|{\bf r}|^p}} \int\limits_{\alpha^{1/2} |{\bf r}|}^\infty 2 s^{\,p- 1} \exp(-s^2)\,{\rm d}s&\cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{\Gamma(p/2)}} \int\limits_{\Re^3} {\rm d}^3{\bf u} \exp(-2 \pi i {\bf u}\cdot {\bf r}) \left({{\pi |{\bf u}| }\over {\alpha^{1/2}}}\right)^{p-3} &\cr &\quad\times \int\limits_{\pi|{\bf u}| / \alpha^{1/2}}^\infty 2 s^{\,2-p} \exp(-s^2) \,{\rm d}s&\cr &= {{f_p(\alpha^{1/2}|{\bf r}|)}\over{|{\bf r}|^p}}&\cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{\Gamma(p/2)}} \int\limits_{\Re^3} \exp(-2 \pi i {\bf u}\cdot {\bf r}) g_p(\pi|{\bf u}|/\alpha^{1/2})\,{\rm d}^3{\bf u} &\cr &= {{f_p(\alpha^{1/2}|{\bf r}|)}\over{|{\bf r}|^p}} &\cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{\Gamma(p/2)}} \sum\limits_{\bf m} \int\limits_{U^{*}} g_p\left({{\pi|{\bf v}+{\bf m}|}\over{\alpha^{1/2}}}\right) &\cr&\quad\times \exp[-2 \pi i ({\bf v}+{\bf m})\cdot {\bf r}]\,{\rm d}^3{\bf v},&(3.5.2.19)} ]where[f_p(x) = {{1}\over{\Gamma(p/2)}} \int\limits_{x}^\infty 2 s^{\,p- 1} \exp(-s^2)\,{\rm d}s \eqno(3.5.2.20) ]and[g_p(x) = x^{\,p-3} \int\limits_x^\infty 2 s^{\,2-p} \exp(-s^2)\,{\rm d}s.\eqno(3.5.2.21) ]Specializing to the case p = 1 we see that [f_p(x) = {\rm erfc}(x) ] and that [g_p(x) = \exp(-x^2) / x^2], which proves identity (3.5.2.16)[link]. Another important case is p = 6, for which [f_p(x)] = [(1 + x^2 + x^4/2) \exp(-x^2)] and [g_p(x)] = [(1 - 2x^2)\exp(-x^2) / 3 + 2 x^3 \pi^{1/2} {\rm erfc}(x) ].

Note that from equation (3.5.2.17)[link], substituting t = s2, λ = r2 and u = p/2 we have[{{1}\over{|{\bf r}|^p}} = {{2}\over{\Gamma(p/2)}} \int\limits_0^{\infty} s^{\,p-1} \exp(-r^2s^2) \,{\rm d}s. ]Also, substituting [s = |{\bf r}|t] into equation (3.5.2.20)[link] we have[{{f_p(\alpha^{1/2} |{\bf r}|)}\over{|{\bf r}|^p}} = {{2}\over{\Gamma(p/2)}} \int\limits_{\alpha^{1/2}}^\infty t^{\,p-1} \exp(-r^2t^2) \,{\rm d}t ]and so[\eqalignno{ &\lim_{{\bf r} \to {\bf 0}} {{1}\over{|{\bf r}|^p}} [1 - f_p(\alpha^{1/2} |{\bf r}|) ] &\cr &\quad= \lim_{{\bf r} \to {\bf 0}} {{2}\over{\Gamma(p/2)}} \int\limits_0^{\alpha^{1/2}} t^{\,p-1} \exp(-r^2t^2)\, {\rm d}t&\cr&\quad= {{2 \alpha^{\,p/2}}\over {p \Gamma(p/2)}}.&(3.5.2.22)} ]

Note that [g_p(x)] is a smooth function of x for [x \ne 0 ], and that its derivatives are uniformly bounded for x bounded away from 0. To show this, substitute s = xt in equation (3.5.2.21)[link] to write[g_p(x) = 2 \textstyle\int\limits_1^\infty t^{2-p} \exp(-x^2 t^2) \,{\rm d}t. ]

In the case that p = 6, [g_p(x)] is smooth with uniformly bounded derivatives for all x, but this is unfortunately not true for all p. However, from the above expression we can show that for p > 3, [g_p(x)] is bounded and continuous as [x \rightarrow 0 ], whereas for [1 \leq p \leq 3] [g_p(x) \rightarrow \infty ] as [x \rightarrow 0]. Thus for p > 3 we see that for any x > 0, [g_p(x) \,\lt\, g_p(0) = 2 / (p-3)]. However for [1 \leq p \leq 3], for any R > 1 if x > 0 sufficiently small [\exp(-x^2R^2) > \textstyle{1\over2}]; thus [g_p(x) >] [2 \textstyle\int_1^R t^{\,2-p}\exp(-x^2 t^2)\,{\rm d}t > \textstyle\int_1^R t^{\,2-p}\,{\rm d}t \geq \ln(R) ] and thus [g_p(x)] is unbounded as [x \rightarrow 0].

Despite this latter result, the integral of [g_p({\bf r})] over [\Re^3] is bounded for [p \ge 1]. This result is important below for the validity of rearranging the order of summation or taking limits inside infinite sums involving [g_p]. To prove boundedness, use the above expression for [g_p] and integrate by parts to write[\eqalign{\textstyle\int\limits_{\Re^3} g_p({\bf r}) \,{\rm d}^3{\bf r} &= 4 \pi \textstyle\int\limits_0^{\infty} x^2 g_p(x) \,{\rm d}x \cr &= 4 \pi \textstyle\int\limits_0^{\infty} \,{\rm d}x \textstyle\int\limits_1^\infty t^{1-p} (2 x^2 t)\exp(-x^2t^2) \,{\rm d}t \cr &= 4\pi \left[\textstyle\int\limits_0^\infty \exp(-x^2) \,{\rm d}x \right. \cr&\quad+ (1-p)\left. \textstyle\int\limits_1^\infty t^{-p} \,{\rm d}t \textstyle\int\limits_0^\infty \exp(-x^2t^2) \,{\rm d}x\right] \cr &= 4\pi \left[\textstyle\int\limits_0^\infty \exp(-x^2) \,{\rm d}x \right.\cr &\quad+ \left. (1-p) \textstyle\int\limits_1^{\infty} t^{-(p+1)} \,{\rm d}t \textstyle\int\limits_0^\infty \exp(-s^2)\,{\rm d}s\right] \cr &= {{2\pi^{3/2}}\over{p}}.} ]

3.5.2.3.1. The case of p > 3

| top | pdf |

We first discuss the lattice sum for the case p > 3. In this case the lattice sum converges absolutely, since the tail of this sum can be bounded from above by a constant times the integral[\int\limits_{|{\bf r}| > K} {{1}\over{|{\bf r}|^p}} \,{\rm d}^3{\bf r} = 4\pi \int\limits_{K}^\infty r^{(2-p)} \,{\rm d}r = {{K^{(3-p)}}\over{p-3}}\eqno(3.5.2.23) ]for some K > 0. Although the sum can thus be evaluated directly, it can be significantly accelerated using the same Ewald direct- and reciprocal-space decomposition. For [{\bf r} \ne {\bf 0}], using equation (3.5.2.19)[link] we have[\eqalign{\sum\limits_{\bf n} {{1}\over{|{\bf r}+{\bf n}|^p}} &= \sum\limits_{\bf n} {{f_p(\alpha^{1/2}|{\bf r}+{\bf n}|)}\over{|{\bf r} + {\bf n}|^p}} + {{\pi^{3/2}\alpha^{(p-3)/2}}\over{\Gamma(p/2)}} \cr &\quad \times \sum\limits_{\bf n}\sum_{\bf m} \int\limits_{U^{*}} h_{p,{\bf m},{\bf r}}({\bf v})\exp(-2\pi i {\bf v}\cdot{\bf n})\,{\rm d}^3{\bf v},} ]where [h_{p,{\bf m},{\bf r}}({\bf v}) = g_p(\pi|{\bf v}+{\bf m}|/\alpha^{1/2})\exp[-2 \pi i ({\bf v}+{\bf m})\cdot {\bf r}] ]. Since [f_p(\alpha^{1/2}|{\bf r}+{\bf n}|)] is bounded, the first of the above sums converges absolutely and we turn to the second. We recognize the integral [\textstyle\int_{U^{*}} h_{p,{\bf m},{\bf r}}({\bf v})\exp(-2\pi i {\bf v}\cdot{\bf n})\,{\rm d}^3{\bf v} ] as (1/V) times [\hat{h}_{p,{\bf m},{\bf r}}({\bf n})], the nth Fourier coefficient of [h_{p,{\bf m},{\bf r}}], where V is the volume of the unit cell U (1/V is the volume of [U^{*}]). From the above results for [g_p] it is clear that for [{\bf m} \ne {\bf 0}], [h_{p,{\bf m},{\bf r}}] are smooth and have uniformly (in v, r and [{\bf m} \ne {\bf 0}]) bounded derivatives. Thus the sum of the Fourier coefficients can be shown to converge absolutely and uniformly in m, which guarantees uniform (over [{\bf m} \ne {\bf 0}]) convergence of the sum of these integrals to [(1 / V)h_{p,{\bf m},{\bf r}}({\bf 0})], which in turn equals [(1/V) g_p(\pi|{\bf m}|/\alpha^{1/2})\exp(-2 \pi i {\bf m}\cdot {\bf r})].

The case of m = 0 is not as straightforward, since the uniform boundedness of derivatives of [h_{p,{\bf m},{\bf r}}] does not generally hold in this case. One way of handling it is as follows. Since everything else in the above equation converges as [|{\bf n}| \to \infty], the sum of the Fourier coefficients for [h_{p,{\bf m}={\bf 0},{\bf r}}] converges to something. That is, if[S(n_1,n_2,n_3) = \textstyle\sum\limits_{k_1 = -n_1}^{n_1} \textstyle\sum\limits_{k_2 = -n_2}^{n_2} \textstyle\sum\limits_{k_3 = -n_3}^{n_3} \hat{h}_{p,{\bf 0},{\bf r}}(k_1,k_2,k_3) ]for [n_\alpha > 0], [\alpha = 1,2,3], then for some s, [S(n_1,n_2,n_3) \to s] as [\max(n_1,n_2,n_3) \to \infty ]. If[\sigma(n) = [{{1}/({n+1})} ]^3 \textstyle\sum\limits_{n_1 = 0}^n \textstyle\sum\limits_{n_2 = 0}^n \textstyle\sum\limits_{n_3 = 0}^n S(n_1,n_2,n_3), ]then it is straightforward to show that [\sigma(n)] also converges to s as [n \to \infty] [see, for example, Chapter 1 of Körner (1988[link]) for the one-dimensional proof]. On the other hand, due to the boundedness and continuity of [h_p], [\sigma(n) \to h_{p,{\bf m},{\bf r}}({\bf 0}) ] as [n \to \infty]. [This involves the three-dimensional version of the Fejér kernel, which converges to the three-dimensional Dirac delta function. See, for example, Chapter 2 of Körner (1988[link]) for the one-dimensional proof.] Thus [S(n_1,n_2,n_3) \to h_{p,{\bf m},{\bf r}}({\bf 0}) ] as [\max(n_1,n_2,n_3) \to \infty]. Hence the above sum of integrals with m = 0 converges to [(1 / V)h_{p,{\bf 0},{\bf r}}({\bf 0}) ].

Thus from the above arguments we have[\eqalignno{ \sum\limits_{\bf n} {{1}\over{|{\bf r}+{\bf n}|^p}} &= \sum\limits_{\bf n} {{f_p(\alpha^{1/2}|{\bf r}+{\bf n}|)}\over{|{\bf r} + {\bf n}|^p}} &\cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{ \Gamma(p/2)}} {{1}\over{V}} \sum\limits_{\bf m} g_p(\pi|{\bf m}|/\alpha^{1/2}) &\cr &\quad\times \exp(-2 \pi i {\bf m}\cdot {\bf r}).&(3.5.2.24)} ]

Using the limiting behaviour given in (3.5.2.22)[link] we can extend this result to the case r = 0. By removing the singularity when n = 0 we get[\eqalignno{\sum\limits_{{\bf n} \ne {\bf 0}} {{1}\over{|{\bf n}|^p}} &= \lim_{{\bf r} \to {\bf 0}} \sum\limits_{{\bf n} \ne {\bf 0}} {{1}\over{|{\bf r}+{\bf n}|^p}} &\cr &= \lim_{{\bf r} \to {\bf 0}} \left( \sum\limits_{\bf n} {{1}\over{|{\bf r}+{\bf n}|^p}} - {{1}\over{|{\bf r}|^p}} \right) &\cr&= \sum\limits_{{\bf n} \ne {\bf 0}} {{f_p(\alpha^{1/2}|{\bf n}|)}\over{|{\bf n}|^p}} - {{2 \alpha^{\,p/2}}\over {p \Gamma(p/2)}} + {{\pi^{3/2}\alpha^{(p-3)/2}}\over{ \Gamma(p/2)}} &\cr&\quad\times {{1}\over{V}} \sum\limits_{\bf m} g_p(\pi|{\bf m}|/\alpha^{1/2}).&(3.5.2.25)} ]

Now assume particles i and j in the unit cell U interact by the potential [\phi_{ij} = C(i,j) / |{\bf r}_j - {\bf r}_i|^p ]. Using equations (3.5.2.24)[link] and (3.5.2.25)[link] (the latter for particles interacting with their own images) the total energy for the unit cell U consisting of N particles interacting via the potential [\varphi] with each other and all of the image cells is given by[\eqalign{E(\underline{\bf r}^{\lbrace \bf N\rbrace}) &= {\textstyle{{1}\over{2}} }\sum_{\bf n} {}^{'} \sum_{i=1}^N \sum_{j=1}^N {{C(i,j)}\over{ | {\bf r}_{i} - {\bf r}_{j} - {\bf n}|^p}} \cr &= {\textstyle{{1}\over{2}}} \sum_{\bf n} {}^{'} \sum_{i=1}^N \sum_{j=1}^N {{C(i,j) f_p(\alpha^{1/2}|{\bf r}_j - {\bf r}_i+{\bf n}|)}\over{|{\bf r}_j - {\bf r}_i + {\bf n}|^p}} \cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{ \Gamma(p/2)}} {{1}\over{2 V}} \sum_{\bf m} g_p(\pi|{\bf m}|/\alpha^{1/2}) \cr &\quad\times \sum_{i=1}^N \sum_{j=1}^N C(i,j) \exp[-2 \pi i {\bf m}\cdot ({\bf r}_j - {\bf r}_i)] \cr &\quad- {{\alpha^{\,p/2}}\over {p \Gamma(p/2)}} \sum_{i=1}^N C(i,i),} ]where as before the prime on the first summation means that terms where i = j and n = 0 are omitted.

In the case when [C(i,j)] factorizes, i.e. [C(i,j) = C(i)C(\,j) ], the above lattice sum can be further simplified and accelerated using the structure factor [S({\bf m})] given by[S({\bf m}) = \textstyle\sum\limits_{j=1}^N C(\,j) \exp(2 \pi i {\bf m}\cdot {\bf r}_i).\eqno(3.5.2.26) ]In this case,[\eqalign{E(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}} }\sum_{\bf n} {}^{'} \sum_{i=1}^N \sum_{j=1}^N {{C(i)C(\,j) f_p(\alpha^{1/2}|{\bf r}_j - {\bf r}_i+{\bf n}|)}\over{|{\bf r}_j - {\bf r}_i + {\bf n}|^p}} \cr &\quad+ {{\pi^{3/2}\alpha^{(p-3)/2}}\over{ \Gamma(p/2)}} {{1}\over{2 V}} \sum_{\bf m} g_p(\pi|{\bf m}|/\alpha^{1/2}) S({\bf m})S(-{\bf m})\cr & \quad- {{\alpha^{\,p/2}}\over {p \Gamma(p/2)}} \sum_{i=1}^N C(i)^2 .} ]

Note that as expected, due to the absolute convergence of the lattice sums, there is no shape-dependent correction term to the Ewald sum. This is connected with the fact that [g_p(x)] is bounded and continuous as [x \to 0]. Next we see what additional assumptions are necessary when this is no longer true, and what form the correction term takes. We cover the most interesting case, namely the Coulombic case p = 1.

3.5.2.3.2. The case p = 1

| top | pdf |

In the case [1 \le p \le 3], the lattice sum does not converge absolutely, since the integral in equation (3.5.2.23)[link] diverges. Thus extra conditions are needed and the order of summation must be specified. We restrict discussion to the case p = 1, that is to a crystal of point charges interacting via Coulomb's law. We consider a large but finite crystal, taking the limit as it grows larger while maintaining a specific shape. To that end, let P be a closed, bounded region in [\Re^3] centred on the origin, such as a cube, sphere or ellipsoid. For a positive integer K, let [\Omega(P,K)] denote the set of lattice vectors such that [|{\bf n}| / K] is in P, that is [\Omega(P,K) = \lbrace{\bf n}: |{\bf n}| / K \in {\it P}\rbrace]. Charges in the central unit cell U0 interact with each other and with image charges in Un for [{\bf n} \in \Omega(P,K) ]. The electrostatic energy of the central cell in this finite crystal is given by[E_{P,K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) = {\textstyle{{1}\over{2}}} \sum\limits_{{\bf n} \in \Omega(P,K)}{}^{'} \sum\limits_{i=1}^N \sum\limits_{j=1}^N {{q_{i}q_{j}}\over{ | {\bf r}_{i} - {\bf r}_{j} - {\bf n}|}}\eqno(3.5.2.27) ]and we want the limit of [E_{P,K}] as [K \to \infty].

Recalling that [g_p(x) = \exp(-x^2)/x^2] for p = 1, if we define [h_{{\bf m},{\bf r}}({\bf v})] by[h_{{\bf m},{\bf r}}({\bf v}) = {{\exp[-\pi^2 ({\bf m}+{\bf v})^2/\alpha]}\over{|{\bf m}+{\bf v}|^2}}\exp[-2 \pi i ({\bf v}+{\bf m})\cdot {\bf r})], ]then for [{\bf m} \ne {\bf 0}], [h_{{\bf m},{\bf r}}({\bf v}) ] is bounded and smooth with uniformly (in [{\bf m} \ne {\bf 0}]) bounded derivatives. For m = 0 we have [h_{{\bf m},{\bf r}}({\bf v})] = [[1 - 2\pi i ({\bf v}\cdot{\bf r}) - 2\pi^2({\bf v}\cdot{\bf r})^2] / {\bf v}^2 -\pi^2/ \alpha+ \theta({\bf v})], where [\theta({\bf v})] is smooth with bounded derivatives, and with [\theta({\bf 0}) = 0]. Then, using the same arguments that lead to equation (3.5.2.24)[link], we can write for [{\bf r} \ne {\bf 0} ][\eqalignno{\sum_{{\bf n}\in \Omega({\it P},K)} {{1}\over{|{\bf r}+{\bf n}|}} &= \sum_{\bf n} {{{\rm erfc}(\alpha^{1/2}|{\bf r}+{\bf n}|)}\over{|{\bf r} + {\bf n}|}} &\cr &\quad+ {{1}\over{\pi V}} \sum_{{\bf m} \ne {\bf 0}} {{\exp(-\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}}\exp(-2 \pi i {\bf m}\cdot {\bf r})&\cr &\quad- {{\pi}\over{\alpha V}} + {{1}\over{\pi}} H_{{\it P},K}({\bf r}) + {\varepsilon}(K)&\cr &= \varphi_{\rm Ewald}({\bf r}) + {{1}\over{\pi}} H_{{\it P},K}({\bf r}) + {\varepsilon}(K)&(3.5.2.28)} ]where here and below [{\varepsilon}(K)] denotes a quantity (the remainder in the equation) that converges to zero as [K \to \infty] and[\eqalign{H_{P,K}({\bf r}) &= \sum_{{\bf n}\in \Omega({\it P},K)} \int\limits_{U^{*}} {{1 - 2\pi i ({\bf v}\cdot{\bf r}) - 2\pi^2({\bf v}\cdot{\bf r})^2}\over{{\bf v}^2}} \cr &\quad\times \exp(-2\pi i {\bf v}\cdot{\bf n})\,{\rm d}^3{\bf v}.} ]

The Ewald potential [\varphi_{\rm Ewald}({\bf r})], given for [{\bf r} \ne {\bf 0}] by[\eqalign{\varphi_{\rm Ewald}({\bf r}) &= \sum\limits_{\bf n} {{{\rm erfc}(\alpha^{1/2}|{\bf r}+{\bf n}|)}\over{|{\bf r} + {\bf n}|}} - {{\pi}\over{\alpha V}}\cr &\quad+ {{1}\over{\pi V}} \sum\limits_{{\bf m} \ne {\bf 0}} {{\exp(-\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}}\exp(-2 \pi i {\bf m}\cdot {\bf r}),} ]has a couple of noteworthy properties. First, note that since the left-hand side of equation (3.5.2.28)[link] is independent of α, as is [H_{P,K}({\bf r})], [\varphi_{\rm Ewald}({\bf r})] is itself independent of α. Secondly, the average of [\varphi_{\rm Ewald}({\bf r}) ] over the unit cell U is zero. To establish this latter property, first note that for [{\bf m} \ne {\bf 0}], [\textstyle\int_{U} \exp(2 \pi {\bf m} \cdot {\bf r})\,{\rm d}^3{\bf r} = 0 ]. Thus the reciprocal sum in the above equation integrates to zero over U. Next, the integral of −π/(αV) over U is clearly −π/α. Finally, note that[\eqalign{&\sum\limits_{\bf n} \int\limits_{U} {{{\rm erfc}(\alpha^{1/2}|{\bf r}+{\bf n}|)}\over{|{\bf r} + {\bf n}|}} \,{\rm d}^3{\bf r} \cr&\quad= \int\limits_{\Re^3} {{{\rm erfc}(\alpha^{1/2}|{\bf r}|)}\over{|{\bf r}|}} \,{\rm d}^3{\bf r} \cr &\quad= 8 \pi^{1/2} \int\limits_{r = 0}^\infty r \,{\rm d}r \int\limits_{\alpha^{1/2} r}^\infty \exp(-t^2)\,{\rm d}t\cr &\quad= 4 \pi^{1/2} \int\limits_0^\infty \exp(-t^2) \,{\rm d}t \int\limits_0^{t / \alpha^{1/2} } 2 r \,{\rm d}r \cr &\quad= {{4 \pi^{1/2}}\over{\alpha}} \int\limits_{0}^{\infty} t^2 \exp(-t^2) \,{\rm d}t \cr &\quad ={{\pi}\over{\alpha}},} ]completing the proof. The self or Wigner potential [q\zeta] with [\zeta] given by[\eqalign{\zeta &= \lim_{|{\bf r}| \to 0} \left[\varphi_{\rm Ewald}({\bf r}) - {{1}\over{|{\bf r}|}}\right] \cr &= -2 \left(\,{{\alpha}\over{\pi}}\,\right)^{1/2} + \sum\limits_{{\bf n} \ne {\bf 0}} {{{\rm erfc}(\alpha^{1/2}|{\bf n}|)}\over{|{\bf n}|}} \cr &\quad- {{\pi}\over{\alpha V}} + {{1}\over{\pi V}} \sum\limits_{{\bf m} \ne {\bf 0}} {{\exp(-\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}}} ]is the potential of a charge q in a periodic array of images of itself together with a uniform neutralizing plasma throughout the crystal. For cubic unit cells with side L, [\zeta \simeq -2.837297 / L ] (Nijboer & Ruijgrok, 1988[link]). Owing to the above properties, the Ewald potential provides a consistent approach to the treatment of non-neutral unit cells, such as when calculating ionic charging free energies (Hummer et al., 1996[link]; Darden et al., 1999[link]).

If we define the structure factor [S({\bf m})] by[S({\bf m}) = \textstyle\sum\limits_{j=1}^N q_j \exp(2 \pi i {\bf m} \cdot {\bf r}_j) ]and the unit-cell dipole moment D by[{\bf D} = \textstyle\sum\limits_{j = 1}^N q_j {\bf r}_j, ]then assuming the unit cell is neutral we can write the energy [E_{P,K} ] as[\eqalignno{E_{{\it P},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'}\sum_{i=1}^N \sum_{j=1}^N q_i q_j {{{\rm erfc}(\alpha^{1/2}|{\bf r}_j - {\bf r}_i+{\bf n}|)}\over{|{\bf r}_j - {\bf r}_i + {\bf n}|}} &\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} {{\exp(-\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}} S({\bf m})S({-\bf m}) &\cr &\quad+ 2 \pi \sum_{{\bf n}\in \Omega({\it P},K)} \int\limits_{U^{*}} {{({\bf v}\cdot {\bf D})^2}\over{{\bf v}^2}} \exp(-2 \pi i {\bf v} \cdot {\bf n}) \,{\rm d}^3{\bf v} &\cr&\quad- \left(\,{{\alpha}\over{\pi}}\,\right)^{1/2} \sum_{j=1}^N q_j^2 + {\varepsilon}(K). &(3.5.2.29)} ]Thus, recalling equation (3.5.2.15)[link], we can write[E_{P,K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) = E_{\rm Ewald}(\underline{\bf r}^{\bf N}) + J({\bf D},P, K) + {\varepsilon}(K),\eqno(3.5.2.30) ]where[J({\bf D},P,K) = 2 \pi \sum_{{\bf n}\in \Omega({\it P},K)} \int\limits_{U^{*}} {{({\bf v}\cdot {\bf D})^2}\over{{\bf v}^2}} \exp(-2 \pi i {\bf v} \cdot {\bf n}) \,{\rm d}^3{\bf v}. ]Note that if D = 0, [J({\bf D},P,K) = 0] and the Coulomb lattice sum converges absolutely, i.e. the result is the same as the Ewald sum regardless of the shape P, i.e. order of convergence. Expressions for [J({\bf D},P) = \lim_{K \to \infty} J({\bf D},P,K) ] have been derived in some cases by Smith (1981[link]). In particular, if P is the unit sphere [J({\bf D},P) = 2\pi {\bf D}^2 / (3 V) ]. Following Smith's approach we show this by first deriving an integral representation for [J({\bf D},P)] and then evaluating it for the case P is the unit sphere.

First note that for any continuous function f, [\textstyle\sum_{{\bf n}\in \Omega(P,K)} f({\bf n} / K) V / K^3 ] converges to [ \textstyle\int_P f(\rho) \,{\rm d}^3 \rho] as [K \to \infty]. Then, in the above integral for [J({\bf D},P,K)] we can substitute w = Kv to get[\eqalign{J({\bf D},P,K) &= {{2 \pi}\over{V}} \int\limits_{K\cdot U^{*}} {{({\bf w}\cdot {\bf D})^2}\over{{\bf w}^2}}\,{\rm d}^3{\bf w} \cr &\quad\times \sum_{{\bf n}\in \Omega(P,K)} \exp\left(-2 \pi i {\bf w} \cdot {{{\bf n}}\over{K}}\right) {{V}\over{K^3}} \cr &= {{2 \pi}\over{V}} \int\limits_{K\cdot U^{*}} {{({\bf w}\cdot {\bf D})^2}\over{{\bf w}^2}} \,{\rm d}^3{\bf w} \cr &\quad\times\int\limits_P \exp(-2 \pi i {\bf w} \cdot \rho) \,{\rm d}^3{\rho} + {\varepsilon}(K)} ]and thus[J({\bf D},P) = {{2 \pi}\over{V}} \int\limits_{\Re^3} {{({\bf w}\cdot {\bf D})^2}\over{{\bf w}^2}} \,{\rm d}^3{\bf w} \int\limits_P \exp(-2 \pi i {\bf w} \cdot {\rho})\,{\rm d}^3{\rho}. ]

Specializing to the case P = S, the unit sphere, we see by using spherical coordinates for [\rho] with the pole parallel to w that [\textstyle\int_{S} \exp(-2 \pi i {\bf w} \cdot {\rho}) \,{\rm d}^3{\rho} ] is a function only of the norm of w, that is [\textstyle\int_{S} \exp(-2 \pi i {\bf w} \cdot \rho) \,{\rm d}^3{\rho} = g(|{\bf w}|) ]. Then, again using spherical coordinates, this time for w with the pole parallel to D, and noting that [\textstyle\int_{\theta = 0}^\pi \cos^2(\theta)\sin(\theta) \,{\rm d}\theta = {1\over3} \textstyle\int_{\theta = 0}^\pi \sin(\theta) \,{\rm d}\theta ], we have[\eqalignno{J({\bf D},{S}) &= {{2 \pi}\over{V}} \int\limits_{\Re^3} {{({\bf w}\cdot {\bf D})^2}\over{{\bf w}^2}} g(|{\bf w}|)\,{\rm d}^3{\bf w} &\cr &= {{2 \pi {\bf D}^2}\over{V}} \int\limits_{w = 0}^\infty w^2 g(w) \,{\rm d}w &\cr &\quad\times \int\limits_{\theta = 0}^{\pi} \cos^2(\theta) \sin(\theta)\,{\rm d} \theta \int\limits_{\varphi = 0}^{2 \pi} \,{\rm d} \varphi&\cr&= {{2 \pi {\bf D}^2}\over{3 V}} \int\limits_{\Re^3} g(|{\bf w}|) \,{\rm d}^3{\bf w} &\cr &= {{2 \pi {\bf D}^2}\over{3 V}} \int\limits_{S} \,{\rm d}^3{\bf \rho}&\cr &\quad\times \int\limits_{\Re^3} \exp(-2 \pi i {\bf w} \cdot {\bf \rho}) \,{\rm d}^3{\bf w},&(3.5.2.31)} ]where we have switched the order of integration in the last step. Finally, since [\textstyle\int_{\Re^3} \exp(-2 \pi i {\bf w} \cdot \rho) \,{\rm d}^3{\bf w} = \delta_3({\rho}) ], the three-dimensional Dirac delta function, this last double integral equals one, completing the derivation of the Ewald correction term in the case of a large but finite crystal having a macroscopically spherical shape.

3.5.2.4. The surface term

| top | pdf |

Above we have derived the form of the electrostatic energy in the case of a macroscopically spherical crystal in terms of the Ewald sum plus a correction term. In this section, following the approach of van Eijck & Kroon (1997[link]), we develop a surface-integral-based expression for the shape-dependent correction term for more general macroscopic shapes.

Given a closed bounded region P that contains the unit sphere S in its interior and a K > 0, and recalling equations (3.5.2.4)[link] and (3.5.2.27)[link], we can write the electrostatic energy of the central unit cell in the finite crystal as[\eqalign{E_{P,K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= E(U,U) + {\textstyle{{1}\over{2}}} \sum_{{\bf n} \in \Omega(P,K), {\bf n} \ne {\bf 0}} E(U,U_{\bf n}) \cr &= E(U,U) + {\textstyle{{1}\over{2}}} \sum_{{\bf n} \in \Omega({\it S},K), {\bf n} \ne {\bf 0}} E(U,U_{\bf n}) \cr &\quad+ {\textstyle{{1}\over{2}}} \sum_{{\bf n} \in \Omega(P,K) \setminus \Omega({\it S},K)} E(U,U_{\bf n}) \cr &= E_{{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) + [E_{P,K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) - E_{{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace})].} ]

The asymptotic form of [E_{{S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) ] has been developed above. Recalling equation (3.5.2.5)[link], for large K and neutral unit cells, we can approximate [\delta E_{P\setminus{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) = [E_{P,K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) - E_{{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace})] ] by[\eqalign{\delta E_{P\setminus{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{({\bf n} / {K}) \in P \setminus {\it S}} E(U,U_{\bf n}) \cr & \simeq {\textstyle{{1}\over{2}}} \sum_{({\bf n} / {K}) \in P \setminus {\it S}} {{{\bf D}^2 }\over{|{\bf n}|^3}} - 3 {{({\bf D} \cdot {\bf n})^2}\over{|{\bf n}|^5}} \cr &\simeq {{1}\over{2 V}} \int\limits_{{\bf r} \in P \setminus {\it S}} {{{\bf D}^2}\over{|{\bf r}|^3}} - 3 {{({\bf D} \cdot {\bf r})^2}\over{|{\bf r}|^5}} \,{\rm d}^3{\bf r} \cr &\simeq {{1}\over{2 V}} \oint\limits_{{\bf r} \in \partial P \setminus \partial{\it S}} {{{\bf D}\cdot{\bf r}}\over{|{\bf r}|^3}} {\bf D}\cdot \,{\rm d}{\sigma},} ]where we have used Gauss' divergence theorem (Arfken & Weber, 2000[link]) together with[\nabla \cdot \left( {{{\bf D}\cdot{\bf r}}\over{|{\bf r}|^3}} {\bf D}\right) = {{{\bf D} \cdot {\bf D}}\over{|{\bf r}|^3}} - 3 {{({\bf D} \cdot {\bf r})^2}\over{|{\bf r}|^5}}. ]

Since, considering the outward unit normal on S,[\eqalign{{{1}\over{2 V}} \oint\limits_{{\bf r} \in \partial{\it S}} {{{\bf D}\cdot{\bf r}}\over{|{\bf r}|^3}} {\bf D}\cdot \,{\rm d}{\sigma} &= {{{\bf D}^2}\over{2 V}} \int \int \cos^2(\theta) \sin(\theta) \,{\rm d}\theta \,{\rm d}\varphi \cr &= {{2 \pi {\bf D}^2}\over{3 V}}} ]and since from the previous results[\lim_{K \to \infty} E_{{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) = E_{\rm Ewald}(\underline{\bf r}^{\bf N}) + {{2 \pi {\bf D}^2}\over{3 V}}, ]we can now generalize the lattice summation to asymptotic shapes P:[\lim_{K \to \infty} E_{{\it S},K}(\underline{\bf r}^{\lbrace\bf N\rbrace}) = E_{\rm Ewald}(\underline{\bf r}^{\bf N}) + {{1}\over{2 V}} \oint\limits_{{\bf r} \in \partial P} {{{\bf D}\cdot{\bf r}}\over{|{\bf r}|^3}} {\bf D}\cdot {\rm d}{\sigma}. ]

The latter surface integral is referred to as the `surface term' and clearly must equal the quantity [J({\bf D},P)] discussed above. It was derived previously by Deem et al. (1990[link]) by different arguments. Van Eijck and Kroon calculate it explicitly for some specific crystal shapes. They consider, however, that in actual crystal growth, surface charge distributions should develop to cancel the surface-integral contribution (see also Smith, 1981[link]). In contrast, Scheraga and co-workers argue for the importance of considering this term, which they refer to as the `Lorentz–Ewald correction', in crystal-structure prediction (Pillardy et al., 2000[link]), although in this they are disputed by van Eijck & Kroon (2000[link]) (see also the response by Wedemeyer et al., 2000[link]).

The results in this section are all derived assuming the crystal is surrounded by vacuum, or at least by a medium with no dielectric screening capability. In the following section we study the effect of a surrounding continuum dielectric medium, assuming a spherical crystal (to make the calculations tractable).

3.5.2.5. The polarization response

| top | pdf |

In this section, following the approach in Smith (1981[link], 1994[link]), we derive the polarization correction to the electrostatic energy of U due to immersing the large but finite crystal in a uniform dielectric with dielectric constant [\epsilon > 1], assuming the macroscopic shape of the crystal is spherical. We first derive the dielectric reaction potential at a point [{\bf r} \in U] due to the charges [q_1,\ldots,q_N \in U ] and their periodic images in [U_{\bf n}, {\bf n} \in \Omega(P,K) ], assuming the unit cell is neutral, i.e. [q_1 + \ldots + q_N = 0 ], that P = S, the unit sphere, and that the sphere KS is immersed in the dielectric.

To begin with, given a point [{\bf r} \in \Re^3] and a point [{\bf r}_0 \in K S], with spherical polar coordinates [(r,\theta,\varphi) ] and [(r_0, \theta_0,\varphi_0)], respectively, the electrostatic potential at r due to a unit charge at [{\bf r}_0] is given by the solution to Laplace's equation. Let β be the angle between r and [{\bf r}_0], that is [\beta = \cos^{-1}({\bf r}\cdot{\bf r}_0 / |{\bf r}||{\bf r}_0|) ]. Owing to the axial symmetry, the general solution [\psi({\bf r};{\bf r}_0) ] to Laplace's equation in this case can be written[\psi({\bf r};{\bf r}_0) = \textstyle\sum\limits_{l=0}^\infty [a_l r^{\,l} + b_l r^{-(l+1)}] P_l[\cos(\beta)] ,]where [P_l[\cos(\beta)]] are the Legendre polynomials (Böttcher, 1973[link]). The contributions [a_l] can be thought of as due to charges further from the origin than r, while [b_l] are due to charges closer to the origin. A particular solution to Laplace's equation is given by [1 / |{\bf r} - {\bf r}_0|]. Thus, we can write the potential due to the unit charge at [{\bf r}_0] in the presence of the dielectric continuum surrounding KS as[\eqalign{V_1({\bf r}) &= {{1}\over{|{\bf r} - {\bf r}_0|}} + \sum_{l=0}^\infty \left[ A_l r^{\,l} + {{B_l}\over{ r^{\,(l+1)}}} \right] P_l[\cos(\beta)] \quad {\rm for} \ r \le K, \cr V_2({\bf r}) &= \sum_{l=0}^\infty \left[ C_l r^{\,l} + {{D_l}\over{ r^{\,(l+1)}}}\right] P_l[\cos(\beta)] \quad {\rm for} \ r \ge K.} ]We now apply boundary conditions, noting that [P_l] are an orthogonal basis of functions, so that boundary conditions can be applied term by term. Then we have

  • (1) [V_1] must be finite as [r \to 0], implying [B_l = 0] for all l;

  • (2) [V_2] must go to zero as [r \to \infty], implying [C_l = 0] for all l > 0.

We expand [1 / |{\bf r} - {\bf r}_0|] (for [r > r_0]) in Legendre polynomials to allow comparison of [V_1] with [V_2] as [r \to K] from below:[{{1}\over{|{\bf r} - {\bf r}_0|}} = \sum_{l = 0}^\infty {{r_0^{\,l}}\over{r^{\,(l+1)}}} P_l[\cos(\beta)]. ]Using this we get:

  • (1) [V_2 = V_1] for r = K, implying [D_l = r_0^{\,l} + A_l K^{(2l +1)} ];

  • (2) [\epsilon ({{\partial V_2}/{\partial r}}) = ({{\partial V_1}/{\partial r}}) ] for r = K, implying [A_l] = [-\{{{(\epsilon - 1)(l+1)}/[{l + \epsilon(l+1)}}]\}[ {{r_0^{\,l}}/{K^{(2l+1)}}}] ].

Finally, since we need the solution for the potential due to each of the charges in [U_{\bf n}], we invoke the spherical harmonic addition theorem (Arfken & Weber, 2000[link]) applying it to the points r, [{\bf r}_0]:[{{2 l + 1}\over{4 \pi}} P_l[\cos(\beta)] = \sum_{m = -l}^{+l} (-1)^m Y_l^m(\theta,\varphi) Y_l^{-m}(\theta_0,\varphi_0), ]where (Arfken & Weber, 2000[link])[Y_l^m(\theta,\varphi) = (-1)^m \left[{{2l+1}\over{4 \pi}} {{(l-m)!}\over{(l+m)!}}\right]^{1/2} P_l^m[\cos(\theta)] \exp({im\varphi}) ]and the associated Legendre polynomials [P_l^m(x)], [-l \le m \le l ] are defined by Rodrigues' formula,[\eqalignno{P_l^m[\cos(\theta)] &= (-1)^l {{1}\over{2^l l!}} \sin(\theta)^m {{d^{l+m}}\over{d[\cos(\theta)]^{l+m}}}[1 - \cos^2(\theta)]^l ,&\cr &&(3.5.2.32)} ]and[P_l^{m} [\cos(\theta)] = (-1)^m {{(l + m)!}\over{(l - m)!}} P_l^{-m}[\cos(\theta)].\eqno(3.5.2.33) ]In particular, the first-order spherical harmonics are given by[\eqalignno{Y_1^1(\theta,\varphi) &= -\left({{3}\over{8\pi}}\right)^{1/2} \sin(\theta) \exp({i\varphi}), &\cr Y_1^0(\theta,\varphi) &= \left({{3}\over{4\pi}}\right)^{1/2} \cos(\theta) ,&\cr Y_1^{-1}(\theta,\varphi) &= \left({{3}\over{8\pi}}\right)^{1/2} \sin(\theta) \exp({-i\varphi}).&(3.5.2.34)} ]The spherical harmonics satisfy an orthogonality relation:[\textstyle\int\limits_{\theta = 0}^{\pi} \textstyle\int\limits_{\varphi = 0}^{2\pi} [Y_l^m(\theta,\varphi)]^{*} Y_{l'}^{m'} \sin(\theta) \,{\rm d}\theta \,{\rm d}\varphi = \delta_{l,l'} \delta_{m,m'}. ]Then, for [{\bf r} \in U] we express the polarization potential due to the dielectric response to a unit charge at [{\bf r}_0 \in KS] as[\psi_{\rm Pol}({\bf r};{\bf r}_0) = -\sum\limits_{l=0}^\infty \sum\limits_{m=-l}^{+l} C(l,m,\epsilon) {{1}\over{K^{2l+1}}} g_{l,m}({\bf r}) g_{l,-m}({\bf r}_0), ]where[C(l,m,\epsilon) = (-1)^m {{4 \pi}\over{2l+1}} {{(\epsilon - 1)(l+1)}\over{l + \epsilon(l+1)}} ]and[g_{l,m}({\bf r}) = r^l Y_l^m(\theta,\varphi).]

We now examine the total polarization potential[\eqalignno{ \psi_{{\rm Pol},K,S}({\bf r}) &= \sum_{{\bf n}\in \Omega(S,K)} \sum_{j=1}^N q_j \psi_{\rm Pol}({\bf r}\semi{\bf r}_j + {\bf n}) &\cr &= -{{1}\over{V}} \sum_{l=0}^\infty \sum_{m=-l}^{+l} C(l,m,\epsilon) {{1}\over{K^{l-1}}} g_{l,m}({\bf r}) \cdot K {{V}\over{K^3}}&\cr &\quad \times \sum_{{\bf n}\in \Omega(S,K)} \sum_{j=1}^N q_j \left[g_{l,-m}\left({{{\bf r}_j + {\bf n}}\over{K}}\right) - g_{l,-m}\left({{{\bf n}}\over{K}}\right)\right] &\cr &= -{{1}\over{V}} \sum_{l=0}^\infty \sum_{m=-l}^{+l} C(l,m,\epsilon) {{1}\over{K^{l-1}}} g_{l,m}({\bf r}) &\cr &\quad \times \bigg[{\bf D} \cdot \int\limits_S \nabla g_{l,-m}(\rho)\,{\rm d}^3\rho \bigg]+ {\varepsilon}(K),&(3.5.2.35)} ]where D is the unit-cell dipole and we have used unit-cell neutrality in the next to last equation.

The gradient [\nabla] with respect to the usual Cartesian directions, expressed in spherical polar coordinates (Arfken & Weber, 2000[link]) can be re-expressed as[\eqalign{{{\partial}\over{\partial x}} + i {{\partial}\over{\partial y}} &= \sin(\theta) \exp({i \varphi} ){{\partial}\over{\partial r}} + {{1}\over{r}} \cos(\theta)\exp({i\varphi}) {{\partial}\over{\partial \theta}} + {{i \exp({i\varphi})}\over{r \sin(\theta)}} {{\partial}\over{\partial \varphi}} \cr {{\partial}\over{\partial x}} - i {{\partial}\over{\partial y}} &= \sin(\theta) \exp({-i \varphi}) {{\partial}\over{\partial r}} + {{1}\over{r}} \cos(\theta)\exp({-i\varphi}) {{\partial}\over{\partial \theta}}\cr&\quad- {{i \exp({-i\varphi})}\over{r \sin(\theta)}} {{\partial}\over{\partial \varphi}} \cr {{\partial}\over{\partial z}} &= \cos(\theta) {{\partial}\over{\partial r}} - {{1}\over{r}} \sin(\theta) {{\partial}\over{\partial \theta}}.} ]

Using spherical polar coordinates to carry out the integrations over the unit sphere in equation (3.5.2.35)[link], we see after the first partial integration over [\varphi] that[\eqalign{\int\limits_S \left({{\partial}\over{\partial x}} + i {{\partial}\over{\partial y}}\right) g_{l,-m}(\rho) \,{\rm d}^3\rho &= 0 \quad {\rm if} \ m \ne 1, \cr \int\limits_S \left({{\partial}\over{\partial x}} - i {{\partial}\over{\partial y}}\right) g_{l,-m}(\rho) \,{\rm d}^3\rho &= 0 \quad {\rm if} \ m \ne -1 \ \ {\rm and} \cr \int\limits_S {{\partial}\over{\partial z}} g_{l,-m}(\rho) \,{\rm d}^3\rho &= 0 \quad {\rm if} \ m \ne 0.} ]We recognize the factors multiplying [({{\partial}/{\partial r}} )] above as being proportional to [Y_1^1], [Y_1^0] and [Y_1^{-1}], respectively, so that, using the above orthogonality of spherical harmonics, these contribute nothing to the integrals unless l = 1.

Next, examining I given by[\eqalign{I &= \int\limits_{\theta=0}^{\pi} \int\limits_{\varphi = 0}^{2\pi} \left[\cos(\theta)\exp({i\varphi}) {{\partial}\over{\partial \theta}} + {{i \exp({i\varphi})}\over{\sin(\theta)}} {{\partial}\over{\partial \varphi}}\right] \cr&\quad\times g_{l,-1}(r,\theta,\varphi) \sin(\theta)\,{\rm d}\theta \,{\rm d}\varphi,} ]integrating over [\varphi] and then applying integration by parts over [\theta] to the first term, expressing [g_{l,-1}(r,\theta,\varphi) ] using the Rodrigues' formulas (3.5.2.32)[link] and (3.5.2.33)[link], and finally substituting [w = \cos(\theta)] we see that I is proportional to the integral J given by[J = \int\limits_{w=-1}^{1} (1 - w^2) {{d^{l+1}}\over{dw^{l+1}}} (1 - w^2)^l \,{\rm d}w, ]which in turn is zero unless l = 1. Putting the above results together we see that[\int\limits_S \left({{\partial}\over{\partial x}} + i {{\partial}\over{\partial y}}\right) g_{l,-m}(\rho) \,{\rm d}^3\rho = 0 \quad {\rm unless} \ \ l=1, \ \ m = 1. ]Similarly[\int\limits_S \left({{\partial}\over{\partial x}} - i {{\partial}\over{\partial y}}\right) g_{l,-m}(\rho) \,{\rm d}^3\rho = 0 \quad {\rm unless} \ \ l=1, \ \ m = -1. ]Finally,[I = \int\limits_{\theta=0}^{\pi} \int\limits_{\varphi = 0}^{2\pi} \sin(\theta) {{\partial}\over{\partial \theta}} g_{l,0}(r,\theta,\varphi)\sin(\theta)\,{\rm d}\theta \,{\rm d}\varphi ]is proportional to[J = \int\limits_{w=-1}^1 w {{{\rm d}^{l}}\over{{\rm d}w^{l}}} (1 - w^2)^l \,{\rm d}w ]and J is zero unless l = 1. Thus[\int\limits_S {{\partial}\over{\partial z}} g_{l,-m}(\rho) \,{\rm d}^3\rho = 0 \quad{\rm unless}\ \ l=1, \ \ m = 0. ]

Expressing [g_{1,m}] in terms of the explicit forms for [Y_1^m ] given in equations (3.5.2.34)[link], we get[\int\limits_S \left({{\partial}\over{\partial x}} - i {{\partial}\over{\partial y}}\right) g_{1,1}(\rho) \,{\rm d}^3\rho = -(8\pi/3)^{1/2} ]and thus[\textstyle\int\limits_S \nabla g_{1,1}(\rho)\,{\rm d}^3\rho = -(8\pi/3)^{1/2}(1/2,i/2,0). ]Similarly[\textstyle\int\limits_S \nabla g_{1,0}(\rho)\,{\rm d}^3\rho = (4\pi/3)^{1/2}(0,0,1) ]and finally[\textstyle\int\limits_S \nabla g_{1,-1}(\rho) \,{\rm d}^3\rho =(8\pi/3)^{1/2}(1/2,-i/2,0). ]Substituting these into equation (3.5.2.35)[link] we get[\lim_{K \to \infty} \psi_{{\rm Pol},K,S}({\bf r}) = {{-4 \pi}\over{3V}} {{2\epsilon - 2}\over{2 \epsilon+1}} ({\bf r} \cdot {\bf D}). ]

Recalling the shape-dependent term [J({\bf D},S) = {{2 \pi {\bf D}^2}/{3 V}} ] from equation (3.5.2.31)[link] and adding the reaction energy [\textstyle\sum_{j=1}^N q_j \psi_{{\rm Pol},K,S}({\bf r}_j)], we see that the electrostatic energy of the unit cell U inside a large spherical crystal which in turn is immersed in a dielectric continuum with dielectric constant [\epsilon] is given by[\eqalign{E_{S,K,\epsilon}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= E_{\rm Ewald}(\underline{\bf r}^{\lbrace\bf N\rbrace}) + {{2 \pi {\bf D}^2}\over{3 V}} - {{2 \pi}\over{3V}} {{2\epsilon - 2}\over{2 \epsilon+1}} {\bf D}^2 + {\varepsilon}(K)\cr &= E_{\rm Ewald}({\underline{\bf r}}^{\lbrace\bf N\rbrace}) + {{1}\over{2\epsilon+1}}\cdot {{2\pi}\over{V}}{\bf D}^2 + {\varepsilon}(K),} ]where [E_{\rm Ewald}] is given in equation (3.5.2.15)[link]. Note that we recover the usual Ewald sum in the limit as [\epsilon \to \infty ] (`tin-foil' boundary conditions).

The force on a charge [q_j] in the central cell is obtained by taking the gradient of the energy with respect to [{\bf r}_j] for the finite spherical crystal, and then taking the limit as [K \to \infty], or by taking the gradient of the limiting energy above (i.e. the limit of the gradient is the gradient of the limit):[{\bf F}_j = -{\nabla}_j E_{\rm Ewald}(\underline{\bf r}^{\bf N}) - {{1}\over{2\epsilon+1}}\cdot {{4\pi}\over{V}} q_j {\bf D}.\eqno(3.5.2.36) ]

Thus we can find tractable analytic expressions for the energy and force on particles within the central unit cell as long as the asymptotic shape of the crystal immersed in the dielectric medium is spherical. In Section 3.5.2.6[link] below we show that we can also arrive at a tractable result for the thermodynamic pressure under these circumstances, although surprisingly it does not seem to agree with the mechanical pressure as given by Smith (1994[link]), which may reflect some not yet fully understood subtleties in the implementation of periodic boundary conditions.

3.5.2.6. Calculating the pressure

| top | pdf |

We first discuss the scalar pressure, following the developments in Smith (1987[link], 1993[link]). According to the thermodynamic definition, the scalar pressure P is defined by[P = -\left({{\partial A}\over{\partial V}}\right)_{T}, ]where A is the Helmholtz free energy, V is the volume and T is the temperature. In turn, A is given by [A = -(1/\beta) \ln Z_N(V,T) ] where [\beta = 1 / kT] ([k] is Boltzmann's constant) and [Z_N] is the canonical ensemble partition function,[Z_N(V,T) = {{1}\over{N! h^{3N}}} \int \int \exp[ -\beta H_N(\underline{\bf r}^{\lbrace\bf N\rbrace},\underline{\bf p}^{\lbrace\bf N\rbrace}) ]\,{\rm d}\underline{\bf r}^{\lbrace\bf N\rbrace}\,{\rm d}\underline{\bf p}^{\lbrace\bf N\rbrace}, ]where [H_N] is the system Hamiltonian, h is Planck's constant and [\underline{\bf p}^{\lbrace\bf N\rbrace}] denotes the set of momenta [\lbrace {\bf p}_1,\ldots,{\bf p}_N\rbrace] with [{\bf p}_i] conjugate to [{\bf r}_i]. Introducing the scaling relations [{\bf r}_i = V^{1/3} {\bf s}_i], [{\boldpi}_i = V^{2/3} m_i \dot{\bf s}_i ] and [{\bf p}_i = V^{-1/3} {\boldpi}_i], we get[\eqalign{P &= - {{1}\over{[N! h^{3N} Z_N(V,T)]}} \cr &\quad\times \int \int \left[{{\partial}\over{\partial V}} H_N(V^{1/3}\underline{\bf s}^{\lbrace\bf N\rbrace}, V^{-1/3} \underline{\boldpi}^{\lbrace\bf N\rbrace}) \right]_T \cr &\quad\times \exp[ -\beta H_N(V^{1/3}\underline{\bf s}^{\lbrace\bf N\rbrace},V^{-1/3}\underline{\boldpi}^{\lbrace\bf N\rbrace}) ]\,{\rm d}\underline{\bf s}^{\lbrace\bf N\rbrace}\,{\rm d}\underline{\boldpi}^{\lbrace\bf N\rbrace} \cr &= - \left \langle \left[ {{\partial}\over{\partial V}} H_N(V^{1/3}\underline{\bf s}^{\lbrace\bf N\rbrace}, V^{-1/3} \underline{\boldpi}^{\lbrace\bf N\rbrace}) \right] \right \rangle,} ]where [\langle \ \rangle] denotes the expectation or ensemble average. For the above system of N point charges [q_1,\ldots,q_N ] having mass [m_1,\ldots,m_N] at positions [{\bf r}_1,\ldots,{\bf r}_N ] in the central unit cell [U_0] of a large spherical crystal (of radius K) immersed in a dielectric continuum, the Hamiltonian can be written[\eqalign{H_N(\underline{\bf r}^{\lbrace\bf N\rbrace},\underline{\bf p}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{j=1}^N {{{\bf p}_j^2}\over{m_j}} + E_{\rm Ewald}(\underline{\bf r}^{\lbrace\bf N\rbrace}) \cr &\quad+ {{1}\over{2\epsilon+1}}\cdot {{2\pi}\over{V}}{\bf D}^2 + {\varepsilon}(K)} ]and so[\eqalign{&H_N(V^{1/3}\underline{\bf s}^{\lbrace\bf N\rbrace}, V^{-1/3} \underline{\boldpi}^{\lbrace\bf N\rbrace})\cr &\quad={\textstyle {{1}\over{2}} }V^{-2/3} \sum_{j=1}^N {{{\boldpi}_j^2}\over{m_j}} + E_{\rm Ewald}(V^{1/3} \underline{\bf r}^{\lbrace\bf N\rbrace})\cr&\quad\quad+ V^{-1/3} {{2 \pi}\over{2\epsilon+1}} \bigg(\sum_{j=1}^N q_j {\bf s}_j \bigg)^2 + {\varepsilon}(K).} ]The direct- and reciprocal-lattice vectors are scaled as [{\bf n} = V^{1/3} {\boldeta} ] and [{\bf m} = V^{-1/3} {\boldmu}]. The derivative with respect to volume of the Ewald direct sum is straightforward using the chain rule, since if [x = |{\bf r}_j + {\bf n}| = V^{1/3}|{\bf s}_j + {\boldeta}|] then [\partial x/\partial V = x /3V] and [({\rm d}/{\rm d}x) [{\rm erfc}(\alpha^{1/2}x)/x]] = [-[(2/\pi^{1/2}) \exp(-\alpha x^2) / x + {\rm erfc}(\alpha^{1/2}x) / x^2 ] ]. Since m and [{\bf r}_j] scale oppositely with V, [\partial S({\bf m})/ \partial V = 0]. Continuing with the other terms in the Ewald reciprocal sum, and taking the limit as [K \to \infty], we get[\eqalign{P &= {{1}\over{3V}} \left\langle \sum_{j=1}^N {{{\bf p}_j^2}\over{m_j}} \right. + {\textstyle{{1}\over{2}} }\sum_{\bf n}{}^{'}\sum_{i=1}^N \sum_{j=1}^N q_i q_j \cr &\quad \times \left\lbrace{{2}\over{\pi^{1/2}}} \exp\left[-\alpha ({\bf r}_j - {\bf r}_i + {\bf n})^2\right] + {{{\rm erfc}(\alpha^{1/2} |{\bf r}_j - {\bf r}_i + {\bf n}|)}\over{|{\bf r}_j - {\bf r}_i + {\bf n}|}}\right\rbrace \cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} {{\exp(-\pi^2 {\bf m}^2 / \alpha)}\over{{\bf m}^2}}\bigg(1 - {{2\pi {\bf m}^2}\over{\alpha}}\bigg) S({\bf m})S({-\bf m}) \cr &\quad + \left. {{1}\over{2\epsilon+1}}\cdot {{2\pi}\over{V}}{\bf D}^2 \right\rangle.} ]

The scalar pressure has been generalized to the pressure tensor (Brown & Neyertz, 1995[link]) using a more general scaling technique based on the expression for r in terms of unit-cell basis vectors and fractional coordinates, treating the standard Ewald sum without the surface correction or dielectric response term. To accommodate the latter, the energy of the central unit cell within a large ellipsoidal crystal immersed in a dielectric continuum must be studied. In Redlack & Grindley (1972[link], 1975[link]) expressions are provided for the surface correction in ellipsoidal crystals in a vacuum. If this result were combined with the polarization response term (we are not aware of any papers giving this term for an ellipsoidal crystal immersed in a dielectric continuum) it is reasonable to assume that the standard Ewald sum would apply in this case as well, in the limits [K \to \infty] and [\epsilon \to \infty]. The trace of the pressure tensor defined by Nosé & Klein (1983[link]) agrees with the above scalar pressure in the limit [\epsilon \to \infty ].

Interestingly, Smith (1994[link]) derives a scalar pressure that differs from the above. In particular he finds a surface term that remains even in the limit of `tin-foil' boundary conditions, i.e. when [\epsilon \to \infty]. He gives a mechanical definition of the pressure tensor in terms of momentum flux across surfaces within and surrounding the crystal. Periodic boundary conditions are enforced using velocity constraints via the techniques of DeLeeuw et al. (1990[link]). He derives forces on particles within the central unit cell (and on their images in other cells interior to the crystal) that agree in the large K limit with our result above in equation (3.5.2.36)[link]. However, combining the forces with particle positions in the scalar virial expression, after summing over particles and their images, and proceeding to the limit as [K \to \infty] leads to a vacuum surface correction involving the reciprocal sum term at m = 0 that, unlike our thermodynamic pressure expression above, is not cancelled by the dielectric response potential in the limit [\epsilon \to \infty]. To reiterate, his mechanical virial tensor, restricted to the standard Ewald sum part of the final expression, agrees completely with the results of Nosé & Klein (1983[link]). Furthermore, his scalar virial, given by the trace of his virial tensor, has an [\epsilon]-dependent part due to the dielectric response that agrees completely with our results above. However, the trace of the vacuum surface correction tensor within his mechanical pressure tensor does not agree with the equivalent term in our thermodynamic pressure derived above, and thus is not cancelled in the infinite dielectric limit.

References

Arfken, G. B. & Weber, H. J. (2000). Mathematical Methods for Physicists, 5th ed. San Diego: Academic Press.
Brown, D. & Neyertz, S. (1995). A general pressure tensor calculation for molecular dynamics simulations. Mol. Phys. 84, 577–595.
Böttcher, C. J. F. (1973). Theory of Electric Polarization. Amsterdam: Elsevier Science Publishers BV.
Darden, T., Pearlman, D. A. & Pedersen, L. G. (1999). Ionic charging free energies: Spherical versus periodic boundary conditions. J. Chem. Phys. 109, 10921–10935.
DeLeeuw, S. W., Perram, J. W. & Petersen, H. G. (1990). Hamilton's equations for constrained dynamical systems. J. Stat. Phys. 61, 1203–1222.
Deem, M. W., Newsam, J. M. & Sinha, S. K. (1990). The h = 0 term in Coulomb sums by the Ewald transformation. J. Phys. Chem. 94, 8356–8359.
Eijck, B. P. van & Kroon, J. (1997). Coulomb energy of polar crystals. J. Phys. Chem. B, 101, 1096–1100.
Eijck, B. P. van & Kroon, J. (2000). Comment on `Crystal structure prediction by global optimization as a tool for evaluating potentials: role of the dipole moment correction terms in successful predictions'. J. Phys. Chem. B, 104, 8089.
Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593.
Eyges, L. (1980). The Classical Electromagnetic Field. New York: Dover Publications Inc.
Hummer, G. (1996). Electrostatic potential of a homogeneously charged square and cube in two and three dimensions. J. Electrostat. 3, 285–291.
Hummer, G., Pratt, L. R. & Garcia, A. E. (1996). On the free energy of ionic hydration. J. Phys. Chem. 100, 1206–1215.
Kittel, C. (1986). Introduction to Solid State Physics. New York: John Wiley and Sons.
Körner, T. W. (1988). Fourier Analysis. Cambridge University Press.
Nijboer, B. R. A. & Ruijgrok, Th. W. (1988). On the energy per particle in three- and two-dimensional Wigner lattices. Mol. Phys. 53, 361–382.
Nosé, S. & Klein, M. L. (1983). Constant pressure molecular dynamics for molecular systems. Mol. Phys. 50, 1055–1076.
Pillardy, J., Wawak, R. J., Arnautova, Y. A., Czaplewski, C. & Scheraga, H. A. (2000). Crystal structure prediction by global optimization as a tool for evaluating potentials: Role of the dipole moment correction term in successful predictions. J. Am. Chem. Soc. 122, 907–921.
Redlack, A. & Grindley, J. (1972). The electrostatic potential in a finite ionic crystal. Can. J. Phys. 50, 2815–2825.
Redlack, A. & Grindley, J. (1975). Coulombic potential lattice sums. J. Phys. Chem. Solids, 36, 73–82.
Smith, E. R. (1981). Electrostatic energy in ionic crystals. Proc. R. Soc. London Ser. A, 373, 27–56.
Smith, E. R. (1994). Calculating the pressure in simulations using periodic boundary conditions. J. Stat. Phys. 77, 449–472.
Smith, W. (1987). Coping with the pressure: How to calculate the virial. CCP5 Inf. Q. 26, 43–50.
Smith, W. (1993). Calculating the pressure. CCP5 Inf. Q. 39, 1–7.
Wedemeyer, W. J., Arnautova, Y. A., Pillardy, J., Wawak, R. A., Czaplewsky, C. & Scheraga, H. A. (2000). Reply to `Comment on Crystal structure prediction by global optimization as a tool for evaluating potentials: role of the dipole moment correction terms in successful predictions'. J. Phys. Chem. B, 104, 8090–8092.








































to end of page
to top of page