A fully analytical integration of properties over the 3D volume of the β sphere in topological atoms.

Atomic multipole moments associated with a spherical volume fully residing within a topological atom (i.e., the β sphere) can be obtained analytically. Such an integration is thus free of quadrature grids. A general formula for an arbitrary rank spherical harmonic multipole moment is derived, for an electron density comprising Gaussian primitives of arbitrary angular momentum. The closed expressions derived here are also sufficient to calculate the electrostatic potential, the two types of kinetic energy, as well as the potential energy between atoms. Some integrals have not been solved explicitly before but through recursion and substitution are broken down to more elementary listed integrals. The proposed method is based on a central formula that shifts Gaussian primitives from one center to another, which can be derived from the well-known plane-wave expansion (or Rayleigh equation).

Introduction

Quantum chemical topology (QCT),1 pioneered2 by the research group of the late Richard Bader, has carved out a space among non‐topological methods in extracting insight from wave functions, both in chemistry (including metal‐metal interactions3) and solid state physics (including high‐resolution X‐ray crystallography4). This imaginative approach started, in quantum chemistry, as (the Quantum Theory of) Atoms in Molecules5 (QT)(AIM) in the early 1970s. After adopting the mathematical language of dynamical systems with an eye on analyzing the electron density, QTAIM proposed a parameter‐free definition of an atom and a chemical bond (the latter continuing to be a topic of contemporary research). The concepts of dynamical systems (e.g., separatrix, critical point, basin, gradient path) have been applied with success to quantum mechanical density functions, other than the electron density, such as the Laplacian of the electron density, the nuclear potential, the electron localization function, the electron localizability indicator, the virial field, to name a few (references and more details see both the Introduction and Box 8.1 in Ref. 6). A detailed history of QTAIM and its QCT follow up has been given elsewhere.6, 7, 8

Although sometimes still mistaken9 for yet another population analysis, QTAIM offers a good number of important atomic properties other than an atomic charge, which all derive from a single universal formula. This formula is the 3D integration of a relevant integrand over the volume of the topological atom, and can yield atomic (kinetic) energy, volume, electrostatic potential, and multipole moments, to name a few common ones. The existence of a single formula that delivers atomic properties in a consistent way should at least be emphasized and perhaps even celebrated. Indeed, this is not true for other methods: for example, there is no Mulliken volume or Hirshfeld volume; similarly, there is no van der Waals charge (as opposed to volume or radius). Better controlling the atomic integration over a spherical volume that resides completely within a topological atom (see below) is a sign of progress.

An important development within the QCT framework is that of an energy partitioning method called interacting quantum atoms (IQA),10 which is based on earlier work.11 In spite of its great computational cost, IQA is increasingly applied to a wide variety of chemical phenomena, non‐exhaustively ranging from hydrogen bond cooperativity12 over metal carbonyl bonds,13 halogen bonds,14 Zn‐complexes,15 excited states,16 and congested molecules,17 to torsional energy barriers in peptides.18 IQA hereby avoids the pitfalls,19 both conceptual and numerical, of older non‐topological energy partitioning schemes such as EDA and SAPT. The surge of IQA has been made possible by improved algorithms (e.g., parallelization and parameter fine‐tuning) such as those found in the program20 AIMAll.

Still, the computational challenge of the topological partitioning remains, compared to cheaper non‐topological methods. This is why research that improves algorithms to integrate topological atoms should continue. This article focuses on the elimination of integration quadrature, but only inside the so‐called β sphere,21 which is the largest sphere that is completely contained within the topological atom at whose nucleus it is centered. In practice, the β sphere is just a sphere with an adequately large radius, typically ∼90% of the distance between the nucleus and the nearest bond critical point. Numerical volume integration over the whole topological atom does not require a β sphere, in principle. Indeed, the radial quadrature can be performed by a single interval stretching from the nucleus to the edge of the atom. However, the very high values of electron density near the nucleus make it numerically advantageous to keep the β sphere, and to work with two radial integration intervals, one inside the β sphere and one outside, each with its own quadrature grid. Heavy elements (fourth period and beyond) benefit from even more than two radial integration intervals. In 2011, a fully analytical 3D volume integration over the β‐sphere was achieved for the first time.22 Here, we propose an alternative analytical integration, which not only avoids the tedious axes rotations occurring in the former integration procedure but also generalizes much better to arbitrary multipole moments, and enables controlled elimination of Gaussian primitive functions far from the nucleus of the atom being integrated.

