Tables for
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 3.5, pp. 458-481
doi: 10.1107/97809553602060000772

Chapter 3.5. Extensions of the Ewald method for Coulomb interactions in crystals

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:

The electrostatic energy per unit cell of a large but finite crystal of point charges that is immersed in a dielectric medium can be efficiently calculated using the Ewald-sum technique (assuming a neutral unit cell). The form of this sum is derived here in a manner that generalizes to other interactions, such as long-range dispersion interactions, that depend on higher inverse powers of distance. The Ewald method is then generalized to finite crystals of model charge densities using Hermite Gaussian basis functions. Finally, efficient methods to accelerate the Ewald sum for these charge densities to near-linear scaling (i.e. scaling as N log N with system size N) are derived and tested.

3.5.1. Introduction

| top | pdf |

High-precision single-crystal X-ray structural analysis of small organic molecules, yielding the space group, the unit-cell parameters and the fractional coordinates of the atoms making up the molecule(s) in the asymmetric unit, has become a routine matter as long as crystals of sufficient quality can be obtained. The thermodynamic stability of the crystal, as described by the enthalpy of sublimation ΔHsub, can also be determined experimentally (although not always to high precision). Theoretical models for calculating intermolecular interaction energies can be used to connect the crystal structure to the molar enthalpy of sublimation using the relationship[\Delta H_{\rm sub} (0\ {\rm K}) = -E_{\rm lattice},]where the lattice or packing energy Elattice is the total (molar) intermolecular interaction energy between all the molecules in the crystal, which are treated as rigid entities with zero-point energies of intra- and intermolecular vibrations neglected. Connection to experimentally accessible heats of sublimation at higher temperatures involves thermodynamic corrections. Methods for calculating thermodynamic quantities of solids are discussed in Gavezzotti (2002a[link]) and (in more detail) in Frenkel & Smit (2002[link]).

Thus, given a parameterized intermolecular potential-energy function, or if computationally affordable a first-principles approach such as density-functional theory (or preferably, when it becomes feasible for crystals, a good-quality post-Hartree–Fock potential-energy surface that describes dispersion interactions), one can sum the intermolecular energies to obtain the lattice energy as a function of the above parameters defining the crystal structure. If such an energy function is used together with a method for systematic search of the crystal-structure parameters, one could in principle predict the minimum-lattice-energy crystal structure for a rigid organic molecule. To extend this approach to flexible molecules one would need to minimize the sum of the intramolecular energy plus the lattice energy. If the experimental crystal structure corresponds to the thermodynamic minimum-energy structure (i.e. it is not a metastable state determined by crystal-growth kinetics), one could in principle predict the experimental crystal structure of an organic compound through this minimization protocol. Moreover, one could ideally predict the additional metastable forms of the crystal.

Prediction of the structure of crystals of an organic molecule from its molecular structure is a difficult problem that has been compared to the protein folding problem (Dunitz, 2003[link]; Dunitz & Scheraga, 2004[link]). Like the protein folding problem, a solution of the crystal prediction problem has significant practical ramifications. A compound is often polymorphic, that is it has more than one crystal structure, and it may be difficult to characterize the conditions under which a particular crystal structure is formed (Dunitz & Bernstein, 1995[link]). Polymorphs may have very different physical properties. An obvious example is diamond versus graphite, but other commercially important examples include food additives, various solid forms of explosives and the bio­availability of various forms of a drug such as ritonavir (Chemburkar et al., 2000[link]). A method for predicting the possible crystal structures of the compound, and ideally for predicting the dominant crystal structure given the experimental conditions, would thus be very valuable. Note that due to the subtleties of crystallization, the lowest-free-energy polymorph at given temperature and pressure may not be the likeliest to form. The kinetics of growth of microcrystals may largely determine which low-energy polymorph appears (Dunitz, 2003[link]). However, it is generally agreed that accurate calculation of the relative free energy of polymorphs is a prerequisite for predicting crystal structures.

To assess progress towards solving this latter problem, a series of blind tests of crystal-structure prediction has been undertaken (Day, Motherwell, Ammon et al., 2005[link]). The results of these tests have highlighted the need for continued improvements in sampling methods and intermolecular energy potentials. Since extensive sampling of the crystal-structure parameters is necessary [between 104 and 105 starting structures, each followed by parameter minimizations (Price & Price, 2005[link])], there is a trade-off in the computational cost versus accuracy of the intermolecular energy functions used. Calculating the work of transforming between polymorphs is yet more ambitious in terms of sampling. Consequently empirical force fields are likely to be needed for the near term at least.

In the remainder of the introduction we outline some of the approaches to empirical potentials used in the calculation of the lattice energy, and then, motivated by these developments, discuss techniques for efficient summation of the electrostatic and other slow-decaying interaction terms that occur in these potential functions.

Methodological developments in the intermolecular force fields used in crystal-structure prediction from early times to the present state of the art have been reviewed (Price & Price, 2005[link]). Until recently, these force fields were made up of atom–atom interactions. The earliest involved only repulsion and dispersion, usually in the `exp-6' form[U^{MN} = U_{\rm rep}^{MN} + U_{\rm disp}^{MN} = \textstyle\sum\limits_{i \in M,j \in N} A_{ij} \exp(-B_{ij}r_{ij}) - C_{i,j}r_{ij}^{-6}, ]where UMN denotes the intermolecular potential energy between molecules M and N and rij is the distance between atoms [i \in M] and [j \in N ]. Sometimes the exponential form in the above equation is replaced by a simpler power law, as in the Lennard–Jones potential. As was pointed out by Dunitz (2003[link]), in comparison with more sophisticated force fields, this repulsion–dispersion form readily allows analysis of the significance of particular atom–atom interactions, since the interactions are short-ranged and thus can be localized. That is, the r−6 form of the attractive dispersion energy means that interaction energies are halved for every 12% increase in distance. In contrast, introduction of long-range Coulombic interactions not only entails subtleties in lattice summation (the subject of this contribution), but greatly complicates the assignment of `key' atom–atom interactions. Gavezzotti and Fillipini systematically explored the use of the exp-6 potential in fitting organic crystal structures with and without hydrogen-bond interactions (Gavezzotti & Fillipini, 1994[link]). They were surprisingly successful in accounting for weak hydrogen bonding in this way, but selective use of point charges improved the directionality of the potential. Earlier, Williams derived exp-6 parameters for the atoms C, H, N, O, Cl, F and polar H for use in organic crystal structures, but found it necessary (Williams & Cox, 1984[link]) to supplement these with selected point charges, both atomic and at off-atom sites. Price and co-workers have used a more sophisticated electrostatic model based on distributed multipole analysis (DMA) with multipoles up to hexadecapole order positioned at atom centres (Stone, 1996[link]). A similar approach has been adopted by Mooji et al., with the addition of dipole polarizabilities to account for intermolecular induction energy (Mooji et al., 1999[link]). The AMOEBA force field has similar capabilities (Ren & Ponder, 2003[link]).

Recently, a systematic study (Day et al., 2004[link]; Day, Motherwell & Jones, 2005[link]) of the effect of the electrostatic model on crystal-structure prediction for rigid organic molecules was carried out. The exp-6 parameters from the W99 force field (Williams, 2001[link]) were used in conjunction with point charges derived either from a bond-increment model or electrostatic potential (ESP) fit, or with a DMA model as used by Price and co-workers. Importantly, they found a systematic and significant improvement in prediction success going from the simple bond-increment model to the ESP charges to point multipoles. Similar results were found by Broderson et al. (2003[link]). This is a positive indicator for the use of theory in attacking this difficult prediction problem. On the other hand, Broderson et al. found that more sophisticated electrostatic models did not help in predicting the packing of flexible molecules. As discussed by Price & Price (2005[link]), crystal-packing energies of flexible molecules are more problematic for at least two reasons. First the intramolecular energy must be fully compatible with the intermolecular energy in order that energy deformation from the gas-phase minimum be on the same scale as the concomitant gain in packing energy. Secondly, the conformational dependence of the electrostatic model must be accurately accounted for, or any gains from a more sophisticated model are immediately negated. Existing force fields are unable to meet these challenges, which may explain the lower success rate for predicting crystal structures of flexible molecules. Price has developed a methodology for flexible molecules wherein intramolecular conformations are systematically explored. For each conformer studied, the gas-phase intramolecular energy and DMA electrostatic model are calculated using quantum-chemical methods. The optimal packing of that conformer is then obtained by treating it as a rigid molecule, and the resulting minimum lattice energy is added to its intramolecular energy to give its total energy. The total energy for that conformer is then compared with that of the other conformers to predict the optimal conformation and packing of the flexible molecule. This approach has proved successful in a number of cases, for example with the successful prediction of a novel second polymorph of 1-hydroxy-7-azabenzotriazole (Nowell et al., 2006[link]).

The above electrostatic models are all derived from gas-phase quantum-chemical calculations. A different line of development has been the parameterization of electrostatic models from condensed-phase experimental or theoretical electron densities, which should better represent the polarization due to the crystal environment. This is accomplished using aspherical multipole refinement of experimental X-ray charge densities (Destro et al., 2000[link]) or of theoretical charge densities from quantum-chemical calculations carried out under the space-group symmetry (Dovesi et al., 2005[link]). Crystal structure analysis is usually carried out assuming the independent-atom model (IAM), where the individual atomic densities are modelled as those of the spherically symmetric isolated atoms, with deviations in the experimental density being fitted to anisotropic temperature factors modelling the thermal motion of the nuclei. However, the IAM breaks down for data of sufficiently high quality, and the residual or deformation density is then fitted by a linear combination of atom-centred density basis elements resembling the well known Slater atomic basis sets from quantum chemistry (i.e. solid harmonics times a decaying exponential in the distance r from the atom). Technical details are provided in Coppens (1997[link]) and Chapter 1.2[link] of this volume. The theoretical density can be similarly approached as a sum of pro-molecule density (superposition of isolated atomic densities in the molecular geometry) plus deformation density modelled by a sum of atomic Slater-like density elements. Until recently, the theoretical densities were thought to be of lower quality than the experimental densities, but continued improvements in the level of theory that can be applied have resulted in much closer agreement (Coppens & Volkov, 2004[link]). Coppens and co-workers as well as others (Spackman et al., 1988[link]) have applied these modelled densities in studies of intermolecular electrostatic energies. The Slater-like functions are typically difficult to work with in this context, so multipole approximations are used except at close range, where numerical integration is used instead (Volkov et al., 2004[link]).

