International
Tables for Crystallography Volume B Reciprocal space Edited by U. Shmueli © International Union of Crystallography 2010 
International Tables for Crystallography (2010). Vol. B, ch. 3.5, pp. 458481
doi: 10.1107/97809553602060000772 Chapter 3.5. Extensions of the Ewald method for Coulomb interactions in crystals^{a}Laboratory of Structural Biology, National Institute of Environmental Health Sciences, 111 T. W. Alexander Drive, Research Triangle Park, NC 27709, USA 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 Ewaldsum technique (assuming a neutral unit cell). The form of this sum is derived here in a manner that generalizes to other interactions, such as longrange 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 nearlinear scaling (i.e. scaling as N log N with system size N) are derived and tested. 
Highprecision singlecrystal Xray structural analysis of small organic molecules, yielding the space group, the unitcell 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 ΔH_{sub}, 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 relationshipwhere the lattice or packing energy E_{lattice} is the total (molar) intermolecular interaction energy between all the molecules in the crystal, which are treated as rigid entities with zeropoint 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) and (in more detail) in Frenkel & Smit (2002).
Thus, given a parameterized intermolecular potentialenergy function, or if computationally affordable a firstprinciples approach such as densityfunctional theory (or preferably, when it becomes feasible for crystals, a goodquality postHartree–Fock potentialenergy 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 crystalstructure parameters, one could in principle predict the minimumlatticeenergy 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 minimumenergy structure (i.e. it is not a metastable state determined by crystalgrowth 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; Dunitz & Scheraga, 2004). 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). 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 bioavailability of various forms of a drug such as ritonavir (Chemburkar et al., 2000). 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 lowestfreeenergy polymorph at given temperature and pressure may not be the likeliest to form. The kinetics of growth of microcrystals may largely determine which lowenergy polymorph appears (Dunitz, 2003). 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 crystalstructure prediction has been undertaken (Day, Motherwell, Ammon et al., 2005). The results of these tests have highlighted the need for continued improvements in sampling methods and intermolecular energy potentials. Since extensive sampling of the crystalstructure parameters is necessary [between 10^{4} and 10^{5} starting structures, each followed by parameter minimizations (Price & Price, 2005)], there is a tradeoff 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 slowdecaying interaction terms that occur in these potential functions.
Methodological developments in the intermolecular force fields used in crystalstructure prediction from early times to the present state of the art have been reviewed (Price & Price, 2005). Until recently, these force fields were made up of atom–atom interactions. The earliest involved only repulsion and dispersion, usually in the `exp6' formwhere U^{MN} denotes the intermolecular potential energy between molecules M and N and r_{ij} is the distance between atoms and . 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), 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 shortranged 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 longrange 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 exp6 potential in fitting organic crystal structures with and without hydrogenbond interactions (Gavezzotti & Fillipini, 1994). 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 exp6 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) to supplement these with selected point charges, both atomic and at offatom sites. Price and coworkers 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). 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). The AMOEBA force field has similar capabilities (Ren & Ponder, 2003).
Recently, a systematic study (Day et al., 2004; Day, Motherwell & Jones, 2005) of the effect of the electrostatic model on crystalstructure prediction for rigid organic molecules was carried out. The exp6 parameters from the W99 force field (Williams, 2001) were used in conjunction with point charges derived either from a bondincrement model or electrostatic potential (ESP) fit, or with a DMA model as used by Price and coworkers. Importantly, they found a systematic and significant improvement in prediction success going from the simple bondincrement model to the ESP charges to point multipoles. Similar results were found by Broderson et al. (2003). 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), crystalpacking 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 gasphase 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 gasphase intramolecular energy and DMA electrostatic model are calculated using quantumchemical 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 1hydroxy7azabenzotriazole (Nowell et al., 2006).
The above electrostatic models are all derived from gasphase quantumchemical calculations. A different line of development has been the parameterization of electrostatic models from condensedphase 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 Xray charge densities (Destro et al., 2000) or of theoretical charge densities from quantumchemical calculations carried out under the spacegroup symmetry (Dovesi et al., 2005). Crystal structure analysis is usually carried out assuming the independentatom 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 atomcentred 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) and Chapter 1.2 of this volume. The theoretical density can be similarly approached as a sum of promolecule density (superposition of isolated atomic densities in the molecular geometry) plus deformation density modelled by a sum of atomic Slaterlike 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). Coppens and coworkers as well as others (Spackman et al., 1988) have applied these modelled densities in studies of intermolecular electrostatic energies. The Slaterlike 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).
A similar densitybased approach using gasphase theoretical density was developed previously (Gavezzotti, 2002b). 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 pixelbased 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).
Cisneros et al. have also begun implementing an intermolecular potential based on the explicit modelling of gasphase theoretical electron density, using the methodology of auxiliary basis set fitting (Cisneros et al., 2005). The density is expanded into atom and offatom 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 nonisotropic site–site repulsion based on density overlap have been implemented (Piquemal et al., 2006) and efficient calculation of terms in periodic boundary conditions has also been implemented (Cisneros et al., 2006), generalizing previous developments for point multipoles (Sagui et al., 2004).
Summarizing the above developments, it appears that accurate evaluation of crystal lattice energies necessitates the use of complex electrostatic interactions beyond the sphericalatom 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 longranged, necessitating an accurate treatment of their lattice summation. Dispersion interactions, modelled by higherorder inverse powers of distance, are shortranged but their calculation can be greatly accelerated through latticesummation techniques.
Ewald summation in ideal infinite crystals, for interactions that depend on inverse powers of site–site distance, is covered in Chapter 3.4 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 latticesummation 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 physicsbased 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 higherorder 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 (3.5.2.19), that separates 1/r^{n} into shortranged and longranged components, is the same as that used in Chapter 3.4 of this volume, highlighting the essential similarity of our approach. In the second section we generalize these pointcharge results to lattice sums of interacting spherical Gaussian and Hermite Gaussianbased 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.
We begin by discussing an idealized infinite crystal C, made up of point charges.
The periodicity of an idealized crystal C is specified by lattice basis vectors , and , 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 are defined by the relations (the Kronecker delta) for . Thus for example is given by = . General lattice vectors n are given by integral combinations of the lattice basis vectors:where are integers. General reciprocallattice vectors m are given by integral combinations of the reciprocallattice basis vectors:
An arbitrary point r in the crystal (r also denotes the vector from the origin of coordinates to the point) is given bywhere are the fractional coordinates of the point. From the above equation and the defining relations for the reciprocal basis we see that for . The unit cell U is the set of points r with associated fractional coordinates satisfying , . 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 at positions for all lattice vector n. The periodic image of U, denoted , is given by the periodic images of all point charges q in U. Thus, the crystal is made up of the union of images , including , or
A point charge in U at position interacts with other unitcell charges at positions as well as with all of their periodic images at positions for all lattice vectors n. It also interacts with its own periodic images at for all nonzero lattice vectors n, that is all with . Suppose the unit cell U is made up of N point charges at positions . 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 aswhere here and below denotes the set of points 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 ). When U is neutral, the sum, if carried out in a naive fashion, can still converge quite slowly (see Chapter 3.4 ). 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 reexpress the above energy. Let E(U, U) denote the Coulomb interaction energy of U with itself and the Coulomb interaction energy of U with , that isandThen we can abbreviate equation (3.5.2.3) byNote that for large n (that is, n with large norm ) can be approximated by a Taylor expansion. Let denote the centroid of points in U, given by , and let denote the centroid of . Then, expanding in terms of and (which are small compared to ) we have
Let and denote the net charge and dipole moment of the unit cell U, and let = denote the quadrupole moment of the unit cell. Using the expansion in the above equation, we can Taylor expand the interaction energy asNote that the firstorder term in the expansion has dropped out, since the dipole moments D of U and are the same. If then for large the leading term in the Taylor expansion of is . Since , then using a standard integral comparison is infinite, and so the resulting electrostatic energy of U diverges. If we see that the leading term in the expansion of is given by , that is a quadratic expression in the dipole moment D times a prefactor that decays like . Since , 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 . If Q = 0 and D = 0, the leading term in the Taylor expansion of depends on the quadrupole moment of the unit cell and decays like for large , and thus since 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 unitcell 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.
In this section we give the standard Ewaldsum derivation found in textbooks such as Kittel (1986). 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 due to an isolated Gaussian density given byOwing to the spherical symmetry of we can define the onedimensional function by for any r with . We next invoke Gauss' law (Eyges, 1980) under spherical symmetry to evaluate the electric field due to at any distance . 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 were collected into a point charge at the origin. More specifically, if S is the sphere of radius with surface , Gauss' law statesand thus, defining the scalar quantity , we getand thus, since (again by symmetry) , we can verify thatwhere C is an arbitrary constant of integration. We obtain a potential with the desired asymptotic behaviour, i.e. as , by setting C = 0. Note then that for large r as expected.
Next we derive the electrostatic potential 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 bywhere the density is defined in equation (3.5.2.6) and the sum is over all general lattice vectors n, defined in equation (3.5.2.1). Since is periodic, i.e. for any general lattice vector , we can expand it in a Fourier series over the unit cell:where the sum is over all reciprocallattice vectors m, defined in equation (3.5.2.2), and the Fourier coefficient is given bywhere V is the volume of the unit cell U. The periodicity of implies that of , so we can also expand the latter in a Fourier series:Next we recall Poisson's equation (Eyges, 1980) for in terms of :Using equations (3.5.2.10) and (3.5.2.8) we can rewrite equation (3.5.2.11) asthus using the uniqueness of Fourierseries coefficients . If this latter equation can be solved for . However, for m = 0 it requires that , whereas from equation (3.5.2.9) . To avoid this contradiction we modify by subtracting 1/V (uniform neutralizing plasma). That is, is now defined byNote that is still periodic. Also note that for is still given by equation (3.5.2.9), whereas now . Although we have avoided the contradiction, the value of is still not determined. Some authors define it by considering m as a real variable and taking the limit as (this depends on which direction the origin is approached from), but most treatments take it to be 0 (however, see Section 3.5.2.3). In the latter case the electrostatic potential due to is given by
Now as before let be point charges located at positions within the unit cell U. From the discussion in Section 3.5.2.1 we assume the unit cell is neutral, i.e. . Let r be a point in U with . We want to obtain the electrostatic potential due to the point charges in U together with their periodic images in the remainder of the crystal C. That is,We next supplement the potential due to the point charge at by the potential due to an equally charged Gaussian counterion, that is times a density of the form in equation (3.5.2.6), centred at . Using equation (3.5.2.7) the resulting direct sum potential is given bywhere . Note that since decays exponentially fast to zero as , 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 times a density of the form , centred at , for all directlattice vectors n. Since the unit cell is neutral, the sum of these reciprocal potentials over and all n is the same as the sum over of times the potential due to centred at where is given by equation (3.5.2.12). Note that the extra term involving 1/V drops out by neutrality. [In the case of nonneutrality, such as discussion of the free energy of charging an ion (Hummer et al., 1996), 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).] Thus, using equation (3.5.2.13) we arrive at the reciprocal sum potential:where the structure factor is given by
We have now expressed as the sum of plus (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 at the position of a charge in U we must remove the Coulomb potential due to and then take the limit as . Noting that and that , we haveThus we can write the electrostatic energy of U from equation (3.5.2.3) aswhere as above the prime indicates that terms with i = j and n = 0 are omitted.
The derivation in Section 3.5.2.2 is incomplete. In addition to the caveats concerning the expression for in terms of plus , there is no mention of a quadratic term involving the unitcell 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, 1994).
The derivation depends on the following identity:where is the reciprocal unit cell, i.e. the set of points with , . Actually we derive (following Essmann et al., 1995) a generalization of identity (3.5.2.16), valid for general inverse powers , p > 0. For this first note that for u > 0where is the Euler gamma function. Next note that for a > 0which is the Fourier transform of the Gaussian, as derived above within equation (3.5.2.9). In equation (3.5.2.17) substitute and u = p/2 to get, for α > 0,In we substitute , whereas in we use equation (3.5.2.18), change the order of integration, and substitute to obtainwhereandSpecializing to the case p = 1 we see that and that , which proves identity (3.5.2.16). Another important case is p = 6, for which = and = .
Note that from equation (3.5.2.17), substituting t = s^{2}, λ = r^{2} and u = p/2 we haveAlso, substituting into equation (3.5.2.20) we haveand so
Note that is a smooth function of x for , and that its derivatives are uniformly bounded for x bounded away from 0. To show this, substitute s = xt in equation (3.5.2.21) to write
In the case that p = 6, 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, is bounded and continuous as , whereas for as . Thus for p > 3 we see that for any x > 0, . However for , for any R > 1 if x > 0 sufficiently small ; thus and thus is unbounded as .
Despite this latter result, the integral of over is bounded for . This result is important below for the validity of rearranging the order of summation or taking limits inside infinite sums involving . To prove boundedness, use the above expression for and integrate by parts to write
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 integralfor some K > 0. Although the sum can thus be evaluated directly, it can be significantly accelerated using the same Ewald direct and reciprocalspace decomposition. For , using equation (3.5.2.19) we havewhere . Since is bounded, the first of the above sums converges absolutely and we turn to the second. We recognize the integral as (1/V) times , the nth Fourier coefficient of , where V is the volume of the unit cell U (1/V is the volume of ). From the above results for it is clear that for , are smooth and have uniformly (in v, r and ) bounded derivatives. Thus the sum of the Fourier coefficients can be shown to converge absolutely and uniformly in m, which guarantees uniform (over ) convergence of the sum of these integrals to , which in turn equals .
The case of m = 0 is not as straightforward, since the uniform boundedness of derivatives of does not generally hold in this case. One way of handling it is as follows. Since everything else in the above equation converges as , the sum of the Fourier coefficients for converges to something. That is, iffor , , then for some s, as . Ifthen it is straightforward to show that also converges to s as [see, for example, Chapter 1 of Körner (1988) for the onedimensional proof]. On the other hand, due to the boundedness and continuity of , as . [This involves the threedimensional version of the Fejér kernel, which converges to the threedimensional Dirac delta function. See, for example, Chapter 2 of Körner (1988) for the onedimensional proof.] Thus as . Hence the above sum of integrals with m = 0 converges to .
Thus from the above arguments we have
Using the limiting behaviour given in (3.5.2.22) we can extend this result to the case r = 0. By removing the singularity when n = 0 we get
Now assume particles i and j in the unit cell U interact by the potential . Using equations (3.5.2.24) and (3.5.2.25) (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 with each other and all of the image cells is given bywhere as before the prime on the first summation means that terms where i = j and n = 0 are omitted.
In the case when factorizes, i.e. , the above lattice sum can be further simplified and accelerated using the structure factor given byIn this case,
Note that as expected, due to the absolute convergence of the lattice sums, there is no shapedependent correction term to the Ewald sum. This is connected with the fact that is bounded and continuous as . 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.
In the case , the lattice sum does not converge absolutely, since the integral in equation (3.5.2.23) 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 centred on the origin, such as a cube, sphere or ellipsoid. For a positive integer K, let denote the set of lattice vectors such that is in P, that is . Charges in the central unit cell U_{0} interact with each other and with image charges in U_{n} for . The electrostatic energy of the central cell in this finite crystal is given byand we want the limit of as .
Recalling that for p = 1, if we define bythen for , is bounded and smooth with uniformly (in ) bounded derivatives. For m = 0 we have = , where is smooth with bounded derivatives, and with . Then, using the same arguments that lead to equation (3.5.2.24), we can write for where here and below denotes a quantity (the remainder in the equation) that converges to zero as and
The Ewald potential , given for byhas a couple of noteworthy properties. First, note that since the lefthand side of equation (3.5.2.28) is independent of α, as is , is itself independent of α. Secondly, the average of over the unit cell U is zero. To establish this latter property, first note that for , . Thus the reciprocal sum in the above equation integrates to zero over U. Next, the integral of −π/(αV) over U is clearly −π/α. Finally, note thatcompleting the proof. The self or Wigner potential with given byis 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, (Nijboer & Ruijgrok, 1988). Owing to the above properties, the Ewald potential provides a consistent approach to the treatment of nonneutral unit cells, such as when calculating ionic charging free energies (Hummer et al., 1996; Darden et al., 1999).
If we define the structure factor byand the unitcell dipole moment D bythen assuming the unit cell is neutral we can write the energy asThus, recalling equation (3.5.2.15), we can writewhereNote that if D = 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 have been derived in some cases by Smith (1981). In particular, if P is the unit sphere . Following Smith's approach we show this by first deriving an integral representation for and then evaluating it for the case P is the unit sphere.
First note that for any continuous function f, converges to as . Then, in the above integral for we can substitute w = Kv to getand thus
Specializing to the case P = S, the unit sphere, we see by using spherical coordinates for with the pole parallel to w that is a function only of the norm of w, that is . Then, again using spherical coordinates, this time for w with the pole parallel to D, and noting that , we havewhere we have switched the order of integration in the last step. Finally, since , the threedimensional 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.
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), we develop a surfaceintegralbased expression for the shapedependent correction term for more general macroscopic shapes.
Given a closed bounded region P that contains the unit sphere S in its interior and a K > 0, and recalling equations (3.5.2.4) and (3.5.2.27), we can write the electrostatic energy of the central unit cell in the finite crystal as
The asymptotic form of has been developed above. Recalling equation (3.5.2.5), for large K and neutral unit cells, we can approximate bywhere we have used Gauss' divergence theorem (Arfken & Weber, 2000) together with
Since, considering the outward unit normal on S,and since from the previous resultswe can now generalize the lattice summation to asymptotic shapes P:
The latter surface integral is referred to as the `surface term' and clearly must equal the quantity discussed above. It was derived previously by Deem et al. (1990) 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 surfaceintegral contribution (see also Smith, 1981). In contrast, Scheraga and coworkers argue for the importance of considering this term, which they refer to as the `Lorentz–Ewald correction', in crystalstructure prediction (Pillardy et al., 2000), although in this they are disputed by van Eijck & Kroon (2000) (see also the response by Wedemeyer et al., 2000).
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).
In this section, following the approach in Smith (1981, 1994), 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 , assuming the macroscopic shape of the crystal is spherical. We first derive the dielectric reaction potential at a point due to the charges and their periodic images in , assuming the unit cell is neutral, i.e. , that P = S, the unit sphere, and that the sphere KS is immersed in the dielectric.
To begin with, given a point and a point , with spherical polar coordinates and , respectively, the electrostatic potential at r due to a unit charge at is given by the solution to Laplace's equation. Let β be the angle between r and , that is . Owing to the axial symmetry, the general solution to Laplace's equation in this case can be writtenwhere are the Legendre polynomials (Böttcher, 1973). The contributions can be thought of as due to charges further from the origin than r, while are due to charges closer to the origin. A particular solution to Laplace's equation is given by . Thus, we can write the potential due to the unit charge at in the presence of the dielectric continuum surrounding KS asWe now apply boundary conditions, noting that are an orthogonal basis of functions, so that boundary conditions can be applied term by term. Then we have
We expand (for ) in Legendre polynomials to allow comparison of with as from below:Using this we get:
Finally, since we need the solution for the potential due to each of the charges in , we invoke the spherical harmonic addition theorem (Arfken & Weber, 2000) applying it to the points r, :where (Arfken & Weber, 2000)and the associated Legendre polynomials , are defined by Rodrigues' formula,andIn particular, the firstorder spherical harmonics are given byThe spherical harmonics satisfy an orthogonality relation:Then, for we express the polarization potential due to the dielectric response to a unit charge at aswhereand
We now examine the total polarization potentialwhere D is the unitcell dipole and we have used unitcell neutrality in the next to last equation.
The gradient with respect to the usual Cartesian directions, expressed in spherical polar coordinates (Arfken & Weber, 2000) can be reexpressed as
Using spherical polar coordinates to carry out the integrations over the unit sphere in equation (3.5.2.35), we see after the first partial integration over thatWe recognize the factors multiplying above as being proportional to , and , respectively, so that, using the above orthogonality of spherical harmonics, these contribute nothing to the integrals unless l = 1.
Next, examining I given byintegrating over and then applying integration by parts over to the first term, expressing using the Rodrigues' formulas (3.5.2.32) and (3.5.2.33), and finally substituting we see that I is proportional to the integral J given bywhich in turn is zero unless l = 1. Putting the above results together we see thatSimilarlyFinally,is proportional toand J is zero unless l = 1. Thus
Expressing in terms of the explicit forms for given in equations (3.5.2.34), we getand thusSimilarlyand finallySubstituting these into equation (3.5.2.35) we get
Recalling the shapedependent term from equation (3.5.2.31) and adding the reaction energy , 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 is given bywhere is given in equation (3.5.2.15). Note that we recover the usual Ewald sum in the limit as (`tinfoil' boundary conditions).
The force on a charge in the central cell is obtained by taking the gradient of the energy with respect to for the finite spherical crystal, and then taking the limit as , or by taking the gradient of the limiting energy above (i.e. the limit of the gradient is the gradient of the limit):
Thus we can find tractable analytic expressions for the energy and force on particles within the central unit cell as long as the asymptotic shape of the crystal immersed in the dielectric medium is spherical. In Section 3.5.2.6 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), which may reflect some not yet fully understood subtleties in the implementation of periodic boundary conditions.
We first discuss the scalar pressure, following the developments in Smith (1987, 1993). According to the thermodynamic definition, the scalar pressure P is defined bywhere A is the Helmholtz free energy, V is the volume and T is the temperature. In turn, A is given by where ( is Boltzmann's constant) and is the canonical ensemble partition function,where is the system Hamiltonian, h is Planck's constant and denotes the set of momenta with conjugate to . Introducing the scaling relations , and , we getwhere denotes the expectation or ensemble average. For the above system of N point charges having mass at positions in the central unit cell of a large spherical crystal (of radius K) immersed in a dielectric continuum, the Hamiltonian can be writtenand soThe direct and reciprocallattice vectors are scaled as and . The derivative with respect to volume of the Ewald direct sum is straightforward using the chain rule, since if then and = . Since m and scale oppositely with V, . Continuing with the other terms in the Ewald reciprocal sum, and taking the limit as , we get
The scalar pressure has been generalized to the pressure tensor (Brown & Neyertz, 1995) using a more general scaling technique based on the expression for r in terms of unitcell 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, 1975) 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 and . The trace of the pressure tensor defined by Nosé & Klein (1983) agrees with the above scalar pressure in the limit .
Interestingly, Smith (1994) derives a scalar pressure that differs from the above. In particular he finds a surface term that remains even in the limit of `tinfoil' boundary conditions, i.e. when . 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). He derives forces on particles within the central unit cell (and on their images in other cells interior to the crystal) that agree in the large K limit with our result above in equation (3.5.2.36). 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 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 . 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). Furthermore, his scalar virial, given by the trace of his virial tensor, has an 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.
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) the electrostatic potential and energy for crystal densities modelled by linear combinations of Slaterlike density elements is discussed. As with quantumchemical programs, there are tradeoffs in the choice of density elements. The use of Slaterlike 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 fastFouriertransformbased methods discussed in Section 3.5.4, 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 densityfunctional calculations, is given in Trickey et al. (2004). The continuous fast multipole method (Challacombe et al., 1997) provides yet a different but highly efficient approach to lattice summation of continuous charge densities.
We begin by discussing the electrostatic interaction energy between a normalized spherical Gaussian charge distribution , centred at a point and having exponent , and a second normalized spherical Gaussian charge distribution , centred at and having exponent , together with the latter's periodic images, centred at , . Their interaction energy can be writtenwhere , j = 1, 2. We now recall equation (3.5.2.16), applying its two righthand components separately, to get , whereand
The direct pair energy involves the damped Coulomb interaction, which has been discussed by Gill & Adamson (1996). We can calculate using the following standard results from quantum chemistry, which are derived in Helgaker et al. (2000): Let , be two normalized spherical Gaussian charge distributions having exponents p and q and centred at points P and Q, i.e. and = . Then the electrostatic potential at a point R due to the charge distribution is given bywhere is the zerothorder Boy's function, with for x > 0. The Coulomb interaction energy between and is given bywhere . The Boy's function can be related to the error function:Substituting equation (3.5.3.4) into (3.5.3.2), we see thata result which we had derived above by another method [see equation (3.5.2.7)]. From equations (3.5.3.3) and (3.5.3.5), substituting for p and for , we have for any point where . Iterating this result with replacing and replacing we get the damped interaction energy between and ,where and . Notice that in the limit , as expected, since as , and that in the limit , , , 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 we first apply the first righthand component of equation (3.5.2.16) together with the Fourier transforms of the normalized Gaussians and from equation (3.5.2.18) to obtainwhere we have used the fact that = , with the threedimensional delta function, and similarly for the integral over . Note that in the limit , , reduces to the equivalent expression for unit point charges at P and Q, and in the limit , reduces to the Fourierspace expression for the direct Coulomb interaction between the Gaussians and .
Recalling the developments leading to equation (3.5.2.28), and using equations (3.5.3.6) and (3.5.3.7), we can write the interaction energy from equation (3.5.3.1) aswhere , and .
Note that although the energy of a finiteexponent 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 (3.5.3.3) and (3.5.3.6), if p = q, then as , . We can now generalize equation (3.5.2.29). Let be normalized spherical Gaussian charge distributions with Gaussian exponents , centred at (point charges are included by taking ), and let satisfy (neutral unit cell). Then the energy of the central unit cell within a large spherical crystal, due to interactions of the Gaussian charge distributions with each other and with their periodic images within the crystal, is given bywhere , , and as before is the unitcell dipole.
Note that from the above derivation, the in is allowed to be different for each pair ij. One choice that leads to an efficient algorithm (particularly when combined with the fastFouriertransformbased methods discussed below) is to separate the Gaussian charge distributions into compact and diffuse sets, based on their exponents , and choose 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 so that is constant for all compact pairs ij. More specifically, given , a Gaussian is classified as compact if and diffuse otherwise. We denote these cases as and , respectively. Then for , choose so that = . Otherwise, choose to be infinite. Then the Coulomb energy of the central unit cell can be reexpressed as
Note that in the limit that for all i, all of the Gaussians become compact and this formula reduces to equation (3.5.2.30) for a large spherical crystal of point charges. The volumedependent term is somewhat surprising, and would not appear in a more straightforward Ewald derivation such as our first Ewald derivation. It vanishes by unitcell neutrality for the case of point charges or multipoles, and would also vanish if we chose a single for all Gaussian pairs.
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 be a normalized spherical Gaussian charge distribution with exponent , centred at , i.e. . The associated normalized Hermite Gaussians , where t, u, v are nonnegative integers, are defined byNote that this definition of Hermite Gaussians differs from that found in some textbooks, e.g. Helgaker et al. (2000), in that we have normalized the Gaussian 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. .
Above we discussed the interaction between spherical Gaussian charge distributions . Here we generalize these to higher angular momentum charge distributions given bywhere c is a vector of coefficients, a finite number of which are nonzero. For example, if is the only nonzero coefficient, we are reduced to the spherical charge distributions discussed above.
The fact that the partial derivatives within 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:
Next we consider the interaction energy between a generalized charge distribution centred at and another generalized charge distribution centred at , together with all the latter's periodic images, centred at , . Examining equation (3.5.3.8), we note that the terms to be differentiated all depend on . Thuswhere , i =1, 2 are given byandNote that by the arguments leading to equation (3.5.2.28), 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 of normalized Hermite Gaussians and a set of expansion points within U that are repeated periodically in the crystal (for example, some of the expansion points could be atomic nuclei). Let 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 be given bywhere, as above, for each i, l the coefficients are nonzero for a finite number of tuv. As above we can separate the exponents into compact and diffuse, i.e. or , according to whether or , respectively. Let denote the net charge at the expansion point .
We are interested in the Coulomb energy of charge densities , , interacting with each other and their periodic images at expansion points , for large K. As above, the direct Coulomb interaction of 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 densityfunctionaltheory approach, but clearly a nuclear charge cannot interact with itself. The unit cell is assumed to be neutral, i.e. . Then, using equations (3.5.3.9) and (3.5.3.10) we can write the Coulomb energy of = within the spherical crystal aswhere the structure factors are given by is given bywhere , and the unitcell dipole components are given by
In the limit that all Gaussians are compact with , reduces to the Coulomb energy of ideal multipoles positioned at the expansion points.
The Ewaldsum 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 speedups, 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.
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 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 pointcharge case or β for the compact–compact Gaussian or Hermite Gaussian interactions) so that directsum terms can be neglected past a fixed cutoff, i.e. when . Given a fixed Ewald parameter α, the cutoff may be fixed independently of the unitcell 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 unitcell volume or number of atoms in the unit cell U.
For the pointcharge 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). Ewald summation for ideal multipoles up to quadrupole have been discussed by Smith (1982, 1998) and Aguado & Madden (2003), using optimized combinations of terms to minimize the cost of the direct sum. A more general scheme, covering higherorder multipoles, was proposed by Sagui et al. (2004). In the latter scheme the McMurchie–Davidson recursion (McMurchie & Davidson, 1978; Helgaker et al., 2000) was used to calculate the higherorder derivatives of for the multipole–multipole interaction tensor. This same approach was used by Cisneros et al. (2006) 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 with , and let g(r) be a smooth function of r. For example, . 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 nonnegative integers t, u, v, let = . Then we have the following recursion:
The proof of the last recursion depends on the fact that and on Leibniz' generalization of the product rule for differentiation, namely , and is as follows:The other two recursions are proved similarly. To apply this to the case we recall from equation (3.5.3.4) that and that for the higherorder Boy's functions , . Efficient evaluation of the Boy's functions is discussed by Helgaker et al. (2000).
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) as well as in Cisneros et al. (2006).
The smooth particle mesh Ewald (PME) (Essmann et al., 1995) and fast Fourier Poisson (FFP) (York & Yang, 1994) methods both rely on the threedimensional 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 (Perram et al., 1988) to , where N is the number of particles and e.g. denotes a quantity such that remains bounded as . For the pointcharge case, the key insight is that if the point charges happen to be located on a regular grid, the structure factors [equation (3.5.2.14)] become discrete Fourier transforms and can be rapidly calculated using FFTs.
Like other variants of the original particle–particle particle–mesh (P^{3}M) method (Hockney & Eastwood, 1981), the smooth PME method interpolates charge onto grid positions. The smooth PME, like the original PME algorithm (Darden et al., 1993), 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 , and becomes progressively less accurate as grows. However, the weight function multiplying the structure factor in the reciprocal sum severely damps the highfrequency (large ) terms, thus preserving overall accuracy. In contrast, the FFP method relies on the fact that is the Fourier transform of a superposition of normalized spherical Gaussian charge distributions , centred at , 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 Xray 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 . 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, P^{3}M methods based on leastsquares approximation of the exact Ewald reciprocal sum gradient (Pollock & Glosli, 1996; Darden et al., 1997) can be made more efficient than the smooth PME method at the same accuracy, despite requiring more FFTs (Deserno & Holm, 1998). 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 Bsplines, does not need more FFTs as the order of differentiation (t + u + v) in equation (3.5.3.12) grows.
We begin by reexpressing the complex exponentials appearing in the structure factors in equations (3.5.2.14), (3.5.2.26) and (3.5.3.12). Recall the form of m from equation (3.5.2.2), and let be positive integers of the scale of the unitcell dimensions. Then for we can writewhere and , , and thus
To motivate the smooth PME approach, suppose we can approximate by some function of the formwhere is a well behaved function with continuous derivatives as needed for the general structure factor in equation (3.5.3.12). We first examine the structure factor for pointcharge Ewald sums, defined in equation (3.5.2.14). This can be approximated aswhere is given byand we have used the fact that = for .
Next we turn to the general structure factors from equation (3.5.3.12). Define , as above for expansion points , and note thatwhere are constants depending on the reciprocal bases and on , . Similar expansions hold for the other partial derivatives. Thus we can approximate the general structure factor bywhere are the transformed Hermite coefficients obtained from by change of variables, using and where is given by
In equation (3.5.4.2) we recognize the approximation to as times the threedimensional discrete Fourier transform of , which can be calculated rapidly by the 3DFFT. If the function B(w) has finite support then can be calculated rapidly [in order(N) time] as well. Similarly we have the 3DFFT of in equation (3.5.4.4). Thus the efficient calculation of structure factors depends on the approximation given in equation (3.5.4.1), which in turn rests on the choice of B(w). The smooth PME relies on Bsplines.
We proceed by deriving the facts needed about Bsplines. The Bsplines of order n, B_{n}(w) are functions of a real variable defined as follows:
We use the following properties of the Bsplines:
The first property to be established is clearly true for n = 1 and if true for n, then in the above integral recursively defining B_{n+1} if or the integrand must be zero for , so . 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 = = . Thus it is true in general by mathematical induction.
Next we show symmetry: . Again this is clearly true for n = 1. Assuming it is true for n, we havewhere we substituted in the second line above.
The fourth property is an efficient recursion for the derivatives of Bsplines. It can be seen directly for n = 2 for . Assume it is true for . Then for ,thus it is true in general by induction.
The final property is an efficient recursion for . It can be seen directly for n = 2. Assume it is true for . Then for ,where we have used property (4) 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 righthand side over to the left side of the equation.
Next we examine properties of the Bspline approximations to the complex exponentials that appear in the structure factors. Let be given byThe complex number is chosen to optimize the approximation of by as shown below. It is convenient to first study the properties of the function given bywhich is an approximation to the constant one. From its definition it is clear that is periodic in w with period one, i.e. . Therefore it can be expanded in a Fourier series:The Fourier coefficients are given bywhere denotes the Fourier transform of . Since are defined recursively by convolutions, its Fourier transform is obtained simply from that of : . Here we are using the fact that if f, g and h are related by convolution, , then . Since can be obtained by straightforward integration, we have
The original smooth PME method set , where was chosen to enforce interpolation: or for integer values of w. This was related to the socalled Euler spline, studied for more general complexvalued interpolation problems by Schoenberg (1973). The special case of interpolating the function is simpler, due to the Fourierseries representation, allowing us to to derive necessary properties in a selfcontained way. From equations (3.5.4.6) and the above expression for , noting that for integral z with , the interpolation condition can be written
Note that except for , the expressions appearing in the interpolation condition are periodic in z with period one. The interpolation condition thus forces to be periodic as well, and so we focus attention on the range . If the sum in the above interpolation condition is bounded, so in order to be able to choose to satisfy the interpolation condition it is sufficient to demonstrate that for . 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 , etc. Similarly the interpolation condition fails for z = −½. However, for , , etc. Similarly for . Thus we can say that for and the interpolation condition holds, and in this case we can write aswhereThe interpolating is obtained from this by multiplying by . Numerically, it is more convenient to use the alternative expression from equation (3.5.4.5), with given byNote that due to the above results for , if and .
Subsequent to the work by Essmann et al. (1995), we found (Darden et al., 1997) that a more accurate approximation was given by leastsquares fitting of by , or equivalently 1 by . Thus, with given by equation (3.5.4.7) we are looking for a complex number that minimizes the meansquare error MSE, given bywhere and denote the real and imaginary parts of . Since , it is clear that we can set to zero, that is, restrict ourselves to real . Then, minimizing over real we havewith given in equation (3.5.4.8).
Thus we arrive at the leastsquares optimal approximation to by setting in equation (3.5.4.5). The meansquare error due to this can be seen to beand thus the error in the approximation is seen to be small for small but increasing as , as discussed above (the reciprocal Ewald weight term compensates for the increased error by rapidly damping the larger terms). For the error decreases to zero as the Bspline order . The accuracy of the PME is thus increased by raising the order of interpolation and/or the grid size , which lowers for a given m. When derivatives with respect to w of are considered in this latter representation, their effect is to effectively lower the order of interpolation in the coefficients . Thus the order of interpolation needs to be higher, the higher the order of differentiation considered.
The fast Fourier Poisson (FFP) method (York & Yang, 1994) accelerates the calculation of the structure factors in equations (3.5.2.14) and (3.5.3.12) in much the same way as FFT methods accelerate structurefactor and densitymap 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 (3.5.2.15) as follows:where the modified structure factor is given by
We recognize as the threedimensional Fourier transform, evaluated at m, of a density given bywhere is a normalized Gaussian with Gaussian exponent 2α.
In turn we approximate this integral by a discrete Fourier transform. As above let be positive integers of the scale of the unitcell dimensions. Let , where V is the volume of the unit cell. Recall the lattice basis vectors , . Given integers define by . Then we can approximate the modified structure factor bywhere is given by
The FFP can be generalized to Hermite Gaussians by sampling that density in place of the superposition of spherical Gaussians. In equation (3.5.3.11) the termcan be handled by the methods just outlined. However the termcan present problems. The difficulty comes in the case when, for example, is compact but 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 FustiMolnar & Pulay (2002). Their solution was to filter out high frequencies in the Fouriertransformed compact densities and back transform them to obtain diffuse realspace densities. Instead, Cisneros et al. (2006) adjusted the exponents of both the compact and diffuse Gaussians. Choose a such that if is diffuse . For example, choose . As before, define by , and byThen, if and ,Thus instead of approximating by sampling the Gaussian followed by applying the 3DFFT, they approximated or . Note that with the above choice of , for any , and for any , . That is, all the modified Gaussians can be sampled on a uniform grid. Finally, note that for the diffuse–diffuse interactionsi.e. the interaction between adjusted diffuse Gaussians can be corrected efficiently in reciprocal space.
To test the efficiency of these reciprocalspace 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). As in that paper, the electron density of the water molecule was modelled by a fivesite 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 densityfunctional literature. These basis functions were obtained from the Environmental Molecular Sciences Laboratory at Pacific Northwest Laboratories (http://www.emsl.pnl.gov/forms/basisform.html ). The A2 basis set expansion uses 21 Hermite primitives with 102 Hermite coefficients per water molecule, resulting in a rootmeansquare relative molecular force error of about 8% compared to 6–31G* B3LYP intermolecular electrostatics (Cisneros et al., 2006). 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 rootmeansquare 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 3.5.4.1. Note that the times given here are considerably faster than those reported in Figs. 13 and 14 of Cisneros et al. (2006). 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 3.5.3.1, 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 directsum or reciprocalsum 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 N^{3/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.

The error in the FFP for a given Gaussian exponent is controlled by the size of the grid 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 Bspline in equation (3.5.4.3). The additional cost due to this has been ameliorated somewhat by the Gaussian split Ewald approach (Shan et al., 2005). Additionally, in rectangular unit cells fast Poisson solvers such as multigrid (Beck, 2000; Sagui & Darden, 2001) 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) 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.
References
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). Realspace 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. & HeadGordon, 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). Xray Charge Densities and Chemical Bonding. New York: Oxford University Press.
Coppens, P. & Volkov, A. (2004). The interplay between experiment and theory in chargedensity 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). Longrange 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. & ZicovichWilson, 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.
FustiMolnar, 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—Hy (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 ElectronicStructure Theory. Chichester: John Wiley and Sons.
Hockney, R. W. & Eastwood, J. W. (1981). Computer Simulation Using Particles. New York: McGrawHill.
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). Oneelectron and 2electron integrals over Cartesian Gaussian functions. J. Comput. Phys. 26, 218–231.
Mooji, W. T. M., van Duijneveldt, F. B., van Duijneveldtvan 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 twodimensional 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 1hydroxy7azabenzotriazole. 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 molecularinteractions 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, highprecision 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.