The Introduction to the 2011 paper contains a very detailed and rather exhaustive history at the time of topological integration algorithms, which will not be repeated here. Instead, we expeditiously review the literature, now extending it with the period 2012–2017, and starting with the first integrations over topological atoms, carried out23 already in 1973 but for linear molecules only. In 1981, the first computer program called OMEGA21 enabled the integration of topological atoms in any polyatomic molecule, followed by a completely different program24 a year later, called PROAIM. The Bader group gathered these programs, together with critical point location software, into the AIMPAC suite. With the help of collaborators, the main author of the atomic integration code, mathematician Biegler‐König, resurrected and perfected25, 26, 27 his effort of the 1980s, leading to the AIM2000 code. Meanwhile, in the 1990s, the program MORPHY28 appeared as the first integration code29 written independently from the AIMPAC suite. MORPHY, also for the first time combined all AIMPAC functionality into a single program. In 1997, AIMAll emerged from an extensive AIMPAC modification, and is now a popular and fast QCT code. Independent code, also developed30 in the 1990s, was implemented in the program GAUSSIAN but later withdrawn from it. In 1997, a research group in Oviedo (Spain, USE) created an independently written topological code called CRITIC,31 which is devoted to solid state systems, and later honed it to CRITIC2.32 In the early 2000s, they also wrote PROMOLDEN, a code specialized in molecular rather than solid state electron densities, and focusing on IQA. The QTREE method was also invented33 in Oviedo, and was based on the OCTREE algorithm34 proposed eight years earlier. The Electron Localization Function, a popular QCT tool35 to extract chemical information from wave functions, triggered the independent development of the topological code TopMoD.36 Gatti's TOPOND code37 determines the topology of an analytically expressed Hartree‐Fock or DFT wavefunction obtained by the program CRYSTAL.38 The high‐resolution X‐ray crystallography community also created its own algorithms early on such as that of TOPXD,39 NEWPROP,40 or WINXPRO.41 However, a topological analysis can also be carried out using grids,42, 43 such as the InteGriTy package.44 Grid methods are also applied in the context of computed electron densities.45, 46, 47 Two more algorithms are based on different principles, such as the “elastic sheet” method48 and a finite element method.49 In contrast, recent developments50, 51 turned back to the original triangulation method of Biegler‐König. Meanwhile, in the 2010s, the Chinese program Multiwfn52 (multifunctional wavefunction analyzer) emerged from a fast catch‐up exercise consisting of rewriting a large collection of pre‐existing algorithms developed over years by other groups. Multiwfn has an impressive functionality, encompassing both topological and non‐topological methods, but encourages crude parameter settings, and typically implements the simplest (i.e., mathematically elementary) algorithms.

In this article, we focus on the involved mathematics of deriving closed formulae for the β sphere's contribution to atomic multipole moments of arbitrary spherical harmonic rank, the two types of kinetic energy (K and G), the electrostatic potential and the potential energy between topological atoms. To keep the Method Section to the point, much material has been siphoned off to the Supporting Information. Due to the mathematical complexity of the analytical approach, the strategy and all details are reported in this article, while an implementation will follow in a future software‐oriented publication.

Method

Imagine a sphere with a given radius β, centered at a given nuclear position. In the context of topological atoms, this sphere is called the β sphere. This is the largest sphere that completely fits within a topological atom. Whichever the shape of this atom's interatomic surfaces, β must be smaller than the distance between the nucleus and the closest‐by bond critical point.

The first aim is to calculate, by analytical integration, the atomic multipole moments Q associated with the volume of the β sphere only. The integrand of the integral leading to Q(β) is the product of the electron density ρ and a regular spherical harmonic R centered on the center of atomic integration with position vector R , where the radius β acts as a parameter and where the regular spherical (or solid) harmonic R is defined by eq. (2),