A similar density-based approach using gas-phase theoretical density was developed previously (Gavezzotti, 2002b[link]). Gavezzotti splits the density into numerical `pixels' (small cubes of volume) within which the density is considered constant, allowing in principle exact numerical calculation of intermolecular electrostatic interactions. He sums the pixel–pixel interactions discretely, using a hierarchical approach to increase efficiency. The electrostatic interactions are supplemented by induction and repulsion terms, also modelled by pixel-based interactions. These developments, in particular the large size of the positive and negative components of the electrostatic interactions treated in terms of actual electron density, have caused a reassessment of the idea that intermolecular interactions can be accurately decomposed even in principle into localized atom–atom interactions (Dunitz & Gavezzotti, 2005[link]).

Cisneros et al. have also begun implementing an intermolecular potential based on the explicit modelling of gas-phase theoretical electron density, using the methodology of auxiliary basis set fitting (Cisneros et al., 2005[link]). The density is expanded into atom and off-atom sites using Gaussian density functions similar to those used as atomic basis sets in quantum chemical codes. More recently, other components of the potential such as induction, charge transfer and non-isotropic site–site repulsion based on density overlap have been implemented (Piquemal et al., 2006[link]) and efficient calculation of terms in periodic boundary conditions has also been implemented (Cisneros et al., 2006[link]), generalizing previous developments for point multipoles (Sagui et al., 2004[link]).

Summarizing the above developments, it appears that accurate evaluation of crystal lattice energies necessitates the use of complex electrostatic interactions beyond the spherical-atom approximation, and ultimately including an account of penetration effects through explicit modelling of electron density. The electrostatic lattice energies are differences between large positive and negative contributions that are long-ranged, necessitating an accurate treatment of their lattice summation. Dispersion interactions, modelled by higher-order inverse powers of distance, are short-ranged but their calculation can be greatly accelerated through lattice-summation techniques.

Ewald summation in ideal infinite crystals, for interactions that depend on inverse powers of site–site distance, is covered in Chapter 3.4[link] of this volume, where there is also a numerical demonstration of the advantages of the Ewald approach over direct numerical summation. In this contribution we focus on the lattice-summation problem for large but finite crystals, highlighting the role of summation order and leading to a derivation of the surface term and polarization response. We also consider extensions to Coulombic interactions between Gaussian density elements, both spherical and Hermite (multipolar), and methods to further accelerate the lattice summation to (almost) linear scaling.

These developments are organized into three sections. The first section treats lattice summation of point charges. We begin with definitions and an explanation of the origin of the conditional convergence problem for these sums, anticipating the nature of the results to follow. Following that we provide for the benefit of less mathematically oriented readers the standard physics-based derivation of the Ewald sum for point charges, minus any accounting for the order of convergence or resulting surface term. This is followed by a more careful derivation focused on the order of summation, including results for higher-order inverse powers as well as the Coulomb interactions, finishing with a discussion of the surface term, the polarization response to an external dielectric and the pressure calculation. Although pitched at a higher mathematical level, we have tried to provide many details in the derivations unavailable elsewhere for the benefit of the reader. The alert reader will note that the key identity, equation ([link], that separates 1/rn into short-ranged and long-ranged components, is the same as that used in Chapter 3.4[link] of this volume, highlighting the essential similarity of our approach. In the second section we generalize these point-charge results to lattice sums of interacting spherical Gaussian and Hermite Gaussian-based densities. The third and final section provides a discussion of two methods that greatly increase the numerical efficiency of these lattice sums in the case of large unit cells by utilizing the fast Fourier transform. Again we have tried to provide mathematical details not available elsewhere.

3.5.2. Lattice sums of point charges

| top | pdf |

We begin by discussing an idealized infinite crystal C, made up of point charges. 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( ]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( ]

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}|}},&(} ]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 ([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( ]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.&(} ]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. 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( ]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( ]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 ([link] and the sum is over all general lattice vectors n, defined in equation ([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( ]where the sum is over all reciprocal-lattice vectors m, defined in equation ([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),&(} ]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( ]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( ]Using equations ([link] and ([link] we can rewrite equation ([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 ([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( ]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 ([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[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( ]

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[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 ([link], centred at [{\bf r}_i+{\bf n} ]. Using equation ([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 ([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 ([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( ]

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 ([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,&(} ]where as above the prime indicates that terms with i = j and n = 0 are omitted. Ewald sum: more complete derivation

| top | pdf |

The derivation in Section[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},&(} ]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 ([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,&(} ]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}, &(} ]which is the Fourier transform of the Gaussian, as derived above within equation ([link]. In equation ([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 ([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},&(} ]where[f_p(x) = {{1}\over{\Gamma(p/2)}} \int\limits_{x}^\infty 2 s^{\,p- 1} \exp(-s^2)\,{\rm d}s \eqno( ]and[g_p(x) = x^{\,p-3} \int\limits_x^\infty 2 s^{\,2-p} \exp(-s^2)\,{\rm d}s.\eqno( ]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 ([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 ([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 ([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)}}.&(} ]

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 ([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}}.} ] 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( ]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 ([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}).&(} ]

Using the limiting behaviour given in ([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}).&(} ]

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 ([link] and ([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( ]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. 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 ([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( ]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 ([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)&(} ]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 ([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). &(} ]Thus, recalling equation ([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( ]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},&(} ]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. 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 ([link] and ([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 ([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). 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 &&(} ]and[P_l^{m} [\cos(\theta)] = (-1)^m {{(l + m)!}\over{(l - m)!}} P_l^{-m}[\cos(\theta)].\eqno( ]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}).&(} ]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),&(} ]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 ([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 ([link] and ([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 ([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 ([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 ([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 ([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( ]

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[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. 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 ([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.

3.5.3. Generalization to Gaussian- and Hermite-based continuous charge distributions

| top | pdf |

In this section we develop the generalization of Ewald sums to interacting molecular densities modelled by linear combinations of spherical or Hermite Gaussians. In Coppens (1997[link]) the electrostatic potential and energy for crystal densities modelled by linear combinations of Slater-like density elements is discussed. As with quantum-chemical programs, there are trade-offs in the choice of density elements. The use of Slater-like densities provides accuracy with the use of fewer elements, but the Gaussian and Hermite Gaussian elements are far more convenient analytically. In particular the Ewald summation, as well as the fast-Fourier-transform-based methods discussed in Section 3.5.4[link], generalize straightforwardly to these elements. We note ahead of time that the case of lattice sums of point multipoles are obtained by specializing the results for Hermite Gaussians, just as the results for spherical Gaussians specialize to the above lattice summation results for point charges.

A somewhat different approach to Ewald summation over Hermite Gaussian elements, focused on application to periodic density-functional calculations, is given in Trickey et al. (2004[link]). The continuous fast multipole method (Challacombe et al., 1997[link]) provides yet a different but highly efficient approach to lattice summation of continuous charge densities. Lattice sums of interacting spherical Gaussians

| top | pdf |

We begin by discussing the electrostatic interaction energy between a normalized spherical Gaussian charge distribution [\rho_1], centred at a point [{\bf R}_1 \in U] and having exponent [\theta_1], and a second normalized spherical Gaussian charge distribution [\rho_2], centred at [{\bf R}_2 \in U ] and having exponent [\theta_2], together with the latter's periodic images, centred at [{\bf R} + {\bf n}], [{\bf n} \in \Omega(S,K) ]. Their interaction energy can be written[\eqalignno{E_{12} &= \textstyle\sum\limits_{{\bf n} \in \Omega(S,K)} \textstyle\int\limits_{\Re^3}\textstyle\int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n})&\cr &\quad\times {{1}/{|{\bf r}_2 - {\bf r}_1|}}\,{\rm d}^3 {\bf r}_1 \,{\rm d}^3 {\bf r}_2,&(} ]where [\rho_j({\bf r}) = (\theta_j/\pi)^{3/2} \exp(-\theta_j {\bf r}^2) ], j = 1, 2. We now recall equation ([link], applying its two right-hand components separately, to get [E_{12} = E_{12}^{\,\rm d} + E_{12}^{\,\rm r} ], where[\eqalign{E_{12}^{\,\rm d} &= \sum\limits_{{\bf n} \in \Omega(S,K)} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n}) \cr &\quad\times {{{\rm erfc}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2} ]and[\eqalign{E_{12}^{\,\rm r} &= {{1}\over{\pi}} \sum\limits_{{\bf n} \in \Omega(S,K)} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n}) \cr &\quad \times \int\limits_{\Re^3} {{\exp(-\pi^2 {\bf u}^2 / \alpha)}\over{{\bf u}^2}} \exp[-2 \pi i {\bf u}\cdot ( {\bf r}_2 - {\bf r}_1)]\,{\rm d}^3 {\bf u} \, {\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2.} ]

The direct pair energy [E_{12}^{\,\rm d}] involves the damped Coulomb interaction, which has been discussed by Gill & Adamson (1996[link]). We can calculate [E_{12}^{\,\rm d}] using the following standard results from quantum chemistry, which are derived in Helgaker et al. (2000[link]): Let [\rho_p({\bf r})], [\rho_q({\bf r}) ] be two normalized spherical Gaussian charge distributions having exponents p and q and centred at points P and Q, i.e. [\rho_p({\bf r}) = (p / \pi)^{3/2} \exp[-p({\bf r} - {\bf P})^2] ] and [\rho_q({\bf r})] = [(q / \pi)^{3/2} \exp[-q({\bf r} - {\bf Q})^2] ]. Then the electrostatic potential [V_p({\bf R})] at a point R due to the charge distribution [\rho_p({\bf r})] is given by[V_p({\bf R}) = \left({{4 p}\over{\pi}}\right)^{1/2} F_0[p ({\bf P} - {\bf R})^2], \eqno( ]where [F_0] is the zeroth-order Boy's function, with [F_0(x) =] [ \textstyle\int_0^1 \exp(-x t^2) \,{\rm d}t ] for x > 0. The Coulomb interaction energy [E_{pq}] between [\rho_p] and [\rho_q] is given by[\eqalignno{E_{pq} &= \int\limits_{\Re^3} \int\limits_{\Re^3} {{\rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q})}\over{|{\bf r}_1 - {\bf r}_2|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr &= \int\limits_{\Re^3} V_p({\bf r}_2) \rho_q({\bf r}_2 - {\bf Q}) \,{\rm d}^3{\bf r}_2&\cr&= \left({{4 \mu_{pq}}\over{\pi}}\right)^{1/2} F_0[\mu_{pq} ({\bf P} - {\bf Q})^2],&(} ]where [ 1 / \mu_{pq} = 1 / p + 1 / q]. The Boy's function can be related to the error function:[{{2}\over{\pi^{1/2}}} F_0(x^2) = {{{\rm erf}(|x|)}\over{|x|}}.\eqno( ]Substituting equation ([link] into ([link], we see that[V_p({\bf R}) = {{{\rm erf}(p^{1/2} |{\bf R} - {\bf P}|) }\over{ |{\bf R} - {\bf P}|}},\eqno( ]a result which we had derived above by another method [see equation ([link]]. From equations ([link] and ([link], substituting [\alpha] for p and [\rho_p ] for [\rho_q], we have for any point [{\bf r} \in \Re^3][\int\limits_{\Re^3} {{ {\rm erf}(\alpha^{1/2} |{\bf r} - {\bf r}_1|)}\over{|{\bf r} - {\bf r}_1|}} \rho_p({\bf r}_1 - {\bf P})\,{\rm d}^3{\bf r}_1 = {{{\rm erf}(\mu_{p\alpha}^{1/2}\ |{\bf r} - {\bf P}|)}\over{|{\bf r} - {\bf P}|}}, ]where [1/\mu_{p \alpha} = 1 / p + 1 / \alpha]. Iterating this result with [\mu_{p \alpha}] replacing [\alpha] and [\rho_q ] replacing [\rho_p] we get the damped interaction energy between [\rho_p] and [\rho_q],[\eqalignno{ E_{pq}^{\,\rm d} &= \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) &\cr &\quad\times {{{\rm erfc}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2&\cr =& \int\limits_{\Re^3} \int\limits_{\Re^3} {{\rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q})}\over{|{\bf r}_1 - {\bf r}_2|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr&\quad- \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) &\cr &\quad\times {{{\rm erf}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2&\cr &= {{{\rm erf}(\mu_{pq}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} - {{{\rm erf}(\mu_{pq\alpha}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} &\cr &= {{{\rm erfc}(\mu_{pq\alpha}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} - {{{\rm erfc}(\mu_{pq}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}},&(} ]where [1 / \mu_{pq} = 1 / p + 1 / q] and [1 / \mu_{pq\alpha} = 1 / p + 1 / q + 1 / \alpha ]. Notice that in the limit [\alpha \to \infty], [E_{pq,\alpha} \to 0 ] as expected, since [{\rm erfc}(x) \to 0] as [x \to \infty ], and that in the limit [p \to \infty], [q \to \infty], [E_{pq,\alpha} \to {\rm erfc}(\alpha^{1/2}|{\bf P} - {\bf Q}|) / |{\bf P} - {\bf Q}| ], the damped interaction between two unit point charges at P and Q, as in the standard Ewald direct sum.

To calculate the reciprocal pair energy [E_{12}^{\,\rm r}] we first apply the first right-hand component of equation ([link] together with the Fourier transforms of the normalized Gaussians [\rho_p ] and [\rho_q] from equation ([link] to obtain[\eqalignno{E_{pq}^{\,\rm r} &= {{1}\over{\pi}} \int\limits_{\Re^3} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) {{\exp(-\pi^2 {\bf u}^2 / \alpha)}\over{{\bf u}^2}}&\cr &\quad \times \exp[-2 \pi i {\bf u}\cdot ({\bf r}_1 - {\bf r}_2)] \,{\rm d}^3 {\bf u} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr&= {{1}\over{\pi}} \int\limits_{\Re^3} \int\limits_{\Re^3} \int\limits_{\Re^3} {{1}\over{{\bf u}^2}} \exp[-\pi^2 ({\bf u}^2/\alpha + {\bf v}^2/p + {\bf w}^2/q)] &\cr&\quad\times \exp(2 \pi i {\bf v} \cdot {\bf P}) \exp(2 \pi i {\bf w} \cdot {\bf Q}) \,{\rm d}^3{\bf u} \,{\rm d}^3{\bf v}\,{\rm d}^3{\bf w}&\cr&\quad\times \int\limits_{\Re^3} \exp[-2 \pi i ({\bf u}+{\bf v})\cdot{\bf r}_1]\,{\rm d}^3{\bf r}_1&\cr &\quad\times \int\limits_{\Re^3}\exp[2 \pi i ({\bf u}-{\bf w})\cdot{\bf r}_2] \,{\rm d}^3{\bf r}_2 &\cr &= {{1}\over{\pi}} \int\limits_{\Re^3} {{ \exp(-\pi^2 {\bf u}^2 / \mu_{pq\alpha}) }\over{{\bf u}^2}} \exp[-2 \pi i {\bf u} \cdot ({\bf P} - {\bf Q})] \,{\rm d}^3{\bf u},&\cr &&(} ]where we have used the fact that [\textstyle\int \exp[-2 \pi i ({\bf u}+{\bf v})\cdot{\bf r}_1] \,{\rm d}^3{\bf r}_1] = [\delta^3({\bf u}+{\bf v}) ], with [\delta^3] the three-dimensional delta function, and similarly for the integral over [{\bf r}_2]. Note that in the limit [p \to \infty ], [q \to \infty], [E_{pq}^{\,\rm r}] reduces to the equivalent expression for unit point charges at P and Q, and in the limit [\alpha \to \infty], [E_{pq}^{\,\rm r}] reduces to the Fourier-space expression for the direct Coulomb inter­action between the Gaussians [\rho_p ] and [\rho_q].

Recalling the developments leading to equation ([link], and using equations ([link] and ([link], we can write the inter­action energy [E_{12}] from equation ([link] as[\eqalignno{E_{12} &= \sum_{\bf n} {{{\rm erfc}(\mu_{12\alpha}^{1/2} |{\bf R}_{12} - {\bf n}|) - {\rm erfc}(\mu_{12}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over{|{\bf R}_{12} - {\bf n}|}}&\cr &\quad+ {{1}\over{\pi V}} \sum_{{\bf m} \ne {\bf 0}} {{ \exp(-\pi^2 {\bf m}^2 / \mu_{12\alpha}) }\over{{\bf m}^2}}\exp[-2 \pi i {\bf m} \cdot ({\bf R}_{12})] &\cr&\quad- {{\pi}\over{\mu_{12\alpha} V}} + {{1}\over{\pi}} H_{S,K}({\bf R}_{12}) + {\varepsilon}(K),&(} ]where [{\bf R}_{12} = {\bf R}_1 - {\bf R}_2], [1 / \mu_{12} = 1 / \theta_1 + 1 / \theta_2 ] and [1 / \mu_{12\alpha}= 1 / \theta_1 ] [+ \ 1 / \theta_2 + 1 / \alpha ].

Note that although the energy of a finite-exponent Gaussian charge distribution interacting with itself is finite, when considering a generalization of lattice sums involving point charges it is necessary to remove the direct Coulomb interaction of the Gaussian with itself. From equations ([link] and ([link], if p = q, then as [{\bf Q} \to {\bf P}], [E_{pq}^{\,\rm d} - E_{pq} \to -2 (\mu_{pq\alpha} / \pi)^{1/2} ]. We can now generalize equation ([link]. Let [\rho_1, \ldots,\rho_N ] be normalized spherical Gaussian charge distributions with Gaussian exponents [\theta_1,\ldots,\theta_N], centred at [\lbrace{\bf R}_1,\ldots,{\bf R}_N\rbrace = \underline{\bf R}^{\lbrace\bf N\rbrace}\in U ] (point charges are included by taking [\theta \to \infty]), and let [q_1,\ldots,q_N] satisfy [q_1+\ldots+q_N = 0] (neutral unit cell). Then the energy of the central unit cell [U] within a large spherical crystal, due to interactions of the Gaussian charge distributions [q_i \rho_i] with each other and with their periodic images within the crystal, is given by[\eqalign{E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'}\sum_{i,j=1}^N q_i q_j \left \lbrace {{{\rm erfc}(\mu_{ij\alpha}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right. \cr &\quad- \left. \ {{{\rm erfc}(\mu_{ij}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right \rbrace \cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} \sum_{i,j=1}^N q_iq_j {{ \exp(-\pi^2 {\bf m}^2 / \mu_{ij\alpha})}\over{{\bf m}^2}} \cr &\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij}) \cr &\quad- {{\pi}\over{2 V}} \sum_{i,j = 1}^N {{q_i q_j}\over{\mu_{ij\alpha}}} - \sum_{i=1}^N q_i^2 \left({{\mu_{ii\alpha}}\over{\pi}}\right)^{1/2} \cr &+\quad {{2 \pi {\bf D}^2}\over{3 V}} + {\varepsilon}(K),} ]where [{\bf R}_{ij} = {\bf R}_i - {\bf R}_j], [1 / \mu_{ij} = 1 / \theta_i + 1 / \theta_j ], [1 / \mu_{ij\alpha}= 1 / \theta_i +  1 / \theta_j] [ +\ 1 / \alpha] and as before [{\bf D} = q_1{\bf R}_1 +\ldots + q_N {\bf R}_N] is the unit-cell dipole.

Note that from the above derivation, the [\alpha] in [\mu_{ij\alpha} ] is allowed to be different for each pair ij. One choice that leads to an efficient algorithm (particularly when combined with the fast-Fourier-transform-based methods discussed below) is to separate the Gaussian charge distributions [q_j \rho_j] into compact and diffuse sets, based on their exponents [\theta_j], and choose [\alpha] to be infinite for ij pairs where at least one of the two Gaussians is diffuse. These pair interactions are then evaluated entirely in reciprocal space. Next choose [\alpha] so that [\mu_{ij\alpha}] is constant for all compact pairs ij. More specifically, given [\beta > 0], a Gaussian [q_j \rho_j] is classified as compact if [\theta_j \ge 2 \beta ] and diffuse otherwise. We denote these cases as [j \in C] and [j \in D], respectively. Then for [i,j \in C], choose [\alpha ] so that [1/\mu_{ij\alpha}] = [1/\theta_i + 1/\theta_j + 1/\alpha = 1/\beta ]. Otherwise, choose [\alpha] to be infinite. Then the Coulomb energy of the central unit cell [E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) ] can be re-expressed as[\eqalignno{E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'}\sum_{(i,j)\in C\times C} q_i q_j \left \lbrace {{{\rm erfc}(\beta^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right. &\cr&\quad- \left. \ {{{\rm erfc}(\mu_{ij}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right \rbrace &\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}}\ \sum_{(i,j)\in C\times C} q_iq_j \ {{ \exp(-\pi^2 {\bf m}^2 / \beta ) }\over{{\bf m}^2}}&\cr&\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij})&\cr &\quad+ {{1}\over{2\pi V}} \sum_{{\bf m} \ne {\bf 0}}\ \sum_{(i,j) \not\in C\times C}q_iq_j \ {{ \exp(-\pi^2 {\bf m}^2 / \mu_{ij})}\over{{\bf m}^2}} &\cr &\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij}) &\cr&\quad- {{\pi}\over{2 V}} \sum_{(i,j)\in C\times C} q_i q_j\left ({{1}\over{\beta}} - {{1}\over{\theta_i}} - {{1}\over{\theta_j}}\right) &\cr&\quad- \sum_{i\in C} q_i^2 \left({{\beta}\over{\pi}}\right)^{1/2} - \sum_{i \not\in C}q_i^2 \left({{\theta_i}\over{2\pi}}\right)^{1/2} + {{2 \pi {\bf D}^2}\over{3 V}} &\cr &\quad+ {\varepsilon}(K).&(} ]

Note that in the limit that [\theta_i \to \infty] for all i, all of the Gaussians become compact and this formula reduces to equation ([link] for a large spherical crystal of point charges. The volume-dependent term [({{\pi}/{2 V}} )\textstyle\sum_{(i,j)\in C\times C} q_i q_j ({{1}/{\beta}} - {{1}/{\theta_i}} - {{1}/{\theta_j}}) ] is somewhat surprising, and would not appear in a more straightforward Ewald derivation such as our first Ewald derivation. It vanishes by unit-cell neutrality for the case of point charges or multipoles, and would also vanish if we chose a single [\alpha] for all Gaussian pairs. Higher angular momentum: Hermite Gaussians

| top | pdf |

More complex crystal charge distributions can be realized by considering normalized Hermite Gaussian distributions in place of the normalized spherical Gaussian distributions discussed above. Let [\rho] be a normalized spherical Gaussian charge distribution with exponent [\theta], centred at [{\bf R} = ({\bf R}_x,{\bf R}_y,{\bf R}_z)], i.e. [\rho({\bf r}) = (\theta/\pi)^{3/2} \exp[-\theta({\bf r} - {\bf R})^2] ]. The associated normalized Hermite Gaussians [\Lambda_{tuv}], where t, u, v are non-negative integers, are defined by[\eqalign{\Lambda_{tuv}({\bf r},\theta,{\bf R}) &= \left({{\theta}\over{\pi}}\right)^{3/2} \left(\partial / \partial{\bf R}_x\right)^t \left(\partial / \partial{\bf R}_y\right)^u \left(\partial / \partial{\bf R}_z\right)^v \cr &\quad\times \exp[-\theta({\bf r} - {\bf R})^2].} ]Note that this definition of Hermite Gaussians differs from that found in some textbooks, e.g. Helgaker et al. (2000[link]), in that we have normalized the Gaussian [\rho] it is derived from in order to allow exponents to tend to infinity to arrive at ideal point charges or point multipoles. We will refer to these Hermite Gaussians as normalized Hermite Gaussians to emphasize this distinction. The Hermite definition reduces to the spherical Gaussian when t, u and v are all zero, i.e. [\Lambda_{000}({\bf r},\theta,{\bf R}) = \rho({\bf r})].

Above we discussed the interaction between spherical Gaussian charge distributions [\rho]. Here we generalize these to higher angular momentum charge distributions [\varrho] given by[\varrho({\bf r},\theta,{\bf R},{\bf c}) = \textstyle\sum\limits_{tuv} {\bf c}_{tuv}\Lambda_{tuv}({\bf r},\theta,{\bf R}), ]where c is a vector of coefficients, a finite number of which are nonzero. For example, if [{\bf c}_{000}] is the only nonzero coefficient, we are reduced to the spherical charge distributions discussed above.

The fact that the partial derivatives within [\Lambda_{tuv}] are with respect to the position of the centre of the distribution allows us to easily calculate the Coulomb interaction energy between two normalized Hermite Gaussian distributions. The partial derivatives can be pulled out of the double integral, so that the interaction energy is given by the corresponding partial derivatives of the interaction energy of the two spherical Gaussian charge distributions:[\eqalign{E_{12} &= \int\limits_{\Re^3} \int\limits_{\Re^3} {{\Lambda_{t_1u_1v_1}({\bf r}_1,\theta_1,{\bf R}_1) \Lambda_{t_2 u_2 v_2}({\bf r}_2,\theta_2,{\bf R}_2)}\over{|{\bf r}_1 - {\bf r}_2|}} \, {\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 \cr &= \left({{\partial }\over{\partial{\bf R}_{1x}}}\right)^{t_1} \left({{\partial}\over{ \partial{\bf R}_{1y}}}\right)^{u_1} \left({{\partial}\over{ \partial{\bf R}_{1z}}}\right)^{v_1} \left({{\partial }\over{ \partial{\bf R}_{2x}}}\right)^{t_2} \cr &\quad\times \left({{\partial}\over{ \partial{\bf R}_{2y}}}\right)^{u_2} \left({{\partial }\over{\partial{\bf R}_{2z}}}\right)^{v_2} \left({{\theta_1\theta_2}\over{\pi^2}}\right)^{3/2} \int\limits_{\Re^3} \int\limits_{\Re^3} {{1}\over{|{\bf r}_1 - {\bf r}_2|}} \cr &\quad\times \exp[-\theta_1({\bf r}_1 - {\bf R}_1)^2] \exp[-\theta_2({\bf r}_2 - {\bf R}_2)^2]\,{\rm d}^3{\bf r}_1\,{\rm d}^3{\bf r}_2.} ]

Next we consider the interaction energy between a generalized charge distribution [\varrho_1] centred at [{\bf R}_1 \in U] and another generalized charge distribution [\varrho_2] centred at [{\bf R}_2 \in U], together with all the latter's periodic images, centred at [{\bf R}_2 + {\bf n} ], [{\bf n} \in \Omega(S,K)]. Examining equation ([link], we note that the terms to be differentiated all depend on [{\bf R}_{12} = {\bf R}_1 - {\bf R}_2 ]. Thus[\eqalignno{E_{12} &= \sum_{\bf n} \sum_{t_1 u_1 v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} {\bf c}_{1,t_1u_1v_1} {\bf c}_{2,t_2u_2v_2} &\cr &\times \left({{\partial }\over{\partial{\bf R}_{12x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R}_{12y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R}_{12z}}} \right)^{v_1+v_2} &\cr&\quad\times \left \lbrace {{{\rm erfc}(\mu_{12\alpha}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over {|{\bf R}_{12} - {\bf n}|}} - {{{\rm erfc}(\mu_{12}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over {|{\bf R}_{12} - {\bf n}|}} \right \rbrace &\cr &\quad+ {{1}\over{\pi V}} \sum_{{\bf m} \ne {\bf 0}} {{ \exp(-\pi^2 {\bf m}^2 / \mu_{12\alpha}) }\over{{\bf m}^2}} C_1({\bf m}) C_2(-{\bf m})&\cr &\quad- {\bf c}_{1,000} {\bf c}_{2,000}{{\pi}\over{\mu_{p12\alpha} V}} + {{1}\over{\pi}} H_{S,K}^{12}({\bf R}_{12}) &\cr &\quad+ {\varepsilon}(K),&(} ]where [C_i({\bf m})], i =1, 2 are given by[\eqalign{C_i({\bf m}) &= \sum_{t_i u_i v_i}{\bf c}_{i,t_iu_iv_i} \cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{ix}}}\right)^{t_i} \left({{\partial}\over{ \partial{\bf R}_{iy}}}\right)^{u_i} \left({{\partial}\over{ \partial{\bf R}_{iz}}} \right)^{v_i} \exp(-2 \pi i {\bf m}\cdot{\bf R}_i) } ]and[\eqalign{H_{S,K}^{12}({\bf R}_{12}) &= \sum_{t_1 u_1 v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} {\bf c}_{1,t_1u_1v_1} {\bf c}_{2,t_2u_2v_2} \cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{12x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R}_{12y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R}_{12z}}} \right)^{v_1+v_2} \cr &\quad\times H_{S,K}({\bf R}_{12}).} ]Note that by the arguments leading to equation ([link], [H_{S,K}({\bf r}) ] is quadratic in r, i.e. derivatives of higher than second order vanish.

We now can discuss the Coulomb energy of the central unit cell U within a large spherical crystal when the charge density in the unit cell is described by a basis [\Theta] of normalized Hermite Gaussians and a set of expansion points [{\bf R}_1,\ldots,{\bf R}_N] within U that are repeated periodically in the crystal (for example, some of the expansion points could be atomic nuclei). Let [\theta_1,\ldots,\theta_L ] be positive Gaussian exponents (or infinite, in the case of nuclear charge or point charges or ideal multipoles) and let the charge density about the expansion point [{\bf R}_i] be given by[\varrho_i({\bf r},{\bf R}_i,\Theta) = \textstyle\sum\limits_{l = 1}^L \textstyle\sum\limits_{tuv} {\bf c}_{i,l,tuv} \Lambda_{tuv}({\bf r},\theta_l,{\bf R}_i), ]where, as above, for each i, l the coefficients [{\bf c}_{i,l,tuv}] are nonzero for a finite number of tuv. As above we can separate the exponents [\theta_l] into compact and diffuse, i.e. [l \in C] or [l \in D], according to whether [\theta_l \ge 2 \beta ] or [\theta \,\lt\, 2 \beta], respectively. Let [q_i = \textstyle\sum_{l = 1}^L {\bf c}_{i,l,000} ] denote the net charge at the expansion point [{\bf R}_i].

We are interested in the Coulomb energy of charge densities [\varrho_i({\bf r},{\bf R}_i,\Theta) ], [i = 1,\ldots,N], interacting with each other and their periodic images at expansion points [{\bf R}_i + {\bf n}], [{\bf n} \in \Omega(S,K) ] for large K. As above, the direct Coulomb interaction of [\varrho_i({\bf r},{\bf R}_i,\Theta)] with itself is not permitted. This can be modified for more sophisticated treatments e.g. in the context of the Coulomb integrals in a periodic density-functional-theory approach, but clearly a nuclear charge cannot interact with itself. The unit cell is assumed to be neutral, i.e. [\textstyle\sum_{i=1}^N q_i = 0]. Then, using equations ([link] and ([link] we can write the Coulomb energy of [\underline\varrho^{\lbrace\bf N\rbrace}] = [\lbrace \varrho_1,\ldots,\varrho_N\rbrace] within the spherical crystal as[\eqalignno{ E_{S,K}(\underline\varrho^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'} \sum_{i=1}^N \sum_{l_i \in C} \sum_{t_i u_i v_i} {\bf c}_{i,l_i,t_i u_i v_i} &\cr &\quad\times \sum_{j=1}^N \sum_{l_j \in C} \sum_{t_j u_j v_j} (-1)^{(t_j+u_j+v_j)} {\bf c}_{j,l_j,t_j u_j v_j} &\cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{ijx}}}\right)^{t_i+t_j} \left({{\partial}\over{ \partial{\bf R}_{ijy}}}\right)^{u_i+u_j} \left({{\partial}\over{ \partial{\bf R}_{ijz}}} \right)^{v_i+v_j} &\cr &\quad\times \left \lbrace {{{\rm erfc}(\beta^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over {|{\bf R}_{ij} - {\bf n}|}} - {{{\rm erfc}(\mu_{l_i l_j}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over {|{\bf R}_{ij} - {\bf n}|}} \right\rbrace&\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}}{{1}\over{{\bf m}^2}} \exp(-\pi^2 {\bf m}^2 /2 \beta) \sum_{l_1\in C}S_{l_1}({\bf m}) &\cr &\quad\times \exp(-\pi^2 {\bf m}^2 /2 \beta) \sum_{l_2\in C} S_{l_2}(-{\bf m})&\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} {{1}\over{{\bf m}^2}} \sum_{(l_1,l_2)\not\in C\times C} \exp(-\pi^2 {\bf m}^2 / \theta_{l_1}) &\cr &\quad\times \exp(-\pi^2 {\bf m}^2 / \theta_{l_2}) S_{l_1}({\bf m}) S_{l_2}(-{\bf m}) &\cr &\quad- {{\pi}\over{2 V}} \sum_{l_1 \in C} \sum_{l_2 \in C} \sum_{i=1}^N \sum_{j=1}^N {\bf c}_{i,l_1,000} {\bf c}_{j,l_2,000} &\cr &\quad\times \left({{1}\over{\beta}} - {{1}\over{\theta_{l_1}}} - {{1}\over{\theta_{l_2}}}\right) &\cr&\quad- \sum_{i=1}^N E_{\rm self}(\varrho_i) + {{2 \pi {\bf D}^2}\over{3 V}} + {\varepsilon}(K), &(} ]where the structure factors [S_{l}({\bf m})] are given by[\eqalignno{S_{l}({\bf m}) &= \sum_{i=1}^N \sum_{tuv} {\bf c}_{i,l,tuv} \left({{\partial }\over{\partial{\bf R}_{ix}}}\right)^t \left({{\partial }\over{\partial{\bf R}_{iy}}}\right)^u \left({{\partial }\over{\partial{\bf R}_{iz}}}\right)^v &\cr &\quad\times\exp(2 \pi i {\bf m}\cdot{\bf R}_i).&(} ][E_{\rm self}(\varrho_i)] is given by[\eqalign{E_{\rm self}(\varrho_i) &= \lim_{{\bf R} \to {\bf 0}} \sum_{t_1u_1v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} \cr &\quad\times \left({{\partial }\over{\partial{\bf R_x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R_y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R_z}}} \right)^{v_1+v_2} \cr &\quad\times \left\lbrace \sum_{(l_1,l_2) \in C\times C}{\bf c}_{i,l_1,t_1u_1v_1}{\bf c}_{i,l_2,t_2u_2v_2} {{{\rm erf}(\beta^{1/2} |{\bf R}|)}\over{|{\bf R}|}} \right. \cr &\quad+ \left. \sum_{(l_1,l_2) \not \in C\times C}{\bf c}_{i,l_1,t_1u_1v_1}{\bf c}_{i,l_2,t_2u_2v_2} {{{\rm erf}(\mu_{12}^{1/2} |{\bf R}|)}\over{|{\bf R}|}} \right\rbrace,} ]where [1 / \mu_{12} = 1 / \theta_{l_1} + 1 / \theta_{l_2}], and the unit-cell dipole components are given by[\eqalign{{\bf D}_x &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_x + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,100}, \cr {\bf D}_y &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_y + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,010}, \cr {\bf D}_z &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_z + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,001}.} ]

In the limit that all Gaussians are compact with [\theta_l \to \infty ], [E_{S,K}(\underline\varrho^{\lbrace\bf N\rbrace})] reduces to the Coulomb energy of ideal multipoles positioned at the expansion points.

3.5.4. Computational efficiency

| top | pdf |

The Ewald-sum methods discussed above lead to calculations that are much faster than straightforward summation, as well as clarifying the nature of the conditional convergence. However, there are opportunities for further significant speed-ups, which we cover in this section. The computational cost of the Ewald sum is dominated by the direct and reciprocal sum. We cover these in turn. The direct sum

| top | pdf |

For all the cases discussed above – point charge, normalized spherical Gaussians and the more complex densities made up of sums over a basis of normalized Hermite Gaussians – the Ewald direct sum involves terms like [|{\bf R}_{12} - {\bf n}| ] and in the Hermite case partial derivatives with respect to the x, y and z components of these terms. Given the efficiency opportunities in the reciprocal sum to be discussed below, it is advantageous to choose an Ewald parameter (referred to above as α in the point-charge case or β for the compact–compact Gaussian or Hermite Gaussian interactions) so that direct-sum terms can be neglected past a fixed cutoff, i.e. when [|{\bf R}_{12} - {\bf n}| > R_C]. Given a fixed Ewald parameter α, the cutoff may be fixed independently of the unit-cell size. At constant or near constant material density (number of atoms per cubic ångstrom) this leads to `linear scaling', i.e. the computational cost of the Ewald direct sum grows linearly with the unit-cell volume or number of atoms in the unit cell U.

