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.4, pp. 449457
https://doi.org/10.1107/97809553602060000771 Chapter 3.4. Accelerated convergence treatment of R^{−n} lattice sums^{a}Department of Chemistry, University of Louisville, Louisville, Kentucky 40292, USA Acceleratedconvergence treatment of the Coulombic lattice sum of ionic crystals and the dispersionenergy sum of ionic or molecular crystals is described. The Coulombic energy is very longrange and convergence of the Coulombic latticeenergy sum is extremely slow. Acceleratedconvergence treatment is indispensable to achieve accuracy with a reasonable amount of computational effort. The dispersionenergy sum has somewhat better convergence properties than the Coulombic sum, but acceleratedconvergence treatment of the dispersion sum is also strongly recommended since its use can yield at least an order of magnitude improvement in accuracy for a given calculation effort. 
The electrostatic energy of an ionic crystal is often represented by taking a pairwise sum between charge sites interacting via Coulomb's law (the sum). The individual terms may be positive or negative, depending on whether the pair of sites have charges of the same or different signs. The Coulombic energy is very longrange, and it is well known that convergence of the Coulombic latticeenergy sum is extremely slow. For simple structure types Madelung constants have been calculated which represent the Coulombic energy in terms of the cubic lattice constant or a nearestneighbour distance. Glasser & Zucker (1980) give tables of Madelung constants and review the subject giving references dating back to 1884. If the ionic crystal structure is not of a simple type usually no Madelung constant will be available and the Coulombic energy must be obtained for the specific crystal structure being considered. In carrying out this calculation, acceleratedconvergence treatment of the Coulombic lattice sum is indispensable to achieve accuracy with a reasonable amount of computational effort. A model of a molecular crystal may include partial net atomic charges or other charge sites such as lonepair electrons. The sum also applies between these site charges.
The dispersion energy of ionic or molecular crystals may be represented by an sum over atomic sites, with possible inclusion of terms for higher accuracy. The dispersionenergy sum has somewhat better convergence properties than the Coulombic sum. Nevertheless, acceleratedconvergence treatment of the dispersion sum is strongly recommended since its use can yield at least an order of magnitude improvement in accuracy for a given calculation effort. The repulsion energy between nonbonded atoms in a crystal may be represented by an exponential function of short range, or possibly by an function of short range. The convergence of the repulsion energy is fast and no acceleratedconvergence treatment is normally required.
This pairwise sum is taken between atoms (or sites) in the reference unit cell and all other atoms (or sites) in the crystal, excluding the self terms. Thus, the second atom (or site) is taken to range over the entire crystal, with elimination of selfenergy terms. If represents an energy, each atom is assigned one half of the pair energy. Therefore, the energy per unit cell iswhere is a given coefficient, is an interatomic distance, and the prime on the second sum indicates that self terms are omitted. In the case of the Coulombic sum, and is the product of the site charges.
Table 3.4.2.1 gives an example of the convergence behaviour of the untreated Coulombic sum for sodium chloride. Even at the rather large summation limit of 20 Å the Coulombic lattice sum has not converged and is incorrect by about 8%. The 20 Å sum included 832 molecules and 2494 individual distances. At various smaller summation limits the truncation error fluctuates wildly and can be either positive or negative. Note that the results shown in the table always refer to summation over whole molecules, that is, over neutral charge units.

If the Coulombic summation is not carried out over neutral charge units the truncation error is even larger. These considerations support the conclusion that acceleratedconvergence treatment of the Coulombic lattice sum should be regarded as mandatory. Table 3.4.2.2 gives an example of the convergence behaviour of the untreated dispersion sum for benzene. In obtaining this sum it is not necessary to consider whole molecules as in the Coulombic case. The exclusion of atoms (or sites) in the portions of molecules outside the summation limit greatly reduces the number of terms to be considered. At the summation limit of 20 Å, 439 benzene molecules and 22 049 individual distances are considered; the dispersionsum truncation error is 0.4%. Thus, if sufficient computer time is available it may be possible to obtain a moderately accurate dispersion sum without the use of accelerated convergence. However, as shown below, the use of accelerated convergence will greatly speed up the calculation, and is in practice necessary if higher accuracy is required.