Here r, θ, and φ represent spherical coordinates centered on R , which typically is the nucleus of (topological) atom Ω and Y represents a spherical harmonic, defined using the so‐called “quantum mechanics convention” and including the Condon‐Shortley factor.

Three comments are in order here. First, the uppercase indices L and M should not be confused with the lowercase indices ℓ and m, which will appear in formulae below. The regular spherical harmonics share their mathematical shape with the familiar s, p, d, f orbitals, or orbitals of higher rank. Second, eq. (2) defines the R functions as complex (for details see Part A in the Supporting Information). It is best to work with complex R functions during a derivation and, if necessary, convert them to real functions only at the end. Third, there is a local frame centered on a given nucleus, which determines the orientation of R in space. This local frame is parallel to the global frame, which renders the local and global frames identical but for a translation. Earlier work22 shows that one can make the nuclear position of the β sphere coincide with the origin of the global frame, without loss of generality.

The electron density appearing in eq. (1) is expanded in Gaussian primitives (denoted G(r – R)) of arbitrary angular momentum and centered at position R, where n MO signifies the number of molecular orbitals, n G is the number of Gaussian primitives used to expand each MO in, and c are the LCAO coefficients, which include the normalization factors of the Gaussian primitives. This notation makes clear that a Gaussian primitive is centered on either nuclear position R or R . Note that R may or may not coincide with the center of atomic integration R . The same is true for R , independently of whether R coincides with R .

A brief digression points out that, in a quadrature scheme, it is possible to evaluate a molecular orbital on its own, as part of the integration procedure calculating the integral of eq. (1). However, in an analytical integration, one cannot do this and one must instead consider the whole electron density in the integrand in eq. (1). This point is elaborated in Appendix A of the Supporting Information of Ref. 22.

Using the Gaussian product theorem, the product of any two Gaussians centered on R and R can be written as a single Gaussian, centered on R . Ignoring the angular part of the Gaussian primitives, the product theorem is formally written as

The new Gaussian exponent α, the new center R and prefactor K are also easily calculated,

Substitution of eq. (3) into eq. (1) leads to eq. (6),

The integration problem is then confined to the contribution from a single Gaussian to a given multipole moment and denoted Q(β), which is centered on a site that is, in principle, different to the site R at which the multipole moment is centered,

Without loss of generality one can position the atomic integration center R at the global origin, which is equivalent to translate the whole relative configuration of vectors r, R , R , and R to the global origin.

Figure 1 illustrates all vectors introduced so far, and how they relate to one another geometrically, for an arbitrary pair of Gaussian primitives centered on R and R , and contributing to the atomic integral calculated in eq. (7). Figure 1 illustrates the case of two Gaussian p‐type primitives, p and p , making a contribution to ρ(r) and hence to Q.

Figure 1

Schematic representation of the vectors involved in eqs. (3), (6), and (7). The atom to be integrated is centered at the origin o of the global frame, which provides coordinates for all vectors involved. Two arbitrary Gaussian primitives, a p and a p function, are respectively centered at R j and R k. The product of these Gaussian primitives is centered at R jk. The position vector r describes the electron density contributing to the atomic (volume) integral. The angle γ is pivotal in the separation of the variables. [Color figure can be viewed at wileyonlinelibrary.com]

For the general product of two Gaussian primitives of arbitrary angular momentum, eq. (7) becomes where the indices l, j, and k with superscript bars mark the angular momenta. One can use the binomial expression to work out each of the three products in eq. (8), one product in x, one in y and one in z. For example, focusing on x only, one obtains eq. (9), where the function f provides the coefficients corresponding to the power λ = r + s in x, which parametrically depend on the position of the centers (X, X) of the Gaussian primitives and their angular momenta (the various powers of X and X). We can then deduce from eq. (9) the two lowest powers of λ (0 and 1), while Part B of the Supporting Information explicitly lists f 2 and f 3 . The functions f are polynomials in the components of the position vectors R and R , which are expressed with respect to the global axis system (which is translated to the center of integration, for each atom in turn, as discussed above).