For the point-charge or normalized spherical Gaussian case, the remaining issue is efficient evaluation of the erfc function. Various table lookup schemes have been advanced, as well as a polynomial expression (Allen & Tildesley, 1987[link]). Ewald summation for ideal multipoles up to quadrupole have been discussed by Smith (1982[link], 1998[link]) and Aguado & Madden (2003[link]), using optimized combinations of terms to minimize the cost of the direct sum. A more general scheme, covering higher-order multipoles, was proposed by Sagui et al. (2004[link]). In the latter scheme the McMurchie–Davidson recursion (McMurchie & Davidson, 1978[link]; Helgaker et al., 2000[link]) was used to calculate the higher-order derivatives of [{\rm erfc}(\alpha^{1/2} r) / r ] for the multipole–multipole interaction tensor. This same approach was used by Cisneros et al. (2006[link]) to calculate the interaction tensor for normalized Hermite Gaussian interactions. Although originally formulated for Hermite Gaussian interactions, it applies to more general partial derivative calculations. Let r = (x, y, z) be a point in [\Re^3] with [|{\bf r}| = r], and let g(r) be a smooth function of r. For example, [g(r) = {\rm erf}(\alpha^{1/2} r)/r ]. We are interested in the partial derivatives with respect to x, y and z of g(r). Let R(0, 0, 0, 0) = g(r), and for n > 0 let R(0, 0, 0, n) = (1/r)(d/dr)R(0, 0, 0, n − 1). For non-negative integers t, u, v, let [R(t,u,v,n)] = [(\partial/\partial x)^t (\partial/\partial y)^u (\partial/\partial z)^v R(0,0,0,n) ]. Then we have the following recursion:[\eqalign{R(t+1,u,v,n) &= x R(t,u,v,n+1) + t R(t-1,u,v,n+1), \cr R(t,u+1,v,n) &= y R(t,u,v,n+1) + u R(t,u-1,v,n+1), \cr R(t,u,v+1,n) &= z R(t,u,v,n+1) + v R(t,u,v-1,n+1).} ]