Ewald (1921) developed a method which modified the mathematical representation of the Coulombic lattice sum to improve the rate of convergence. This method was based on partially transforming the lattice sum into reciprocal space. Bertaut (1952) presented another method for derivation of the Ewald result which used the concept of the crystallographic structure factor. His formula extended the Ewald treatment to a composite lattice with more than one atom per lattice point. Nijboer & DeWette (1957) developed a general Fourier transform method for the evaluation of sums in simple lattices. Williams (1971) extended this treatment to a composite lattice and gave general formulae for the sums for any crystal. A review article, on which this chapter is based, appeared later (Williams, 1989a,b).
Consider a function, W(R), which is unity at and smoothly declines to zero as R approaches infinity. If each term of the lattice sum is multiplied by W(R), the rate of convergence is increased. However, the rate of convergence of the remainder of the original sum, which contains the difference terms, is not increased.
In the acceleratedconvergence method the difference terms are expressed as an integral of the product of two functions. According to Parseval's theorem (described below) this integral is equal to an integral of the product of the two Fourier transforms of the functions. Finally, the integral over the Fourier transforms of the functions is converted to a sum in reciprocal (or Fouriertransform) space. The choice of the convergence function W(R) is not unique; an obvious requirement is that the relevant Fourier transforms must exist and have correct limiting behaviour. Nijboer and DeWette suggested using the incomplete gamma function for W(R). More recently, Fortuin (1977) showed that this choice of convergence function leads to optimal convergence of the sums in both direct and reciprocal space:where and are the gamma function and the incomplete gamma function, respectively:and
3.4.4. Preliminary derivation to obtain a formula which accelerates the convergence of an R^{−n} sum over lattice points X(d)
The threedimensional directspace crystal lattice is specified by the origin vectors , and . A general vector in direct space is defined aswhere are the fractional cell coordinates of X. A lattice vector in direct space is defined aswhere are integers (specifying particular values of ) designating a lattice point. is the directcell volume which is equal to . A general point in the direct lattice is X(x); the contents of the lattice are by definition identical as the components of x are increased or decreased by integer amounts.
The reciprocallattice vectors are defined by the relationsA general vector in reciprocal space H(r) is defined asA reciprocallattice vector H(h) is defined by the integer triplet (specifying particular values of ) so thatIn other sections of this volume a shortened notation h is used for the reciprocallattice vector. In this section the symbol H(h) is used to indicate that it is a particular value of H(r).
The threedimensional Fourier transform of a function is defined byThe Fourier transform of the set of points defining the direct lattice is the set of points defining the reciprocal lattice, scaled by the directcell volume. It is useful for our purpose to express the lattice transform in terms of the Dirac delta function which is defined so that for any function We then writeFirst consider the lattice sum over the directlattice points X(d), relative to a particular point , with omission of the origin lattice point.The special case with will also be needed:Now define a sum of Dirac delta functionsThen S′ can be represented as an integralin which a term is contributed to S′ whenever the directspace vector X coincides with the lattice vector X(d), except for . Now apply the convergence function to S′:
The first integral is shown here only for the purpose of giving a consistent representation of S′; in fact, the first integral will be reconverted back into a sum and evaluated in direct space. The second integral will be transformed to reciprocal space using Parseval's theorem [see, for example, Arfken (1970)], which states thatConsidering only the second integral in the formula for S′ and explicitly introducing the term we havewhere the unprimed f includes the term which was earlier omitted from f′:Using Parseval's theorem, and evaluating the origin term, we haveThe Fourier transform of the complement of the incomplete gamma function divided by is (Nijboer & DeWette, 1957)If there is a change of origin and the point is used instead of X the transform isEvaluation of the two Fourier transforms in the first term givesBecause of the presence of the Dirac delta function in each integral, we can convert the integrals with h unequal to zero into a sumThe term needs to be evaluated in the limit. Clearly, the complex exponential goes to unity. If n is greater than 3 the limit of the indeterminate form infinity/infinity is needed:The limit can be found by L'Hospital's rule [see, for example, Widder (1961)] which states that if and both approach infinity as x approaches a constant, c, and the limit of the ratio of the first derivatives and exists, that limit is also true for the limit of the ratio of the functions:To differentiate the definite integral function, Leibnitz's formula may be used [see, for example, Arfken (1970)]. This formula states thatIn our case, x becomes ; f becomes which is independent of ; g becomes ; and h is infinite. Thus only the last term of Leibnitz's formula is nonzero and we obtain for the ratio of the first derivativesso that the limiting value for the term for n greater than 3 isThe final result for S′ is
The significance of the terms is as follows. The first term represents the convergenceaccelerated direct sum, which does not include the origin term; the next term, also in direct space, corrects for the remainder resulting from the subtraction of the origin term; the third term comes from Parseval's theorem and is a sum over the nonzero h reciprocallattice points; and the last term is the reciprocallattice term.
If the second term becomes an indeterminate form 0/0. The limit can be found with use of L'Hospital's rule again, this time for the 0/0 form. We need the limit of , where and . To differentiate the incomplete gamma function, we can again use Leibnitz's formula. In this case only the second term of the formula is nonzero and we obtain for the ratio of the first derivativesso that the limiting value for this term as approaches zero isTherefore, the value of the sum when is
Define a general lattice sum over directspace points which interact with pairwise coefficients , where :where the prime indicates that when the selfterms with are omitted. For convenience the terms may be divided into three groups: the first group of terms has , where j is unequal to k; the second group has d not zero and j not equal to k; and the third group had d not zero and . (A possible fourth group with and is omitted, as defined.)By expanding this expression we obtainThis expression for V has nine terms, which are numbered on the righthand side. Term (3) can be expressed in terms of Γ rather than γ:It is seen that cancellation occurs with term (1) so thatwhich is the , j unequal to k portion of the treated directlattice sum. The d unequal to 0, j unequal to k portion corresponds to term (2) and the d unequal to 0, portion corresponds to term (6). The directlattice terms may be consolidated asNow let us combine terms (4) and (8), carrying out the h summation first:Terms (5) and (9) may be combined:The final formula is shown below. The significance of the four terms is: (1) the treated directlattice sum; (2) a correction for the difference resulting from the removal of the origin term in direct space; (3) the reciprocallattice sum, except ; and (4) the term of the reciprocallattice sum.
As taken above, the limit of the reciprocallattice term of or existed only if n was greater than 3. The corresponding contributions to were terms (5) and (9) of Section 3.4.5. To extend the method to we will show in this section that these terms vanish if conditions of unitcell neutrality and zero dipole moment are satisfied.
The integral representation of the term (5) isand for term (9) isCombining these two sums of integrals into one integral sum gives
For , suppose are net atomic charges so that the geometric combining law holds for . Then the double sum over j and k can be factored so that the limit that needs to be considered isIf the unit cell does not have a net charge, the sum over the q's goes to zero in the limit and this is a 0/0 indeterminate form. Let approach zero along the polar axis so that , where subscript 3 indicates components along the polar axis. To find the limit with L'Hospital's rule the numerator and denominator are differentiated twice with respect to . Represent the numerator of the limit by the product and note thatIt is seen that in addition to cell neutrality the product of the first derivatives of the sums must exist. These sums areandwhich vanish if the unit cell has no dipole moment in the polar direction, that is, if . Since the second derivative of the denominator is a constant, the desired limit is zero under the specified conditions. Now the polar direction can be chosen arbitrarily, so the unit cell must not have a dipole moment in any direction for the limit of the numerator to be zero. Thus we have the formula for the Coulombic lattice sumwhich holds on conditions that the unit cell be electrically neutral and have no dipole moment. If the unit cell has a dipole moment, the limiting value discussed above depends on the direction of H. For methods of obtaining the Coulombic lattice sum where the unit cell does have a dipole moment, the reader is referred to the literature (DeWette & Schacher, 1964; Cummins et al., 1976; Bertaut, 1978; Massidda, 1978).
If the denominator considered for the limit in the preceding section is linear in H so that only one differentiation is needed to obtain the limit by L'Hospital's method. Since a term of the type is always a factor, the requirement that the unit cell have no dipole moment can be relaxed. For the zerocharge condition is still required: . When the expression becomes determinate and no differentiation is required to obtain a limit. In addition, factoring the sums into sums is not necessary so that the only remaining requirement for this term to be zero is , which is a further relaxation beyond the requirement of cell neutrality.
The structure factor with generalized coefficients is defined byThe corresponding Patterson function is defined byThe physical interpretation of the Patterson function is that it is nonzero only at the intersite vector points . If the origin point is removed, the lattice sum may be expressed as an integral over the Patterson function. This origin point in the Patterson function corresponds to intersite vectors with and :Using the incomplete gamma function as a convergence function, this formula expands into two integralsThe first integral is shown only for a consistent representation; actually it will be reconverted to a sum and evaluated in direct space. The first part of the second integral will be evaluated with Parseval's theorem and the second part in the limit as approaches zero:The first Fourier transform (of the Patterson function) is the set of amplitudes of the structure factors and the second Fourier transform has already been discussed above; the method for obtaining the limit (for n equal to or greater than 1) was also discussed above. The result obtained isThe integral can be converted into a sum, since is nonzero only at the reciprocallattice points:The term with is evaluated in the limit, for n greater than 3, asSince , this term is identical with the third term of as derived earlier. The case of is handled in the same way as previously discussed, where the limit of this term is zero provided the unit cell has no net charge or dipole moment.
The incomplete gamma function may be expressed in terms of commonly available functions such as the exponential integral and the complement of the error function. The definition of the exponential integral isThe definition of the complement of the error function isNumerical approximations to these functions are given, for example, by Hastings (1955). The recursion formula for the incomplete gamma function (Davis, 1972)may be used to obtain working formulae starting from the special values of and which are defined above. Also we note that .
Let us consider the case where the unit cell contains Z molecules which are related by Z symmetry operations, and it is desired to include only intermolecular distances in the summation. In the direct sum (1) the indices j and k will then run only over the asymmetric unit, and all terms with are eliminated. The calculated energy refers then to one molecule (or mole) rather than to one unit cell. The correction term (2) also refers to one molecule according to the range of j and k. Since the reciprocallattice sum refers to the entire unit cell, terms (3) and (4) need to be divided by Z to refer the energy to one molecule.
Both the direct and reciprocal sums must be corrected for the elimination of intramolecular terms. Using the convergence function , we haveAs mentioned above, the second summation term, which is the intramolecular term in direct space, is simply left out of the calculation. When using the acceleratedconvergence method the third and fourth summation terms are always obtained, evaluated in reciprocal space. The undesired inclusion of the intramolecular term (fourth term above) in the reciprocalspace sum may be compensated for by explicit subtraction of this term from the sum.
In this section let and . Let ; . If the geometric mean combining law holds, ; letThenwhereandThe formulae below describe in terms of and ; the distance is simply represented by .
Consider the case of the sodium chloride crystal structure (a facecentred cubic structure) as a simple example for evaluation of the Coulombic sum. The sodium ion can be taken at the origin, and the chloride ion halfway along an edge of the unit cell. The results can easily be generalized for this structure type by using the unitcell edge length, a, as a scaling constant.
First, consider the nearest neighbours. Each sodium and each chloride ion is surrounded by six ions of opposite sign at a distance of . The Coulombic energy for the first coordination sphere is . Table 3.4.2.1 shows that the converged value of the lattice energy is . Thus the nearestneighbour energy is over three times more negative than the total lattice energy. In the second coordination sphere each ion is surrounded by 12 similar ions at a distance of . The energy contribution of the second sphere is . Thus, major cancellation occurs and the net energy for the first two coordination spheres is which actually has the wrong sign for a stable crystal. The third coordination sphere again makes a negative contribution. Each ion is surrounded by eight ions of opposite sign at a distance of . The energy contribution is , now giving a total so far of . In the fourth coordination sphere each ion is surrounded by six others of the same sign at a distance of a. The energy contribution is +(1/2)(12)(1/a)(1389.3654) = +8336.19/a to yield a total of .
It is seen immediately by examining the numbers that the Coulombic sum is converging very slowly in direct space. Madelung (1918) devised a method for accurate evaluation of the sodium chloride lattice sum. However, his method is not generally applicable for more complex lattice structures. Evjen (1932) emphasized the importance of summing over a neutral domain, and replaced the sum with an integral outside of the first few shells of nearest neighbours. But the method of Ewald remained as the only completely general and accurate method of evaluating the Coulombic sum for a general lattice. Although it was derived in a somewhat different way, Ewald's method is equivalent to accelerated convergence for the special case of .
In the reciprocal lattice of sodium chloride only points with indices (hkl) all even or all odd are permitted by the facecentred symmetry. The reciprocal cell has edge length and the reciprocalaxis directions coincide with the directlattice axis directions. The closest points to the origin are the eight (111) forms at a distance of . For sodium chloride this distance is 0.3078 Å^{−1}.
Table 3.4.12.1 shows the effect of convergence acceleration on the directspace portion of the sum for the sodium chloride structure. The constant term is included in the values given. This constant term is always large if w is not zero; for instance, when this term is (Table 3.4.12.2). For the reciprocallattice sum is zero to six figures. Thus, only the direct sum (plus the constant term) is needed, evaluated out to 20 Å in direct space, to obtain sixfigure accuracy. As shown in Table 3.4.2.1 above, the same summation effort without the use of accelerated convergence gave 8% error, or only slightly better than onefigure accuracy. The acceleratedconvergence technique therefore yielded nearly five orders of magnitude improvement in accuracy, even without evaluation of the reciprocallattice sum.