Equation (8) can then be rearranged as where λ is the integer power of the coordinate x, μ that of y, and ν that of z. Part B of the Supporting Information also applies eqs. (11) and (10) to the specific case illustrated in Figure 1, which shows two different p‐type Gaussian primitives (p and p) contributing to Q , .

To work out eq. (11) analytically, we recognize that R and the product Gaussian G have a different origin: the former is centered at the global origin o while the latter is centered at R . The integral of eq. (11) is naturally carried out in spherical coordinates centered at o because the electron density is essentially spherical, even for an atom inside a molecule. In other words, the deviations from sphericity are too small to abandon spherical coordinates in favor of alternative but inappropriate coordinates, such as Cartesian coordinates. The challenge is therefore to re‐express G in terms of spherical coordinates centered at o. This can be achieved using a key equation proposed before,53 where i ℓ(z) is the modified spherical Bessel function of integer order ℓ, and r 1 and r 2 are the respective magnitudes of vector r 1 and r 2, while Pℓ is the Legendre function of order ℓ and γ is the angle between these two vectors. Six remarks are in order here. First, eq. (12) can be derived from the well‐known plane‐wave expansion (or Rayleigh equation) as demonstrated in Part C of the Supporting Information. Second, no convergence condition accompanying this series expansion is reported. Third, the sum over ℓ needs to be truncated to a finite integer, which offers an opportunity to make the integral evaluation (see below) more efficient. Fourth, on application in the current context, which vector is r 1 and which is r 2 is not important because eq. (12) remains valid if these vectors are swapped. Fifth, the expansion of eq. (12) is pointless if either r 1 or r 2 vanishes. In that case, we trivially recover that 1 = 1. Indeed, the only surviving term in the sum is that of ℓ = 0, and i 0(0) = 1 and P 0(x) = 1 (even when γ is indeterminate, as is the case here, with an undefined angle γ resulting between a null vector and non‐null vector). Sixth and finally, the argument of an exponential function must be dimensionless and hence, if r has the dimension of length then r must have the dimension of reciprocal length, which is the case if it were a wave factor (as indeed it is in the Rayleigh equation).

To make use of eq. (12), the exponential function's argument in eq. (11) needs to be rewritten. This is done in eq. (13), such that eq. (12) applies to the second factor on the right‐hand side of eq. (13), setting r 1 = 2α R and r 2 = r where |R | = R and |r| = r. In terms of dimensionality (see 6th remark just above), the dimension of α is reciprocal length squared, and hence, the argument of the exponential in eq. (13) is again dimensionless. Substitution of eqs. (12) and (13) into eq. (11), and dropping the by now tedious indices j and k, one obtains eq. (14),

Two remarks must be made. If either R = 0 (center of product Gaussian coincides with the center of integration) or r = 0 (lower boundary of volume integration), then the expansion of eq. (12) becomes trivial but it still holds. From the latter statement, it follows that there is no need to treat the case of R = 0 separately, at least at this stage. Later in the derivation a further comment on R = 0 will be needed. In summary, although the case of R = 0 can be solved directly using spherical coordinates centered on the global origin in eq. (11) we proceed anyway with the expansion that is eq. (12). Second, the quantities R, K, α, and f in eq. (14) are parameters, depending on the fixed characteristics of the product Gaussian. The volume integral in eq. (14) is best described via spherical coordinates, denoted r, θ, and φ, and centered at the global origin o (see Fig. 1). The quantity γ is the angle between r and R and hence is a mixture of a parameter (R) and a variable (r). We aim to separate all variables from parameters, so the P l(cos γ) factor needs to be expanded. The well‐known addition theorem involving two spherical harmonics, given in eq. (15), achieves this aim, where vectors r and R are, respectively, described by spherical coordinates (r, θ, φ) and (R, θ R, φ R), which are both centered on o. The choice as to which variables describe the complex conjugate spherical harmonic is steered by the application of an advantageous orthonormalization relationship discussed below, the special case where λ = ν = μ = 0. Substituting eq. (15) into eq. (14), and rearranging delivers

Substituting eq. (2) into eq. (16) and making integration variables explicit yields,

Further rearrangement yields the following master equation, eq. (18),

In the special case where λ = ν = μ = 0, the angular integration in eq. (18) has been prepared for the use of the orthonormalization relationship announced above, such that eq. (18) becomes