The proof of the last recursion depends on the fact that [\partial / \partial z = (z/r) \partial / \partial r ] and on Leibniz' generalization of the product rule for differentiation, namely [(fg)^{(k)} = \textstyle\sum_{j=0}^k {k \choose j} f^{(j)}g^{(k-j)}], and is as follows:[\eqalign{R(t,u,v+1,n) &= \left(\partial / \partial x\right)^t \left(\partial / \partial y\right)^u \left(\partial / \partial z\right)^v \left(\partial / \partial z\right) R(0,0,0,n) \cr &= \left(\partial / \partial x\right)^t \left(\partial / \partial y\right)^u \left(\partial / \partial z\right)^v z (1/r) \cr &\quad\times ({\rm d}/{\rm d}r) R(0,0,0,n) \cr &= \left(\partial / \partial x\right)^t \left(\partial / \partial y\right)^u \left(\partial / \partial z\right)^v z R(0,0,0,n+1) \cr &= \left(\partial / \partial x\right)^t \left(\partial / \partial y\right)^u [z \left(\partial / \partial z\right)^v R(0,0,0,n+1) \cr &\quad + v \left(\partial / \partial z\right)^{v-1} R(0,0,0,n+1) ]\cr &= z R(t,u,v,n+1) + v R(t,u,v-1,n+1).} ]The other two recursions are proved similarly. To apply this to the case [g(r) = {\rm erf}(\alpha^{1/2} r)/r] we recall from equation ([link] that [(2 /\pi^{1/2}) F_0(x^2) = {\rm erf}(|x|) / |x|] and that for the higher-order Boy's functions [F_n(x)], [({\rm d}/{\rm d}x )F_n(x) = -F_{n+1}(x) ]. Efficient evaluation of the Boy's functions is discussed by Helgaker et al. (2000[link]).