The column showing shows an example of how the reciprocallattice sum can also be neglected if lower accuracy is required. Table 3.4.12.2 shows that the reciprocallattice sum is still only 0.003. But now the directlattice sum only needs to be evaluated out to 14 Å, with further savings in calculation effort. For w values larger than 0.15 the reciprocal sum is needed. For this sum must be evaluated out to 0.8 Å^{−1} to obtain sixfigure accuracy.
Table 3.4.12.3 illustrates an application for the dispersion sum. When five figures of accuracy can be obtained without consideration of the reciprocal sum. The direct sum is required out to 18 Å. If , better than fourfigure accuracy can still be obtained without evaluating the reciprocallattice sum. In this case, the direct lattice needs to be summed only to 12 Å, and there is a saving of an order of magnitude in the length of the calculation. As with the Coulombic sum, if w is greater than 0.15 the reciprocallattice summation is needed; Table 3.4.12.4 shows the values.


The time required to obtain a lattice sum of given accuracy will vary depending on the particular structure considered and of course on the computer and program which are used. An example of timing for the benzene dispersion sum is given in Table 3.4.12.5 for the PCK83 program (Williams, 1984) running on a VAX11/750 computer. In this particular case direct terms were evaluated at a rate of about 200 terms s^{−1} and reciprocal terms, being a sum themselves, were evaluated at a slower rate of about 5 terms s^{−1}.