This pleasing result shows that there is no need to evaluate any angular integrals in the case of only s‐like Gaussian primitives contributing to Q(β). Second, the multipole moment Q(β) is simply proportional to the value of Y(θ R, ϕ R). Note that this multipole moment Q(β) can be a complex number because Y can be complex and the other factors in eq. (20) are real. Third, the contribution of a given Gaussian primitive drops off quickly with the distance of its center from the origin o. The radial integral in eq. (20), at any order L, can be reduced to a linear combination of the basic integral of an exponential with quadratic argument. An efficient implementation to any order L is possible through a recursion relation for modified spherical Bessel functions.

We now return to the master equation, eq. (18), which presents a radial and a double angular integral to be solved. The radial integral is tackled first, and can be written as where p = 2αR and S = λ + μ + υ + L + 2 simplify the notation of the integral. This integral does not appear to have been solved based on ready inspection of published solved integrals. Therefore, the integral had to be creatively reduced to known integrals. The explicit derivation and validation has been siphoned off to Part D in the Supporting Information, but we quote the final result here, where Γ(a,x) is the upper incomplete gamma function (or of the second kind). Part D of the Supporting Information shows how this solved integral cascades up, via a recursive relation, to the solution of . Note that if a is an integer then Γ(a,x) returns essentially an exponential, while if a is a half‐integer then the error function emerges. Also note that the sign must be explicitly preserved when taking the square, because if a square root is taken of this square product it is important that the original sign is preserved.

The second sizeable task in solving the master integral in eq. (18) focuses on , the angular part counterpart of the radial integral . To shorten this article even more, the explicit derivation and validation has been siphoned off to Part E in the Supporting Information. This derivation solves the angular integral analytically, splitting it into an integral in θ and one in φ, that is,

The result for the integral in φ is

It is very important to realize that I vanishes if m > λ + μ + M.

In the calculation of θ, it proved crucial to use a finite expansion of the associated Legendre polynomials in terms of cosines and sines, to obtain a closed and more importantly, universal expression, which is

In summary, by substituting eqs. (24) and (25) back into eq. (23), the integral in the master equation eq. (18) has been solved. Thus, we obtained a general formula for an arbitrary spherical harmonic multipole moment generated by the electron density, within a sphere with radius β, and constructed from Gaussian of arbitrary angular momentum.

The next derivation focuses on the kinetic energy, starting with one type of kinetic energy density, usually denoted G(r) but here written as E G(r) to avoid confusion with the Gaussian primitive function. Because E G(r) contains only first derivatives the expressions are less complex than for another type of kinetic energy density called K(r), now denoted E K(r) for consistency. The latter contains second derivatives and is discussed later. The contribution of the β sphere to the integrated kinetic energy, denoted E G(β), is defined as follows: where all symbols in common with eq. (3) mean the same. The derivation leading to a solution for eq. (26) is given in Part F of the Supporting Information. Thus, the equivalent of eq. (11), which formulates the key integral for Q , then becomes

This equation concludes the proof that the G‐type kinetic energy can be calculated by closed expressions because the integral in eq. (27) is a special case (i.e., L = M = 0 in R(r)) of the integral in eq. (11). Hence, the treatment of the integral in eq. (27) is identical of that in eq. (11). Finally, it should be pointed out that the gradient operator never appeared expressed in spherical polar coordinates. Indeed, its differentiation was carried out at Cartesian level, prior to the treatment of the complex integrals in spherical polar coordinates. The same strategy is followed for the second type of kinetic energy denoted E K(r).

The next derivation (or rather outline thereof) focuses on E K(r), where the contribution of the β sphere to the integrated kinetic energy, denoted E K(β), is defined as follows:

One can prove that the Laplacian operator operating on a Gaussian primitive again results in a linear combination of Gaussian primitives:

This equation allows one to write the equivalent of eq. (F4) of the Supporting Information: where the explicit product between G and is now not written out but seen to give rise to a 7‐term sum with its own pre‐factors F, and pattern of powers (l, m, n). As for E G, we are now re‐assured that we make contact with the previously established machinery of solving the arising integrals. This completes the treatment of kinetic energy. Note that the β‐sphere contribution of the Laplacian of the electron density can be trivially obtained from the two kinetic energies because