Further economies for point multipoles can be achieved by using symmetry. For example, it is only necessary to keep track of six quadrupolar terms, with a similar reduction for Hermite Gaussians. A general discussion, including rotation of multipoles or Hermite Gaussians from local molecule frame to laboratory frame, is provided in Sagui et al. (2004[link]) as well as in Cisneros et al. (2006[link]). Reciprocal-sum speed-ups: the smooth particle mesh Ewald and fast Fourier Poisson methods

| top | pdf |

The smooth particle mesh Ewald (PME) (Essmann et al., 1995[link]) and fast Fourier Poisson (FFP) (York & Yang, 1994[link]) methods both rely on the three-dimensional fast Fourier transform (3DFFT) for computational efficiency. In combination with the above use of cutoffs in the direct sum, they lower the computational cost of the Ewald sum from [O (N^{3/2})] (Perram et al., 1988[link]) to [O [N \log(N)]], where N is the number of particles and e.g. [O (N^{3/2})] denotes a quantity such that [O (N^{3/2})/N^{3/2}] remains bounded as [N \to \infty]. For the point-charge case, the key insight is that if the point charges happen to be located on a regular grid, the structure factors [equation ([link]] become discrete Fourier transforms and can be rapidly calculated using FFTs.

Like other variants of the original particle–particle particle–mesh (P3M) method (Hockney & Eastwood, 1981[link]), the smooth PME method interpolates charge onto grid positions. The smooth PME, like the original PME algorithm (Darden et al., 1993[link]), formulates this problem as one of interpolation of the complex exponentials appearing in the structure factor. This interpolation process is very accurate at low frequency, i.e. for small [|{\bf m}|], and becomes progressively less accurate as [|{\bf m}|] grows. However, the weight function [\exp(-\pi^2 {\bf m}^2 / \alpha) / {\bf m}^2] multiplying the structure factor in the reciprocal sum severely damps the high-frequency (large [|{\bf m}|]) terms, thus preserving overall accuracy. In contrast, the FFP method relies on the fact that [\exp(-\pi^2 {\bf m}^2 / \alpha) / {\bf m}^2 S({\bf m}) ] is the Fourier transform of a superposition of normalized spherical Gaussian charge distributions [q_j \rho_j], centred at [{\bf r}_j ], as discussed above. This Fourier transform is approximated by the FFT of the gridded charge density obtained by direct sampling of the Gaussians onto the grid (similar to the use of the FFT in macromolecular X-ray crystallography). Note that the smooth PME method has been easily extended to the treatment of London dispersion interactions, since these involve similar structure factors with a different weight function [g_p, p=6,8,10]. The extension is not so straightforward with the FFP approach. However, as discussed below, the FFP approach has advantages in the case of continuous charge density when more than one Gaussian exponent is involved.

For the case of point charges, P3M methods based on least-squares approximation of the exact Ewald reciprocal sum gradient (Pollock & Glosli, 1996[link]; Darden et al., 1997[link]) can be made more efficient than the smooth PME method at the same accuracy, despite requiring more FFTs (Deserno & Holm, 1998[link]). However, for the case of ideal multipoles or normalized Hermite Gaussians, the number of such FFTs becomes prohibitive, whereas the smooth PME, which relies on analytic differentiation of B-splines, does not need more FFTs as the order of differentiation (t + u + v) in equation ([link] grows. The smooth PME approach

| top | pdf |

We begin by re-expressing the complex exponentials appearing in the structure factors in equations ([link], ([link] and ([link]. Recall the form of m from equation ([link], and let [M_1,M_2,M_3] be positive integers of the scale of the unit-cell dimensions. Then for [{\bf r} \in U ] we can write[\eqalign{{\bf m}\cdot{\bf r} &= m_1 {\bf a}_1^{*} \cdot {\bf r} + m_2 {\bf a}_2^{*} \cdot {\bf r} + m_3 {\bf a}_3^{*} \cdot {\bf r} \cr &= z_1 w_1 + z_2 w_2 + z_3 w_3,} ]where [w_\alpha = M_\alpha {\bf a}_\alpha^{*} \cdot {\bf r}] and [z_\alpha = m_\alpha / M_\alpha], [\alpha = 1,2,3], and thus[\exp(2 \pi i {\bf m} \cdot {\bf r}) = \exp(2 \pi i z_1 w_1)\exp(2 \pi i z_2 w_2) \exp(2 \pi i z_3 w_3) . ]

To motivate the smooth PME approach, suppose we can approximate [\exp(2 \pi i z w) ] by some function [g(w;z)] of the form[g(w;z) = \lambda(z)\textstyle \sum\limits_{k=-\infty} ^ \infty B(w-k) \exp(2 \pi i z k),\eqno( ]where [B(w)] is a well behaved function with continuous derivatives as needed for the general structure factor in equation ([link]. We first examine the structure factor for point-charge Ewald sums, defined in equation ([link]. This can be approximated as[\eqalignno{S({\bf m}) &= \textstyle\sum\limits_{j=1}^N q_j \exp(2 \pi i z_1 w_{1j})\exp(2 \pi i z_2 w_{2j})\exp(2 \pi i z_3 w_{3j}) & \cr &\simeq \lambda(z_1)\lambda(z_2) \lambda(z_3) & \cr &\quad\times \textstyle\sum\limits_{k_1=-\infty} ^ \infty \textstyle\sum\limits_{k_2=-\infty} ^ \infty \textstyle\sum\limits_{k_3=-\infty} ^ \infty \textstyle\sum\limits_{j=1}^N q_j B(w_{1j} - k_1) B(w_{2j} - k_2) & \cr &\quad\times B(w_{3j} - k_3) \exp[2\pi i (z_1 k_1+z_2 k_2 + z_3 k_3)]& \cr &= \lambda(z_1)\lambda(z_2) \lambda(z_3) & \cr &\quad\times \textstyle\sum\limits_{k_1 = -{({M_1}/{2})} + 1}^{{{M_1}/{2}}}\ \textstyle\sum\limits_{k_2 = -{({M_2}/{2})}+1}^{{{M_2}/{2}}}\ \textstyle\sum\limits_{k_3 = -{({M_3}/{2})}+1}^{{{M_3}/{2}}} Q(k_1,k_2,k_3) & \cr &\quad\times \exp\left[2 \pi i \left({{m_1 k_1}\over{M_1}} + {{m_2 k_2}\over{M_2}} + {{m_3 k_3}\over{M_3}}\right)\right],&(} ]where [Q(k_1,k_2,k_3)] is given by[\eqalignno{Q(k_1,k_2,k_3) &= \textstyle\sum\limits_{j=1}^N q_j \textstyle\sum\limits_{l_1=-\infty}^\infty B(w_{1j} - k_1 - l_1 M_1)&\cr &\quad\times \textstyle\sum\limits_{l_2=-\infty}^\infty B(w_{2j} - k_2 - l_2 M_2) & \cr &\quad\times \textstyle\sum\limits_{l_3=-\infty}^\infty B(w_{3j} - k_3 - l_3 M_3) &(} ]and we have used the fact that [\exp[2 \pi i m_\alpha (k_\alpha + M_\alpha)/M_\alpha]] = [\exp(2 \pi i m_\alpha k_\alpha/M_\alpha)] for [\alpha = 1,2,3].

Next we turn to the general structure factors [S_l({\bf m})] from equation ([link]. Define [w_{\alpha j}], [\alpha=1,2,3] as above for expansion points [{\bf R}_j \in U], and note that[\eqalign{{{\partial}\over{\partial {\bf R}_{jx}}} &= {{\partial}\over{\partial w_{1j}}}\ {{\partial w_{1j}}\over{\partial {\bf R}_{jx}}} + {{\partial}\over{\partial w_{2j}}}\ {{\partial w_{2j}}\over{\partial {\bf R}_{jx}}} + {{\partial}\over{\partial w_{3j}}}\ {{\partial w_{3j}}\over{\partial {\bf R}_{jx}}} \cr &= b_1\ {{\partial}\over{\partial w_{1j}}} + b_2\ {{\partial}\over{\partial w_{2j}}} + b_3\ {{\partial}\over{\partial w_{3j}}},} ]where [b_1, b_2, b_3] are constants depending on the reciprocal bases [{\bf a}_\alpha^{*}] and on [M_\alpha], [\alpha = 1,2,3 ]. Similar expansions hold for the other partial derivatives. Thus we can approximate the general structure factor [S_l({\bf m})] by[\eqalignno{ S_l({\bf m}) &= \sum_{j=1}^N \sum_{tuv} {\bf c}_{j,l,tuv} \left({{\partial }\over{\partial{\bf R}_{jx}}}\right)^t \left({{\partial }\over{\partial{\bf R}_{jy}}}\right)^u \left({{\partial }\over{\partial{\bf R}_{jz}}}\right)^v & \cr &\quad\times \exp(2 \pi i {\bf m}\cdot{\bf R}_j) & \cr &= \sum_{j=1}^N \sum_{tuv} {\bf d}_{j,l,tuv} \left({{\partial }\over{\partial w_{1j}}}\right)^t \left({{\partial }\over{\partial w_{2j}}}\right)^u \left({{\partial }\over{\partial w_{3j}}}\right)^v & \cr &\quad\times \exp(2 \pi i z_1 w_{1j})\exp(2 \pi i z_2 w_{2j})\exp(2 \pi i z_3 w_{3j}) & \cr &\simeq \lambda(z_1)\lambda(z_2) \lambda(z_3) & \cr &\quad\times \sum_{k_1 = -{({M_1}/{2})} + 1}^{{{M_1}/{2}}} \ \sum_{k_2 = -{({M_2}/{2})}+1}^{{{M_2}/{2}}} \ \sum_{k_3 = -{({M_3}/{2})}+1}^{{{M_3}/{2}}} Q_l(k_1,k_2,k_3) & \cr &\quad\times \exp\left[2 \pi i \left({{m_1 k_1}\over{M_1}} + {{m_2 k_2}\over{M_2}} + {{m_3 k_3}\over{M_3}}\right)\right], &(} ]where [{\bf d}_{j,l,tuv}] are the transformed Hermite coefficients obtained from [{\bf c}_{j,l,tuv}] by change of variables, using [b_1,b_2,b_3] and where [Q_l(k_1,k_2,k_3)] is given by[\eqalign{Q_l(k_1,k_2,k_3) &= \sum_{j=1}^N \sum_{tuv} {\bf d}_{j,l,tuv} \cr &\quad\times \sum_{l_1=-\infty}^\infty \left({{\partial }\over{\partial w_{1j}}}\right)^t B(w_{1j} - k_1 - l_1 M_1) \cr &\quad\times \sum_{l_2=-\infty}^\infty \left({{\partial }\over{\partial w_{2j}}}\right)^t B(w_{2j} - k_2 - l_2 M_2) \cr &\quad\times \sum_{l_3=-\infty}^\infty \left({{\partial }\over{\partial w_{3j}}}\right)^t B(w_{3j} - k_3 - l_3 M_3) .} ]

In equation ([link] we recognize the approximation to [S({\bf m}) ] as [\lambda(z_1)\lambda(z_2) \lambda(z_3)] times the three-dimensional discrete Fourier transform of [Q(k_1,k_2,k_3)], which can be calculated rapidly by the 3DFFT. If the function B(w) has finite support then [Q(k_1,k_2,k_3)] can be calculated rapidly [in order(N) time] as well. Similarly we have the 3DFFT of [Q_l(k_1,k_2,k_3)] in equation ([link]. Thus the efficient calculation of structure factors depends on the approximation given in equation ([link], which in turn rests on the choice of B(w). The smooth PME relies on B-splines. Properties of B-splines

| top | pdf |

We proceed by deriving the facts needed about B-splines. The B-splines of order n, Bn(w) are functions of a real variable defined as follows:[\eqalign{B_1(w) &= \left\{\let\normalbaselines\relax\openup1pt\matrix{ 1\hfill & \hbox{ if }0 \le w \le 1\hfill \cr 0\hfill& \hbox{otherwise}\hfill} \right. \cr&\cr B_2(w) &= \left\{\let\normalbaselines\relax\openup1pt\matrix{w\hfill & \hbox{ if }0 \le w \le 1\hfill\cr 2-w\hfill & \hbox{ if }1 \le w \le 2\hfill \cr 0 \hfill & \hbox{ otherwise}\hfill}\right. \cr &\cr B_{n+1}(w) &= \textstyle\int\limits_{-\infty}^\infty B_n(w-v) B_1(v) \,{\rm d}v \cr &= \textstyle\int\limits_0^1 B_n(w-v) B_1(v) \,{\rm d}v.} ]

We use the following properties of the B-splines:

  • (1) for [n \ge 1] and all w, [0 \le B_n(w) ] and [B_n(w) = 0] for [w \,\lt\, 0] or [w > n];

  • (2) for [n \ge 1] and all w, [\textstyle\sum_{j = -\infty}^\infty B_n(w-j) = 1 ];

  • (3) for [n \ge 1, \ 0 \le w \le n], [B_n(n-w) = B_n(w) ];

  • (4) for [n \ge 2], [0 \le w \le n], [({\rm d}/{\rm d}w) B_{n}(w) = B_{n-1}(w)\ -] [B_{n-1}(w-1)];

  • (5) for [n \ge 2], [0 \le w \le n], [B_n(w) = [w B_{n-1}(w)] + [(n - w) B_{n-1}(w-1)] / (n-1) ].

The first property to be established is clearly true for n = 1 and if true for n, then in the above integral recursively defining Bn+1 if [w \,\lt\, 0] or [w > (n+1)] the integrand must be zero for [0 \le v \le 1], so [B_{n+1}(w) = 0 ]. Thus it is true in general by mathematical induction.

The second property is also clearly true for n = 1. Assume it is true for n. Then [\textstyle\sum_{j = -\infty}^\infty B_{n+1}(w-j) ] = [\textstyle\int_0^1 [\textstyle\sum_{j=-\infty}^\infty B_n(w-v -j)] B_1(v) \,{\rm d}v ] = [\textstyle\int_0^1 B_1(v) \,{\rm d}v = 1]. Thus it is true in general by mathematical induction.

Next we show symmetry: [B_n(n-w) = B_n(w), \ 0 \le w \le n]. Again this is clearly true for n = 1. Assuming it is true for n, we have[\eqalign{B_{n+1}(n+1 - w) &= \textstyle\int\limits_0^1 B_n(n+1 - w - v) B_1(v) \,{\rm d}v \cr &= \textstyle\int\limits_0^1 B_n[n+1 - w - (1-v')] B_1(1-v') \,{\rm d}v' \cr &= \textstyle\int\limits_0^1 B_n[n - (w-v')] B_1(1-v')\,{\rm d}v' \cr &= \textstyle\int\limits_0^1 B_n(w - v')B_1(v') \,{\rm d}v' \cr &= B_{n+1}(w), } ]where we substituted [v = (1-v')] in the second line above.

The fourth property is an efficient recursion for the derivatives of B-splines. It can be seen directly for n = 2 for [w \ne 0,1, \hbox{or} \ 2 ]. Assume it is true for [n \ge 2]. Then for [0 \le w \le n+1 ],[\eqalign{{{\rm d}\over{{\rm d}w}} B_{n+1}(w) &= \int_0^1 {{\rm d}\over{{\rm d}w}} B_n(w-v) B_1(v) \,{\rm d}v \cr &= \int_0^1 B_{n-1}(w-v) B_1(v)\,{\rm d}v \cr &\quad- \int_0^1 B_{n-1}(w-1-v) B_1(v) \,{\rm d}v \cr &= B_n(w) - B_n(w-1),} ]thus it is true in general by induction.

The final property is an efficient recursion for [B_n(w)]. It can be seen directly for n = 2. Assume it is true for [n \ge 2]. Then for [0 \le w \le n+1],[\eqalign{B_{n+1}(w) &= \int_0^1 \left[{{w-v}\over{n-1}} B_{n-1}(w-v) \right.\cr &\quad + \left. {{n-(w-v)}\over{n-1}} B_{n-1}(w-v-1)\right] B_1(v)\,{\rm d}v \cr &= {{n}\over{n-1}} B_n(w-1) \cr &\quad + {{1}\over{n-1}} \int_0^1 (w-v) {{\rm d}\over{{\rm d}w}} B_n(w-v) B_1(v)\,{\rm d}v \cr &= {{n}\over{n-1}} B_n(w-1) \cr &\quad - {{1}\over{n-1}} \int_0^1 (w-v) {{\rm d}\over{{\rm d}v}} B_n(w-v) B_1(v)\,{\rm d}v \cr &= {{n}\over{n-1}} B_n(w-1) \cr &\quad- {{1}\over{n-1}}[(w-1)B_n(w-1) - w B_n(w)] \cr &\quad + {{1}\over{n-1}} \int_0^1 B_n(w-v) {{\rm d}\over{{\rm d}v}}(w-v)B_1(v) \,{\rm d}v \cr &= {{n+1 - w}\over{n-1}} B_n(w-1) + {{w}\over{n-1}} B_n(w) \cr &\quad- {{1}\over{n-1}} B_{n+1}(w),} ]where we have used property (4[link]) in the third line and integration by parts in the fifth and sixth lines. The induction step follows by moving the last term on the right-hand side over to the left side of the equation. B-spline approximation to trigonometric functions

| top | pdf |

Next we examine properties of the B-spline approximations to the complex exponentials that appear in the structure factors. Let [g_n(w;z)] be given by[g_n(w;z) = \lambda_n(z) \textstyle\sum\limits_{j=-\infty}^\infty B_n(w - j) \exp(2 \pi i z j).\eqno( ]The complex number [\lambda_n(z)] is chosen to optimize the approximation of [\exp(2 \pi i z w)] by [g_n(w;z)] as shown below. It is convenient to first study the properties of the function [h_n(w;z)] given by[\eqalign{h_n(w\semi z) &= \exp(-2 \pi i z w) g_n(u\semi z) \cr &= \lambda_n(z)\textstyle\sum\limits_{j=-\infty}^\infty B_n(w - j) \exp[-2 \pi i z (w-j)],} ]which is an approximation to the constant one. From its definition it is clear that [h_n(w;z)] is periodic in w with period one, i.e. [h_n(w+1;z) = h_n(w;z)]. Therefore it can be expanded in a Fourier series:[h_n(w;z) = \textstyle\sum\limits_{k=-\infty}^\infty c_n(k;z) \exp(2 \pi i k w).\eqno( ]The Fourier coefficients are given by[\eqalign{c_n(k\semi z) &= \textstyle\int\limits_0^1 h_n(w\semi z) \exp(-2 \pi i k w) \,{\rm d}w \cr &= \lambda_n(z)\textstyle\int\limits_0^1 \textstyle\sum\limits_{j=-\infty}^\infty B_n(w - j) \cr &\quad\times \exp[-2 \pi i z (w-j)] \exp(-2 \pi i k w) \,{\rm d}u\cr &= \lambda_n(z)\textstyle\sum\limits_{j=-\infty}^\infty \textstyle\int\limits_0^1 B_n(w - j) \cr &\quad\times \exp[-2 \pi i (z+k) (w-j)]\, {\rm d}u\cr &= \lambda_n(z)\textstyle\int\limits_{-\infty}^\infty B_n(v) \exp[-2 \pi i (z+k) v]\, {\rm d}v \cr &= \lambda_n(z) \ \hat{B}_n (z+k),} ]where [\hat{B}_n] denotes the Fourier transform of [B_n]. Since [B_n] are defined recursively by convolutions, its Fourier transform is obtained simply from that of [B_1]: [\hat{B}_n = (\hat{B}_1)^n ]. Here we are using the fact that if f, g and h are related by convolution, [h(w) = (f\star g)(w) = \textstyle\int_{-\infty}^\infty f(w-v) g(v)\,{\rm d}v ], then [\hat{h} = \hat{f}\hat{g}]. Since [\hat{B}_1] can be obtained by straightforward integration, we have[c_n(k;z) = \left\{\let\normalbaselines\relax\openup2pt\matrix{ \lambda_n(z)\hfill& \hbox{if }z + k = 0\hfill \cr \lambda_n(z) {\displaystyle\left[{{1 -\exp(-2 \pi i z)}\over{2 \pi i}}\right]^n \left({{1}\over{z+k}}\right)^n} \hfill& \hbox{otherwise.}\hfill}\right. ]

The original smooth PME method set [\lambda_n(z) = \bar{\lambda}_n(z) ], where [\bar{\lambda}_n(z)] was chosen to enforce interpolation: [g_n(w;z) = \exp(2\pi i zw)] or [h_n(w;z) = 1] for integer values of w. This was related to the so-called Euler spline, studied for more general complex-valued interpolation problems by Schoenberg (1973[link]). The special case of interpolating the function [\exp(2\pi i zw) ] is simpler, due to the Fourier-series representation, allowing us to to derive necessary properties in a self-contained way. From equations ([link] and the above expression for [c_n(k;z)], noting that [c_n(k;z) = 0] for integral z with [z+k \ne 0], the interpolation condition can be written[1 = \left\{\let\normalbasesline\relax\openup3pt\matrix{\bar{\lambda}_n(z)\hfill & z \hbox{ an integer} \hfill\cr \bar{\lambda}_n(z) \displaystyle{\left[{{1 - \exp(-2 \pi i z)}\over{2 \pi i}}\right]^n\sum_{k=-\infty}^\infty \left({{1}\over{z+k}}\right)^n}\hfill & \hbox{otherwise.} \hfill}\right. ]

Note that except for [\bar{\lambda}_n(z)], the expressions appearing in the interpolation condition are periodic in z with period one. The interpolation condition thus forces [\bar{\lambda}_n(z)] to be periodic as well, and so we focus attention on the range [-\textstyle{1\over2} \le z \le {1\over 2}]. If [n \ge 2] the sum in the above interpolation condition is bounded, so in order to be able to choose [\bar{\lambda}_n(z)] to satisfy the interpolation condition it is sufficient to demonstrate that [\textstyle\sum_{k=-\infty}^\infty [{{1}/({z+k})}]^n \ne 0 ] for [-\textstyle{1\over2} \le z \le \textstyle{1\over2}]. If n is even, the terms in this sum are all positive so it is always possible to interpolate. If n is odd and z = ½ then the sum is zero, since [(1/z)^n + [1 / (-1 +z)]^n = 0 ], [[1/(1+z)^n + 1 / (-2+z)]^n = 0] etc. Similarly the interpolation condition fails for z = −½. However, for [0 \,\lt\, z \,\lt\, \textstyle{1\over2}], [(1/z)^n + [1 / (-1 +z)]^n > 0], [[1/(1+z)^n + 1 / (-2+z)]^n > 0] etc. Similarly for [-\textstyle{1\over2} \,\lt\, z \,\lt\, 0 ]. Thus we can say that for [n \ge 2] and [-\textstyle{1\over2} \,\lt\, z \,\lt\, \textstyle{1\over2} ] the interpolation condition holds, and in this case we can write [h_n] as[h_n(w;z) = \textstyle\sum\limits_{k = -\infty}^\infty \gamma_n(k;z) \exp(2 \pi i k w),\eqno( ]where[\gamma_n(k;z) = \left({{1}\over{z + k}}\right)^n \bigg/ \sum_{j= -\infty}^\infty \left({{1}\over{z + j}}\right)^n.\eqno( ]The interpolating [g_n(w;z)] is obtained from this [h_n] by multiplying by [\exp(2 \pi i z w)]. Numerically, it is more convenient to use the alternative expression from equation ([link], with [\bar{\lambda}_n(z) ] given by[\bar{\lambda}_n(z) = \left[\textstyle\sum\limits_{j=0}^n B_n(j) \exp(-2 \pi i z j)\right]^{-1}. ]Note that due to the above results for [h_n], [\textstyle\sum_{j=0}^n B_n(j) \exp(-2 \pi i z j) \ne 0 ] if [n \ge 2] and [-\textstyle{1\over2} \,\lt\, z \,\lt\, \textstyle{1\over2}].

Subsequent to the work by Essmann et al. (1995[link]), we found (Darden et al., 1997[link]) that a more accurate approximation was given by least-squares fitting of [\exp(2 \pi i z w) ] by [g_n(w;z)], or equivalently 1 by [h_n(w;z)]. Thus, with [h_n] given by equation ([link] we are looking for a complex number [\tilde{\lambda} = \tilde{\lambda}_n(z)] that minimizes the mean-square error MSE, given by[\eqalign{{\rm MSE} &= \textstyle\int\limits_0^1 [1 - \tilde{\lambda} h_n(w\semi z) ][1 - \tilde{\lambda} h_n(w\semi z) ]^{*} \,{\rm d}w \cr &= 1 - 2 Re(\tilde{\lambda}) \gamma_n(0\semi z)+ [Re(\tilde{\lambda})^2 + Im(\tilde{\lambda})^2]\textstyle\sum\limits_{k=-\infty}^\infty \gamma_n^2(k\semi z),} ]where [Re(\tilde{\lambda})] and [Im(\tilde{\lambda})] denote the real and imaginary parts of [\tilde{\lambda}]. Since [\textstyle\sum_{k=-\infty}^\infty \gamma_n^2(k;z) > 0 ], it is clear that we can set [Im(\tilde{\lambda})] to zero, that is, restrict ourselves to real [\tilde{\lambda}]. Then, minimizing over real [\lambda] we have[\tilde{\lambda}_n(z) = \gamma_n(0\semi z) \bigg/ \textstyle\sum\limits_{k=-\infty}^\infty \gamma_n^2(k\semi z), ]with [\gamma_n(k;z)] given in equation ([link].

Thus we arrive at the least-squares optimal approximation to [\exp(2 \pi i z u) ] by setting [\lambda_n(z) = \tilde{\lambda}_n(z)\bar{\lambda}_n(z) ] in equation ([link]. The mean-square error due to this can be seen to be[\eqalign{{\rm MSE} &= 1 - \gamma_n^2(0\semi z) \bigg/ \textstyle\sum\limits_{k=-\infty}^\infty \gamma_n^2(k\semi z) \cr &= \textstyle\sum\limits_{k \ne 0} \gamma_n^2(k\semi z) \bigg/ \textstyle\sum\limits_{k=-\infty}^\infty \gamma_n^2(k\semi z)} ]and thus the error in the approximation is seen to be small for small [|z|] but increasing as [|z| \to \textstyle{1\over2}], as discussed above (the reciprocal Ewald weight term compensates for the increased error by rapidly damping the larger [|z|] terms). For [0 \,\lt\, |z| \,\lt\, \textstyle{1\over2}] the error decreases to zero as the B-spline order [n \to \infty]. The accuracy of the PME is thus increased by raising the order of interpolation and/or the grid size [M_1,M_2,M_3], which lowers [|z|] for a given m. When derivatives with respect to w of [g_n(w;z)] are considered in this latter representation, their effect is to effectively lower the order of interpolation in the coefficients [\gamma_n]. Thus the order of interpolation needs to be higher, the higher the order of differentiation considered. The fast Fourier Poisson approach

| top | pdf |

The fast Fourier Poisson (FFP) method (York & Yang, 1994[link]) accelerates the calculation of the structure factors in equations ([link] and ([link] in much the same way as FFT methods accelerate structure-factor and density-map calculations in macromolecular structure determination. Originally the method was used to calculate structure factors appearing in Ewald sums over point charges. It is easier to see the basic idea if we rewrite the Ewald reciprocal sum from equation ([link] as follows:[\eqalign{E_{\rm r}(\underline{\bf r}^{\lbrace\bf N\rbrace}) &= {{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 &= {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} {{1}\over{{\bf m}^2}} \tilde{S}({\bf m}) \tilde{S}(-{\bf m}),} ]where the modified structure factor [\tilde{S}({\bf m})] is given by[\tilde{S}({\bf m}) = \exp(-\pi^2 {\bf m}^2 / 2 \alpha)\textstyle\sum\limits_{j=1}^N q_j \exp(2 \pi i {\bf m}\cdot{\bf r}_j). ]

We recognize [\tilde{S}({\bf m})] as the three-dimensional Fourier transform, evaluated at m, of a density [\tilde{\rho}({\bf r})] given by[\tilde{\rho}({\bf r}) = \textstyle\sum\limits_{j=1}^N q_j \rho_{2 \alpha}({\bf r} - {\bf r}_j), ]where [\rho_{2 \alpha}] is a normalized Gaussian with Gaussian exponent 2α.

In turn we approximate this integral by a discrete Fourier transform. As above let [M_1,M_2,M_3] be positive integers of the scale of the unit-cell dimensions. Let [\delta v = V / (M_1 M_2 M_3)], where V is the volume of the unit cell. Recall the lattice basis vectors [{\bf a}_\alpha ], [\alpha=1,2,3]. Given integers [k_1,k_2,k_3] define [\tilde{\rho}_0(k_1,k_2,k_3)] by [\tilde{\rho}_0(k_1,k_2,k_3) = \tilde{\rho}[(k_1/M_1) {\bf a}_1 + (k_2/M_2) {\bf a}_2 + (k_3/M_3) {\bf a}_3] ]. Then we can approximate the modified structure factor [\tilde{S}({\bf m}) ] by[\eqalign{\tilde{S}({\bf m}) &\simeq \sum_{k_1 = -{({M_1}/{2})} + 1}^{{{M_1}/{2}}} \ \sum_{k_2 = -{({M_2}/{2})}+1}^{{{M_2}/{2}}} \ \sum_{k_3 = -{({M_3}/{2})}+1}^{{{M_3}/{2}}} Q(k_1,k_2,k_3) \cr &\quad\times \exp\left[2 \pi i \left({{m_1 k_1}\over{M_1}} + {{m_2 k_2}\over{M_2}} + {{m_3 k_3}\over{M_3}}\right)\right],} ]where [Q(k_1,k_2,k_3)] is given by[Q(k_1,k_2,k_3) = \textstyle\sum\limits_{l_1,l_2,l_3=-\infty}^\infty \tilde{\rho}_0(k_1+l_1M_1,k_2+l_2M_2,k_3+l_3M_3) \cdot \delta v. ]

The FFP can be generalized to Hermite Gaussians by sampling that density in place of the superposition of spherical Gaussians. In equation ([link] the term[\eqalign{&({1}/{2 \pi V}) \textstyle\sum\limits_{{\bf m} \ne {\bf 0}}({1}/{{\bf m}^2})\exp(-\pi^2 {\bf m}^2 /2 \beta) \textstyle\sum\limits_{l_1\in C}S_{l_1}({\bf m}) \cr&\quad\times\exp(-\pi^2 {\bf m}^2 /2 \beta) \textstyle\sum\limits_{l_2\in C} S_{l_2}(-{\bf m})} ]can be handled by the methods just outlined. However the term[\eqalign{& ({1}/{2 \pi V})\textstyle\sum\limits_{{\bf m} \ne {\bf 0}} ({1}/{{\bf m}^2})\textstyle\sum\limits_{(l_1,l_2)\not\in C\times C} \exp(-\pi^2 {\bf m}^2 / \theta_{l_1}) S_{l_1}({\bf m}) \cr &\quad\times \exp(-\pi^2 {\bf m}^2 / \theta_{l_2}) S_{l_2}(-{\bf m})} ]can present problems. The difficulty comes in the case when, for example, [l_1] is compact but [l_2] is not. Then the combined Gaussian weights are sufficiently damping to efficiently converge the series, but it can be very inefficient to directly sample the superposition of Hermite Gaussians with a compact exponent. A similar problem was faced by Fusti-Molnar & Pulay (2002[link]). Their solution was to filter out high frequencies in the Fourier-transformed compact densities and back transform them to obtain diffuse real-space densities. Instead, Cisneros et al. (2006[link]) adjusted the exponents of both the compact and diffuse Gaussians. Choose a [\bar{\theta}] such that if [\theta_l] is diffuse [1/\theta_l > 1 / \bar{\theta}]. For example, choose [\bar{\theta} = 1 / (4 \beta) ]. As before, define [\tilde{S}_l({\bf m})] by [\tilde{S}_l({\bf m})\ = ] [\exp(-\pi^2 {\bf m}^2 / \theta_l) S_l({\bf m}) ], and [\tilde{S}_{l+}({\bf m}),\tilde{S}_{l-}({\bf m})] by[\eqalign{\tilde{S}_{l+}({\bf M}) &= \exp\left[-\pi^2 {\bf m}^2 \left({{1}\over{\theta_l}}+{{1}\over{4 \beta}}\right)\right] S_l({\bf m})\cr \tilde{S}_{l-}({\bf M}) &= \exp\left[-\pi^2 {\bf m}^2 \left({{1}\over{\theta_l}}-{{1}\over{4 \beta}}\right)\right] S_l({\bf m}).} ]Then, if [l_1 \in C] and [l_2 \in D],[\textstyle\sum\limits_{l_1 \in C} \textstyle\sum\limits_{l_2 \in D} \tilde{S}_{l_1}({\bf m}) \tilde{S}_{l_2}({\bf m}) = \textstyle\sum\limits_{l_1 \in C} \textstyle\sum\limits_{l_2 \in D} \tilde{S}_{l_1+}({\bf m}) \tilde{S}_{l_2-}({\bf m}). ]Thus instead of approximating [\tilde{S}_l({\bf m})] by sampling the Gaussian followed by applying the 3DFFT, they approximated [\tilde{S}_{l+}({\bf m}) ] or [\tilde{S}_{l-}({\bf m})]. Note that with the above choice of [\bar{\theta}], for any [l \in C], [\theta_{l+} \,\lt\, 4 \beta ] and for any [l \in D], [\theta_{l-} \,\lt\, 4 \beta]. That is, all the modified Gaussians can be sampled on a uniform grid. Finally, note that for the diffuse–diffuse interactions[\eqalign{&\textstyle\sum\limits_{l_1 \in D} \textstyle\sum\limits_{l_2 \in D} \tilde{S}_{l_1}({\bf m}) \tilde{S}_{l_2}({\bf m})\cr &\quad= \textstyle\sum\limits_{l_1 \in D} \textstyle\sum\limits_{l_2 \in D} \tilde{S}_{l_1-}({\bf m}) \tilde{S}_{l_2-}({\bf m}) \exp(-\pi^2 {\bf m}^2 / 2 \beta),} ]i.e. the interaction between adjusted diffuse Gaussians can be corrected efficiently in reciprocal space.

To test the efficiency of these reciprocal-space approaches, the electrostatic energy and intermolecular forces were calculated for a series of water boxes (64, 128, 256, 512 and 1024 water molecules) under periodic boundary conditions. Similar calculations were performed by Cisneros et al. (2006[link]). As in that paper, the electron density of the water molecule was modelled by a five-site expansion consisting of the oxygen, the two hydrogens and the two bond midpoints, and using standard auxiliary basis functions, denoted A2 and P1, from the density-functional literature. These basis functions were obtained from the Environmental Molecular Sciences Laboratory at Pacific Northwest Laboratories ( ). The A2 basis set expansion uses 21 Hermite primitives with 102 Hermite coefficients per water molecule, resulting in a root-mean-square relative molecular force error of about 8% compared to 6–31G* B3LYP intermolecular electrostatics (Cisneros et al., 2006[link]). The more complex P1 basis set expansion uses 38 Hermite primitives with 200 Hermite coefficients per water molecule, reducing the relative molecular force error to about 2%. For each water box the electrostatic energy and forces for the model A2 or P1 densities were calculated to a relative root-mean-square molecular force error of 10−3 or 10−4, compared to exact Ewald summation of the respective model density.

The time to perform these calculations is given in Table[link]. Note that the times given here are considerably faster than those reported in Figs. 13 and 14 of Cisneros et al. (2006[link]). This is partly due to the faster processor used here (3.6 GHz Intel Xeon versus 3.2 GHz Xeon used previously; Intel compiler with similar options used in both cases) but mainly due to the more efficient algorithms described herein. Note that the PME and FFP timings increase in an approximately linear fashion as a function of box size whereas the Ewald sum timings behave quadratically. For each method the timings are a continuous function of the compact–diffuse separating exponent 2β described above in Section[link], except when 2β equals a Gaussian exponent in the model basis set. When this happens, a set of Hermite primitives change from compact to diffuse or vice versa and the relative balance of direct-sum or reciprocal-sum calculation changes discontinuously. For the smooth PME calculations it was optimal to choose the separating exponent at one of these boundaries. This separating exponent value was then retained for the FFP and Ewald calculations, although the FFP method could be made slightly more efficient by further varying this parameter, and the Ewald summation could be considerably accelerated to an order N3/2 method, where N is the number of water molecules, by making the separating exponent a decreasing function of system size. Neither of these optimizations would change the relative efficiencies, however.

Table| top | pdf |
Timing versus accuracy for smooth PME, FFP and regular Ewald sum reciprocal-space approaches

(a) Time in seconds for A2 model density.

Box sizeRMS force error 10−3RMS force error 10−4
64 0.106 0.142 0.365 0.144 0.274 0.520
128 0.218 0.271 1.336 0.287 0.380 1.869
256 0.387 0.528 5.239 0.517 0.846 7.256
512 0.837 1.100 17.881 1.104 1.534 25.511
1024 1.701 2.236 71.513 2.249 3.152 108.158

(b) Time in seconds for P1 model density

Box sizeRMS force error 10−3RMS force error 10−4
64 0.310 0.399 2.321 0.478 0.688 3.858
128 0.591 0.832 7.706 0.923 1.576 11.056
256 1.186 1.920 35.178 1.805 2.736 49.107
512 2.549 3.825 119.863 3.794 5.684 183.487
1024 4.953 6.890 486.384 6.947 11.307 714.644

The error in the FFP for a given Gaussian exponent is controlled by the size of the grid [M_1,M_2,M_3] as in the PME method, and also by the completeness in the Gaussian sampling (the cutoff in the tails of the Gaussian). Typically, for equivalent accuracy, the FFP requires more extensive sampling onto grid points than does the smooth PME. That is, each Gaussian centred at an expansion point must be sampled over a considerably larger set of points than the tensor product B-spline in equation ([link]. The additional cost due to this has been ameliorated somewhat by the Gaussian split Ewald approach (Shan et al., 2005[link]). Additionally, in rectangular unit cells fast Poisson solvers such as multigrid (Beck, 2000[link]; Sagui & Darden, 2001[link]) can be applied, which in principle should parallelize more efficiently than the 3DFFT. Finally, although the smooth PME is typically more efficient for a single Gaussian exponent structure factor, the FFP approach has the strong advantage that the various Gaussian superpositions having possibly different exponents can be sampled and then fast Fourier transformed together, whereas they must be kept separate in the PME approach. Thus, as here Cisneros et al. (2006[link]) found that the smooth PME approach was more efficient for smaller auxiliary basis sets, whereas they found that the FFP approach became more efficient for large auxiliary basis sets involving many different Gaussian exponents.


Aguado, A. & Madden, P. A. (2003). Ewald summation of electrostatic multipole interactions up to the quadrupolar level. J. Chem. Phys. 119, 7471–7483.
Allen, M. P. & Tildesley, D. J. (1987). Computer Simulation of Liquids. Oxford: Clarendon Press.
Arfken, G. B. & Weber, H. J. (2000). Mathematical Methods for Physicists, 5th ed. San Diego: Academic Press.
Beck, T. L. (2000). Real-space mesh techniques in density functional theory. Rev. Mod. Phys. 72, 1041–1080.
Böttcher, C. J. F. (1973). Theory of Electric Polarization. Amsterdam: Elsevier Science Publishers BV.
Broderson, S., Wilke, S., Leusen, F. J. L. & Engel, G. (2003). A study of different approaches to the electrostatic interaction in force field methods for organic crystals. Phys. Chem. Chem. Phys. 5, 4923–4931.
Brown, D. & Neyertz, S. (1995). A general pressure tensor calculation for molecular dynamics simulations. Mol. Phys. 84, 577–595.
Challacombe, M., White, C. & Head-Gordon, M. (1997). Periodic boundary conditions and the fast multipole method. J. Chem. Phys. 107, 10131–10140.
Chemburkar, S. R., Bauer, J., Deming, K., Spiwek, H., Patel, K., Morris, J., Henry, R., Spanton, S., Dziki, W., Porter, W., Quick, J., Bauer, P., Donabauer, J., Narayanan, B. A., Soldani, M., Riley, D. & McFarland, K. (2000). Dealing with the impact of ritonavir polymorphs on the late stages of bulk drug process development. Org. Process Res. Dev. 4, 413–417.
Cisneros, G. A., Piquemal, J. P. & Darden, T. A. (2005). Intermolecular electrostatics using density fitting. J. Chem. Phys. 123, 044109.
Cisneros, G. A., Piquemal, J. P. & Darden, T. A. (2006). Generalization of the Gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods. J. Chem. Phys. 125, 184101.
Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.
Coppens, P. & Volkov, A. (2004). The interplay between experiment and theory in charge-density analysis. Acta Cryst. A60, 357–364.
Darden, T., Pearlman, D. A. & Pedersen, L. G. (1999). Ionic charging free energies: Spherical versus periodic boundary conditions. J. Chem. Phys. 109, 10921–10935.
Darden, T. A., Toukmaji, A. & Pedersen, L. (1997). Long-range electrostatic effects in biomolecular simulations. J. Chim. Phys. 94, 1346–1364.
Darden, T. A., York, D. M. & Pedersen, L. G. (1993). Particle mesh Ewald: An N log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092.
Day, G. M., Chisholm, J., Shan, N., Motherwell, W. D. S. & Jones, W. (2004). An assessment of lattice energy minimization for the prediction of molecular organic crystal structures. Cryst. Growth Des. 4, 1327–1340.
Day, G. M., Motherwell, W. D. S., Ammon, H. L., Boerrigter, S. X. M., Della Valle, R. G., Venuti, E., Dzyabchenko, A., Dunitz, J. D., Schweizer, B., van Eijck, B. P., Erk, P., Facelli, J. C., Bazterra, V. E., Ferrarro, M. B., Hofmann, D. W. M., Leusen, F. J. J., Liang, C., Pantelides, C. C., Karamertzanis, P. G., Price, S. L., Lewis, T. C., Nowell, H., Torrisi, A., Scheraga, H. A., Arnautova, Y. A., Schmidt, M. U. & Verwer, P. (2005). A third blind test of crystal structure prediction. Acta Cryst. B61, 511–527.
Day, G. M., Motherwell, W. D. S. & Jones, W. (2005). Beyond the isotropic atom model in crystal structure prediction of rigid molecules: Atomic multipoles versus point charges. Cryst. Growth Des. 5, 1023–1033.
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.
DeLeeuw, S. W., Perram, J. W. & Petersen, H. G. (1990). Hamilton's equations for constrained dynamical systems. J. Stat. Phys. 61, 1203–1222.
Deserno, M. & Holm, C. (1998). How to mesh up Ewald sums I: A theoretical and numerical comparison of various particle mesh routines. J. Chem. Phys. 109, 7678–7693.
Destro, R., Roversi, P., Barzaghi, M. & Marsh, R. E. (2000). Experimental charge density of α-glycine at 23 K. J. Phys. Chem. A, 104, 1047–1054.
Dovesi, R., Orlando, R., Civalleri, B., Roetti, C., Saunders, V. R. & Zicovich-Wilson, C. M. (2005). Crystal: a computational tool for the ab initio study of the electronic properties of crystals. Z. Kristallogr. 220, 571–573.
Dunitz, J. D. (2003). Are crystal structures predictable? Chem. Commun. 5, 545–548.
Dunitz, J. D. & Bernstein, J. (1995). Disappearing polymorphs. Acc. Chem. Res. 28, 193–200.
Dunitz, J. D. & Gavezzotti, A. (2005). Molecular recognition in organic crystals: Directed intermolecular bonds or nonlocalized bonding? Angew. Chem. 44, 1766–1787.
Dunitz, J. D. & Scheraga, H. A. (2004). Exercises in prognistication: Crystal structures and protein folding. Proc. Natl Acad. Sci. USA, 101, 14309–14311.
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.
Frenkel, D. & Smit, B. (2002). Understanding Molecular Simulation. San Diego: Academic Press.
Fusti-Molnar, L. & Pulay, P. (2002). The Fourier transform Coulomb method: Efficient and accurate calculation of the Coulomb operator in a Gaussian basis. J. Chem. Phys. 117, 7827–7835.
Gavezzotti, A. (2002a). Structure and intermolecular potentials in molecular crystals. Modelling Simul. Mater. Sci. Eng. 10, R1–R29.
Gavezzotti, A. (2002b). Calculation of intermolecular interaction energies by direct numerical interaction energies I. Electrostatic and polarization energies in molecular crystals. J. Phys. Chem. B, 106, 4145–4154.
Gavezzotti, A. & Fillipini, G. (1994). The geometry of the intermolecular x—H...y (x, y = N, O) hydrogen bond and the calibration of empirical hydrogen bond potentials. J. Phys. Chem. 98, 4831–4837.
Gill, P. M. W. & Adamson, R. D. (1996). A family of attenuated Coulomb operators. Chem. Phys. Lett. 261, 105–110.
Helgaker, T., Jorgensen, P. & Olsen, J. (2000). Molecular Electronic-Structure Theory. Chichester: John Wiley and Sons.
Hockney, R. W. & Eastwood, J. W. (1981). Computer Simulation Using Particles. New York: McGraw-Hill.
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.
McMurchie, L. E. & Davidson, E. R. (1978). One-electron and 2-electron integrals over Cartesian Gaussian functions. J. Comput. Phys. 26, 218–231.
Mooji, W. T. M., van Duijneveldt, F. B., van Duijneveldt-van de Rijdt, J. G. C. M. & van Eijck, B. P. (1999). Transferable ab initio intermolecular potentials. 1. Derivation from methanol dimer and trimer calculations. J. Phys. Chem. A, 103, 9872–9882.
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.
Nowell, H., Frampton, C. S., Waite, J. & Price, S. L. (2006). Blind crystal structure prediction of a novel second polymorph of 1-hydroxy-7-azabenzotriazole. Acta Cryst. B62, 642–650.
Perram, J. W., Petersen, H. G. & DeLeeuw, S. W. (1988). An algorithm for the simulation of condensed matter that grows as the 3/2 power of the number of particles. Mol. Phys. 65, 875–893.
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.
Piquemal, J. P., Cisneros, G. A. & Darden, T. A. (2006). Towards a force field based on density fitting. J. Chem. Phys. 124, 104101.
Pollock, E. L. & Glosli, J. (1996). Comments on p(3)m, fmm and the Ewald method for large periodic Coulombic systems. Comput. Phys. Comm. 95, 93–110.
Price, S. L. & Price, L. S. (2005). Modelling intermolecular forces for organic crystal structure prediction. Struct. Bond. 115, 81–123.
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.
Ren, P. Y. & Ponder, J. W. (2003). Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B, 107, 5933–5947.
Sagui, C. & Darden, T. (2001). Multigrid methods for classical molecular dynamics simulations of biomolecules. J. Chem. Phys. 114, 6578–6591.
Sagui, C., Pedersen, L. G. & Darden, T. A. (2004). Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulations. J. Chem. Phys. 120, 73–87.
Schoenberg, I. J. (1973). Cardinal Spline Interpolation. Philadelphia: Society for Industrial and Applied Mathematics.
Shan, Y. B., Klepeis, J. L., Eastwood, M. P., Dror, R. O. & Shaw, D. E. (2005). Gaussian split Ewald: A fast Ewald mesh method for molecular simulations. J. Chem. Phys. 122, 054101.
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. (1982). Point multipoles in the Ewald summation. CCP5 Inf. Q. 4, 13–25.
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.
Smith, W. (1998). Point multipoles in the Ewald summation (revisited). CCP5 Inf. Q. 46, 18–30.
Spackman, M. A., Weber, H. P. & Craven, B. M. (1988). Energies of molecular-interactions from Bragg diffraction data. J. Am. Chem. Soc. 110, 775–782.
Stone, A. J. (1996). The Theory of Intermolecular Forces. New York: Oxford University Press.
Trickey, S. B., Alford, J. A. & Boettger, J. C. (2004). Methods and implementation of robust, high-precision Gaussian basis DFT calculations for periodic systems: the GTOFF code. In Theoretical and Computational Chemistry, edited by J. Leszcynski, Vol. 15, ch. 6. Amsterdam: Elsevier.
Volkov, A., Koritsansky, T. & Coppens, P. (2004). Combination of the exact potential and multipole methods (ep/mm) for evaluation of intermolecular interaction energies with pseudoatom implementation of electron densities. Chem. Phys. Lett. 391, 170–175.
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.
Williams, D. E. (2001). Improved intermolecular force field for molecules containing H, C, N and O atoms, with application to nucleoside and peptide crystals. J. Comput. Chem. 22, 1154–1166.
Williams, D. E. & Cox, S. R. (1984). Nonbonded potentials for azahydrocarbons: the importance of the Coulombic interactions. Acta Cryst. B40, 404–417.
York, D. & Yang, W. (1994). The fast Fourier Poisson (FFP) method for calculating Ewald sums. J. Chem. Phys. 101, 3298–3300.

to end of page
to top of page