Table 3.4.12.5 shows the time required for evaluation of the dispersion sum using various values of the convergence constant, w. The timing figures show that there is an optimum choice for w; for the PCK83 program the optimum value indicated is 0.15–0.2 Å^{−1}. In the program of Pietila & Rasmussen (1984) values in the range 0.15–0.2 Å^{−1} are also suggested. For the WMIN program (Busing, 1981) a slightly higher value of 0.25 Å^{−1} is suggested. Trial calculations can be used to determine the optimum value of w for the situation of a particular crystal structure, program and computer.
References
Arfken, G. (1970). Mathematical Methods for Physicists, 2nd ed. New York: Academic Press.Bertaut, E. F. (1952). L'énergie électrostatique de réseaux ioniques. J. Phys. (Paris), 13, 499–505.
Bertaut, E. F. (1978). The equivalent charge concept and its application to the electrostatic energy of charges and multipoles. J. Phys. (Paris), 39, 1331–1348.
Busing, W. R. (1981). WMIN, a computer program to model molecules and crystals in terms of potential energy functions. Oak Ridge National Laboratory Report ORNL5747. Oak Ridge, Tennessee 37830, USA.
Cummins, P. G., Dunmur, D. A., Munn, R. W. & Newham, R. J. (1976). Applications of the Ewald method. I. Calculation of multipole lattice sums. Acta Cryst. A32, 847–853.
Davis, P. J. (1972). Gamma function and related functions. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, edited by M. Abramowitz & I. A. Stegun, pp. 260–262. London, New York: John Wiley. [Reprint, with corrections of 1964 Natl Bur. Stand. publication.]
DeWette, F. W. & Schacher, G. E. (1964). Internal field in general dipole lattices. Phys. Rev. 137, A78–A91.
Evjen, H. M. (1932). The stability of certain heteropolar crystals. Phys. Rev. 39, 675–694.
Ewald, P. P. (1921). Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. (Leipzig), 64, 253–287.
Fortuin, C. M. (1977). Note on the calculation of electrostatic lattice potentials. Physica (Utrecht), 86A, 574–586.
Glasser, M. L. & Zucker, I. J. (1980). Lattice sums. Theor. Chem. Adv. Perspect. 5, 67–139.
Hastings, C. Jr (1955). Approximations for digital computers. New Jersey: Princeton University Press.
Madelung, E. (1918). Das elektrische Feld in Systemen von regelmässig angeordneten Punktladungen. Phys. Z. 19, 524–532.
Massidda, V. (1978). Electrostatic energy in ionic crystals by the planewise summation method. Physica (Utrecht), 95B, 317–334.
Nijboer, B. R. A. & DeWette, F. W. (1957). On the calculation of lattice sums. Physica (Utrecht), 23, 309–321.
Pietila, L.O. & Rasmussen, K. (1984). A program for calculation of crystal conformations of flexible molecules using convergence acceleration. J. Comput. Chem. 5, 252–260.
Widder, D. V. (1961). Advanced Calculus, 2nd ed. New York: PrenticeHall.
Williams, D. E. (1971). Accelerated convergence of crystal lattice potential sums. Acta Cryst. A27, 452–455.
Williams, D. E. (1984). PCK83, a crystal molecular packing analysis program. Quantum Chemistry Program Exchange, Department of Chemistry, Indiana University, Bloomington, Indiana 47405, USA.
Williams, D. E. (1989a). Accelerated convergence treatment of R^{−n} lattice sums. Crystallogr. Rev. 2, 3–23.
Williams, D. E. (1989b). Accelerated convergence treatment of R^{−n} lattice sums. Corrections. Crystallogr. Rev. 2, 163–166.