To complete the treatment of IQA energies by their β‐sphere contributions, the final section focuses on the potential energy, both within the same topological atom (intra) and between two different topological atoms (inter). An efficient strategy to calculate potential energies builds on two insights: (i) the calculation of the electrostatic potential (generated by a given electron density) as a useful intermediate physical quantity, and (ii) the exploitation of the exact convergence of the multipole expansion (due to the finite volume of the β‐sphere).

Figure 2 sets the scene, introducing symbols that will appear in the equations below. A topological atom, denoted Ω, consists of two regions: the β‐sphere and the remaining region outside the β‐sphere, which we call “α.” All integration in the α‐region is done by numerical quadrature, typically by a Gauss‐Legendre for the radial part and a Lebedev grid for the angular part. The integration of the complex shape of the α‐region could be possible analytically as well, by introducing so‐called natural coordinates in eq. (18), but actually achieving this is perhaps more of a retirement project. The point P marks a possible position where one wants to know the electrostatic potential generated by the electron density in β A.

Figure 2

Electrostatic potential and potential energies within and between two topological atoms ΩA and ΩB. Note that for convenience the global origin coincides with the nuclear position of ΩA. [Color figure can be viewed at wileyonlinelibrary.com]

The electrostatic potential at a given point P, generated by the total charge density inside the β sphere, is defined by where ρ tot(r) is the sum of the electronic and nuclear charge density, described by the position vector r, ρ(r) the purely electronic density featuring throughout this article [starting with eq. (1)], and Z is the charge of the nucleus inside the β sphere. We take advantage of the following well known expansion

In this expression, is the smaller and is the larger value and the common choice is to assign r = |r| to , and to assign P = |P| (where both vectors r and P refer to the same global origin) to . This expansion will then converge, provided r < P, which means that the point P must lie outside the β sphere. This is perfectly possible because the electron density is confined to the β sphere. In other words, formal convergence is only possible if the point (P) at which the potential is evaluated lies further from the nucleus than any element of electron density within the β‐sphere. Substituting eq. (33) into eq. (32) leads to which proves that the electrostatic potential can be calculated directly from the previously computed multipole moments Q(β), without introducing new integrals. The calculation of the potential energy between two topological atoms benefits from this straightforwardly calculated electrostatic potential because this energy can be written as follows: where is the electrostatic potential generated by atom B in atom A.

It is important to fully exploit the fact that the electrostatic potential generated by the β sphere is analytically computed in any point outside it. Because each atom is a union of β sphere and an α region, or Ω = α U β, there are four possible interactions, between two different atoms (A ≠ B): β A – β B, β A – α B, α A – β B, and α A – α B. This division can be summarized as follows: where is a purely geometric interaction tensor and R AB is the internuclear vector. In the second and third term, the quadrature grid over the respective α‐regions of atoms B and A will ask for values of the analytically evaluated electrostatic potential in each of its points. The fourth term is purely numerical and hence does not benefit from the work of the current article. There remains only one case to be mentioned, which is that of A = B, because it causes its own computational regime, but only for β A – β A. Indeed, the interactions β A – α A and α A – α A can be calculated without extra knowledge. For the β A – β A case, one option is to follow the method “inverse multipole moments,”54 another is to introduce a Fourier transform that is used in the evaluation of nuclear‐electron attraction integrals.

Conclusions

We present a detailed derivation of a fully analytical 3D integration over the volume bounded by the β sphere inside a topological atom. The formulae are sufficiently general such that a Gaussian primitive of arbitrary angular momentum can contribute to the electron density within the β sphere. We obtained a general formula for an arbitrary spherical harmonic multipole moment, the kinetic energies and the potential energy between atoms. We showed that closed expressions can be obtained involving the exponential, error function, and gamma functions. The implementation of the results of the derivations given here will follow in a subsequent publication. This implementation will make a choice of quadrature grid superfluous, but it is not yet clear how the current formulae will perform computationally, compared with quadrature schemes.

Supporting Information

Click here for additional data file.