A multipolar approach to the interatomic covalent interaction energy

Interatomic exchange‐correlation energies correspond to the covalent energetic contributions to an interatomic interaction in real space theories of the chemical bond, but their widespread use is severely limited due to their computationally intensive character. In the same way as the multipolar (mp) expansion is customary used in biomolecular modeling to approximate the classical Coulomb interaction between two charge densities ρA(r) and ρB(r) , we examine in this work the mp approach to approximate the interatomic exchange‐correlation (xc) energies of the Interacting Quantum Atoms method. We show that the full xc mp series is quickly divergent for directly bonded atoms (1–2 pairs) albeit it works reasonably well most times for 1– n (n > 2) interactions. As with conventional perturbation theory, we show numerically that the xc series is asymptotically convergent and that, a truncated xc mp approximation retaining terms up to l1+l2=2 usually gives relatively accurate results, sometimes even for directly bonded atoms. Our findings are supported by extensive numerical analyses on a variety of systems that range from several standard hydrogen bonded dimers to typically covalent or aromatic molecules. The exact algebraic relationship between the monopole‐monopole xc mp term and the inter‐atomic bond order, as measured by the delocalization index of the quantum theory of atoms in molecules, is also established. © 2017 Wiley Periodicals, Inc.


Introduction
The role of the quantum mechanical exchange-correlation (xc) energy as the basic glue binding together atoms and molecules has been clearly stressed in the past. [1] In the chemical literature, however, this insight is less well-known. Although exchange-correlation functionals, for instance, are the essential ingredients in modern implementations of Density Functional Theory (DFT), [2] not much work has been devoted to examine the importance of the exchange-correlation energy itself in the theory of chemical bonding from the DFT viewpoint. [3] Actually, almost all that is known about the chemical relevance of the xc energy has been derived in the last decade through the study of bonding in real or position space. [4,5] With this term, we gather together a number of techniques that are being actively explored [6][7][8] which use orbital invariant reduced densities (or density matrices) to develop a new paradigm that may one day replace the standard molecular orbital (MO) approach. [9] Usually, these techniques use a partition of real space into regions endowed with chemical meaning, be them atoms, bonds, cores, lone pairs, and so forth. In many cases, the space is divided using the topology induced by the gradient field of an orbital invariant scalar, like the electron density (which gives rise to the atomic partioning of the quantum theory of atoms in molecules (QTAIM) developed by Bader, [10] or the electron localization function (that isolates core, bond and lone pair regions). [11,12] When this topological tools are used, we say that we are under the Quantum Chemical Topology umbrella. [13] In the context of the QTAIM/QCT, we proposed a number of years ago an exact, general decomposition of the total molecular energy E into atomic and inter-atomic terms that we called the interacting quantum atoms (IQA) approach. [4,5] All the expectation values of the standard Coulomb Hamiltonian that make up E are written in IQA as a sum of domain contributions, and E is obtained by adding atomic self-energies, which tend to the free atomic energies when the atoms that interact are sufficiently far apart, and pairwise additive interaction energies. The latter are composed of a classical term that depends only on classical electrostatic contributions, and an exchange-correlation energy, V xc which accounts for the quantum mechanical effects. As we and others have shown over the years, [14,15] the classical part of the interaction measures its ionic component, while the V xc energy is to be associated with its covalent counterpart.
In these years, the interatomic xc energy has become an important ingredient of any quantitative account of chemical bonding in position space. [16,17] For instance, it has been shown to be intimately related to the appearance of the bond critical points of the QTAIM, leading to the concept of priviledged exchange-correlation channels. [18] It has also been used to reconstruct molecular graphs from purely energetic quantities, [19] to shed light on new concepts like halogen bonding, [20,21] to recover stereolectronic effects, [22] or to find new long-range electronic anomalies. [23] Interatomic V xc energies are intimately linked to the delocalization or shared electron delocalization indices (DIs) used in the QTAIM, defined almost 40 years ago by Bader and Stephens. [24] These are obtained by directly integrating the xc density of very two different atomic domains and measure the number of shared pairs of electrons between them. They have been successfully used as real space generalization of the bond order concept, reducing to the Wiberg-Mayer [25,26] bond orders if atomic domains are imagined to collapse onto their nuclei. In a sense, V xc 's are the energetic counterparts of DIs, and both have been empirically found to correlate very well when a given couple of atoms is examined in different molecular environments.
The computational complexity of obtaining DIs is considerably smaller than that of calculating V xc 's, as the former may be factorized into sums of products of atomic overlap matrices (AOM; 3D numerical integrals), while the latter need, in principle, very costly 6D quadratures. Thus, if we are not interested in very accurate results, but only in semi-quantitative estimations of covalent energies, any procedure that might approximate the V xc values in terms of cheaper to compute quantities like the DIs should be welcome. That procedure was initially examined by Rafat and Popelier, [27] that wrote each interatomic Hartree-Fock (HF) V xc interaction as a multipolar series, exploring the convergence of this series in different closedshell molecules computed at this level of theory. In this work, we generalize their algebraic formalism to multi-determinant wavefunctions. This generalization is possible thanks to the use of the monadic diagonalization of the exchangecorrelation density, [28] customary used within the IQA methodology. [4,5] Our expressions converge to those of Ref. [27] when HF exchange-correlation densities are used in the calculation. We will show that, regardless the type of calculation, the monopole-monopole term of the multipolar xc interaction between two atoms A and B of the molecule coincides with that of Rafat and Popelier, being equal to 2d AB =ð2RÞ (R is the AB internuclear distance), and that the series is usually divergent although many times asymptotically convergent. Moreover, our results clearly establish in what conditions V xc can be safely approximated by a truncated series, and how in some situations retaining up to the charge-quadrupole terms may give reasonable results even for directly bonded atoms. In the latter cases, the use of the crudest approach V AB xc $ 2d AB =ð2RÞ is even preferable to using the multipolar expansion up to very high order.
We will first consider the multipolar expansion of V xc , including a short account of the IQA methodology. Then, we will turn to examine how the series converges or diverges for a number of selected systems.

Multipolar Expansion of V AB xc
In this section, we briefly describe the IQA method and the role played by the exchange-correlation (xc) interaction in this energy partition method (Subsection), the exact computation of this interaction (Subsection), and its multipolar approximation (MP) with or without truncating the expansion of the angular momentum series (Subsection). It is worth noting that the experience gained to date with the IQA method, both by us and by other groups, clearly indicates that the magnitude of V AB xc correlates very well with the degree of covalency between the pair of atoms A and B as measured by means of the DI defined by Bader and Stephens, and weighted through the inverse of the distance between both atoms. As we will see, this correlation would be perfect as long as the crudest MP to V AB xc (consisting in truncating the multipolar series in the term l 1 5m 1 5l 2 5m 2 50) were exact.

The IQA method
The IQA method [4,5] is a real space energetic partition inspired in the QTAIM that focuses on domain-averaged integrated quantities. The total energy in this approach is given by where A runs over all the atoms in the molecule, V AB nn 5Z A Z B =R AB is the repulsion between the nuclei A and B, V AB en 52Z B Ð XA d r 1 qðr 1 Þr 21 1B is the nuclear attraction of the electrons within the basin of A (X A ) to the nucleus B, and V AB ee is the total electron repulsion between X A and X B . The latter is given by is the classical or Coulomb electron-electrons repulsion, and where q xc ðr 1 ; r 2 Þ is the exchange-correlation (xc) density, is the purely quantum-mechanical electron-electron xc interaction, which is the main subject of this work. In this way, The term E A self in eq. (2) collects all the energetic components affecting exclusively to the atom A while E AB int represents the full interaction energy between atoms A and B, that is made of the full electrostatic or classical interaction (V AB cl ) and the quantum-mechanical part (V AB xc ). The expression 2 is valid, not only for the IQA methodology but also for other energetic partitions, such as a recently proposed one inspired in the IQA method, although using a fuzzy partition of the space and localized MO. [29] Mayer and Hamza have also dealt with the exchange component in eq. (4) in the framework of a Hilbert space partition instead of the real space QTAIM partition we use here. [30] The exact xc interaction energy Over the years, it has become clear that the magnitude of V AB xc measures the degree of covalency of the chemical bond between the atoms A and B. The more negative its value, the bigger the bond order between the two atoms and vice versa. [6,14,18] Their values have been recently proposed as a novel solution to the problem of assigning a molecular graph to a collection of nuclei [23] (i.e., how to draw a molecular structure). In the IQA approach, this term is exactly computed as follows. First, we use the fact that for both single-(1-det) and multi-determinant wavefunctions built in with real MOs / i , q xc ðr 1 ; r 2 Þ can be written as where M is the number of partially or fully occupied MOs, and k ijkl is a symmetric matrix in the (i, j) and (k, l) pairs. Defining a set of coefficients, where d ij is the Kronecker symbol (d ij 51 for i 5 j, d ij 50 for i 6 ¼ j) we may write an (i, j),(k, l) symmetric simpler expression, Using the basis of products of MOs, f/ i ðrÞ/ j ðrÞ; i ! jg, that contains MðM11Þ=2 members, we diagonalize eq. (8), and get [28] : where the f ij eigenfunctions are linear combinations of the above products. The E matrix may be easily computed from the explicit form of a given calculated wavefunction. For closed-shell 1-det wavefunctions (and formally also for a Kohn-Sham determinant) M5N=2, where N is the number of electrons, the E matrix is already diagonal in the (i, j) and (k, l) pairs, each eigenvector is the product of two MOs, f ij 5/ i / j , and the g ij eigenvalues are simply g ii 522 and g ij 524 (i 6 ¼ j). Using eq. (9) in the expression of V AB xc one gets The integrals 3 and 11 can be computed numerically and (in principle) exactly, that is, without invoking any approximation such as the multipolar expansion, by means of the bipolar expansion as described in Ref. [31]. Notice that using the Fock-Dirac exchange from Kohn-Sham determinants is an approximation that has no rigorous justification.
The multipolar approach for V AB xc Comparing eq. (11) with eq. (3) for the Coulomb repulsion it is evident that if J AB is approximated making use of physically reasonable arguments that are also valid for K AB ij , the steps to approximate the latter will be the same used for J AB . The longrange or MP to J AB , given by where m 1 (m 2 ) runs from 2l 1 (2l 2 ) to 1l 1 (1l 2 ), Fig. 14) is the position vector of the B center with respect to the A center, C l1m1l2m2 ðRÞ are known coefficients, Q X lm are the spherical atomic multipoles, defined as Q X lm 5N l ð X r l S lm ðrÞ qðrÞ dr; N l 5 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4p=ð2l11Þ p , and S lm ðrÞ are real spherical harmonics (see Appendix) is exact when the basins X A and X B are nonoverlapping (see Fig. 1 and the definition of overlapping and non-overlapping regions below). Equation (12) is the same  used by Popelier et al. [33][34][35] in their discussion of the multipolar expansion for the diatomic Coulomb repulsion. Retaining only terms with l 1 1 and l 2 1 in this equation one has J AB lr;cd ' where Q X 5 Ð X qðrÞdr andl X 5 Ð X rqðrÞdr are the total electron charge and the dipole moment of the X region, respectively. The first, second plus third, and fourth terms of 14 correspond to the charge-charge (cc), charge-dipole (cd), and dipole-dipole (dd) interactions, respectively. We should note that the second and third terms have opposite signs.
If the same approximation is used for It is important to stress that, similarly to ðJ AB Þ lr , the expression 15 for ðV AB xc Þ lr provides the exact xc interaction when the atomic basins X A and X B do not overlap (Fig. 1). In the present context, these two basins are non-overlapping because the two spheres of radii R A and R B , centered at the origin of X A and X B , respectively, do not intersect each other, being R A (R B ) the maximum distance from the origin of the basin to the surface of X A (X B ). On the contrary, X A and X C are overlapping despite that none point inside X A belongs also to X C and vice versa. When the non-overlapping condition is not met, the current expressions for ðJ AB Þ lr and ðV AB xc Þ lr are only conditionally convergent. We will see different examples of this in section. The function N l r l S lm ðrÞ is 1 for l5m50, (y, z, x) for l 5 1 and m5ð21; 0; 11Þ, and ð ffiffi ffi ffiffi ffi 3 p xz; ffiffi 3 p 2 ðx 2 2y 2 ÞÞ for l 5 2 and m5ð22; 21; 0; 11; 12Þ. If, as in the case of J AB , only terms with l 1 1 and l 2 1 are included, ðV AB xc Þ lr becomes where q X ij q X ij;00 5 ð X f ij ðrÞdr, and If terms with ðl 1 50; l 2 52Þ and ðl 1 52; l 2 50Þ are also included, the extra contribution must be added to 18. The cq subscript in eq. (20) stands for charge-quadrupole interactions. The improved expression for ðV AB xc Þ lr is then ðV AB xc Þ lr;cdq 5ðV AB xc Þ lr;cd 1ðV AB xc Þ lr;cq : The physical meaning of q X ij andl X ij are easy to grasp. If we consider the particular case of their diagonal expressions (i 5 j) for a 1-det wavefunction, f ii ðrÞ5/ 2 i ðrÞ, so that q X ii is the electron charge of the orbital distribution / 2 i ðrÞ within the X region, andl X ii the dipole moment of X due to this distribution. For this reason, q X ij andl X ij may be called orbital overlap charge and orbital overlap dipole, respectively. At the HF level,  the q X ij 's coincide with the AOM elements of the QTAIM, However, given that f ij ðrÞ at the correlated level is a linear combination of / i ðrÞ/ j ðrÞ products, q X ij in this case is a linear combination of AOM elements. Nevertheless, for both types of wavefunctions 22d AB 00;00 coincides with d AB , the so-called DI between the atoms A and B 22d AB 00;00 5d AB 52 so that the leading term of ðV AB xc Þ lr ðl 1 5l 2 5m 1 5m 2 50Þ can be written as The above equation is behind the good existing correlation between the values of V AB xc and d AB for a large collection of AB couples in many systems. The present derivation shows that the proportionality between V AB xc and d AB is modulated by the inverse of the distance between the nuclei of both atomic basins.

Systems and Computational Details
All the calculations of this work have performed with our PRO-MOLDEN code. [36] This program allows the exact computation [28,31] (i.e., without suffering the convergence problems inherent to the multipolar series expansion) of V AB xc as well as the full (lr) and truncated (lr,cd) and (lr,cdq) MP described in section. For brevity, only the exact, and the (lr) and (lr,cdq) numbers will be given in the tables. The errors plotted in the figures are defined as ½ðV AB xc Þ method 2ðV AB xc;exact Þ=jV AB xc;exact j3100, where method5(lr), (lr,cd), or (lr,cdq). The studied systems include several standard hydrogen bonded (HB) dimers (Fig. 2), the staggered BH 3 NH 3 , eclipsed BH 3 NH 3 , N 1 5 , and Li 9 H 9 molecules ( Fig. 3), 11 molecules derived from saturated hydrocarbons by substituting C or H atoms by Be, B, N, O, F atoms, plus the benzene molecule (Fig. 4), the saturated hydrocarbons ethane, propane, butane, and pentane (Fig. 5), and the phenol dimer (Fig. 6). The labels of the atoms in the tables are those defined in these figures. For simplicity, the MO required for evaluating the exchange-correlation density of eq. (8) have been obtained through restricted Hartree-Fock (RHF) calculations at the corresponding equilibrium geometries with basis sets of quality 6-311G(d,p) or higher. However, as our results in this article stem from the algebraic properties of the multipolar expansion, we do not expect significant changes neither in the numerical results nor in the subsequent discussion when using more accurate wavefunctions or the approximate data coming from Kohn-Sham determinants in the computation of the xc interactions. To prove the validity of this assertion, we will compare the V AB x energies obtained for staggered ethane in a CAS [14,14] calculation (Complete Active Space calculation with all except the carbon 1s electrons distributed into 14  valence orbitals) with the RHF results. All other xc interactions except those of the above CAS calculation lack a correlation energy component, being thus pure exchange contributions that should be more properly labelled V AB x . However, as all the expressions in section are valid for general wavefunctions the original name will be used hereinafter. The sums over l 1 and l 2 in eq. (15) were truncated at l max 1 5l max 2 58, so that terms up to a range L5l max 1 1l max 2 11517 were included in the multipolar expansion. As QTAIM domains are usually finite and quite irregular, very fine radial and angular grids are needed to carry out the 6D numerical integrations. Here, we have systematically considered a b2sphere around each atom, with a radius equal to 60-90% the distance of its nucleus to the closest bond critical point, and used high quality Lebedev angular and radial grids, with (5810, 512) and (194, 400) points outside and inside the spheres, respectively. The errors in the total energy of the studied molecules attributable to these numerical integrations, necessarily approximate, are of the order of 1.0 kJ/mol. Our accumulated experience in IQA calculations makes us believe that the accuracy achieved in each interatomic interaction is even higher. Despite this issue regarding the full numerical accuracy of our integrations, once Table 1. xc interaction energies ! 0:1 kJ/mol for the HB dimers of Figure 2. the computational conditions of a given calculation have been chosen the convergence of the bipolar expansion (the exact benchmark) is ordinarily well below the 1 kJ/mol barrier for the xc contributions.

Results and Discussion
The more representative results regarding the approximate V AB xc values, as well as their errors for the systems listed in Section are gathered in Tables 1-3 and Figures 7-11. We can see in Table 1, where the V AB xc 's for the HB systems of Figure 2 are collected, that the full MP (V AB xc ) lr [eq. (15)] fails miserably for all intramolecular A-H pairs (A5N,O,F). Surprisingly, the crude (lr,cdq) approximation gives xc interactions with relative errors of about 10% or smaller for the intramolecular directly bonded atoms. Regarding the intermolecular interactions, the xc energy between the two atoms involved in the HB is well represented by (V AB xc ) lr , with differences with respect to the exact values smaller than 0.3 kJ/mol in all the cases. We note again that the (lr,cdq) values differ only by 0.3-0.6 kJ/mol from the exact ones, confirming that the multipolar expansion for these interactions is practically converged at this level of calculation. Intermolecular A-H and H-H xc energies other than the above ones are given by the lr approximation with errors smaller than 0.1 and 0.01 kJ/mol, respectively. For the intermolecular H-H energies, the same is true in the (lr,cdq) approximation. However, the xc interaction between the A atom of the proton donor (PD) and the B atom of the proton acceptor (PA) molecule is predicted with errors as large as 1.4 kJ/mol (FHÁ Á ÁNH 3 ) when the (lr,cdq) approximation is used, which clearly indicates that multipolar interactions higher than the chargequadrupole ones included in this approximation are required to represent this type of interaction with accuracy.
The relative errors of the A-B, A-H (A,B5N,O,F), and H-H xc interaction energies for all the intramolecular and intermolecular pairs of the HB systems are represented in Figure 7. We observe in Figure 7 that cd and cdq intramolecular H-H energies are, in general, more accurate than the ones for the A-H interactions, which is clearly due to the 1-3 (1-2) character of all the intramolecular H-H (A-H) pairs. It is also striking that, with a couple of exceptions, cdq relative errors are negative whereas the contrary happens with the cd approximation. Moreover, as previously commented, only a single lr relative error appears in the figure, the remaining ones having errors greater than 20%. Regarding the intermolecular xc energies, we observe in right Figure 7 the progressive decreasing of relative errors in passing from cd to cdq, and from cdq to lr.
The discussion for the systems in Figure 3 runs parallel to that of the HB dimers. The (V AB xc ) lr value for the B-N pair in eclipsed and staggered BH 3 2NH 3 has no sense. Similarly, the lr xc interaction between the directly bonded (i.e., 1-2) B-H and N-H pairs is quite absurd. Not only that, but also the (V AB xc ) lr 's for the 1-3 pairs H 3 2H 4 and H 6 2H 7 are several orders of magnitude greater than the exact values. Contrarily to this, In Li 9 H 9 , the directly bonded Li and H atoms are signalled with (1).

FULL PAPER
WWW.C-CHEM.ORG Table 3. Representative xc interaction energies (kJ/mol) for the systems of Figure 4. the cdq approximation works relatively well for the B-N pair and the 1-2 B-H and N-H pairs (relative error <5%). The H 3 2H 4 interaction is also extremely well reproduced by this approximation (error <0:2%), whereas the H 6 2H 7 is slightly worse (error $4%). The xc interaction between the B atom and a H atom of the NH 3 unit (B 1 2H 6 ) is fairly accurate in both the lr and cdq approximations. This is not so with the symmetric interaction N 2 2H 3 , with errors about 4-7% in both cases.
In the N 1 5 molecule, the xc energy for the 1-2 pairs is again badly represented by the lr approximation but is reasonable in the cdq approach, particularly for N 1 2N 2 . The cdq and lr values for the 1-3 N 1 2N 4 interaction differs from the exact value by about 0.9 and 0.3 kJ/mol, respectively. The error in the other 1-3 interaction (N 2 2N 3 ) in considerably higher in both approaches. Finally, the lr and cdq xc energies for the 1-4 N 2 2N 5 and 1-5 N 4 2N 5 pairs are practically the same and coincident with the exact value. This result highlights two important facts: (i) the atomic basins of N 2 and N 5 (or N 4 and N 5 ) atoms fullfil almost exactly the non-overlapping criterion displayed in Figure 1, and (ii) the multipolar series 15 converges very quickly in this particular case.
Finally, the results for Li 9 2H 9 reinforce what was said in the above three paragraphs. The full lr expansion fails completely in predicting xc interaction energies for 1-2 Li-H pairs, while the cdq values are pretty accurate. All Li-Li xc energies are well represented in the lr and cdq approximations, with the exception of the lr Li 2 2Li 3 interaction. This is probably related with the almost spherical character of Li atomic basins. According to this, the 1-3 Li-H lr xc energies and, more importantly, the 1-3 H-H xc energies are less accurately computed due to the far from the spherical character of H atomic basins. This is exacerbated in the lr approximation, where higher angular number l values are involved [see eq. (17)].
The xc pair interaction energies of the systems in Figure 4 are collected in Table 3, and the relative errors of the cd, cdq and lr approximate values displayed in Figure 8, for the 1-2 (left-top), 1-3 (right-top), and 1-4 (bottom) pairs, respectively. Virtually all of the above comments also apply here: the 1-2 xc interactions cannot be represented at all using the full lr  expansion. However, they are given with reasonable accuracy by the cdq approximation. The 1 -n (n > 2) interactions are gradually better reproduced as n increases in both the lr and cdq approximations. It is very satisfactory to check that the cdq approach, a severe truncation of the full mp expansion, is perfectly suited to simulate the xc interaction between pairs of atoms beyond the directly bonded ones. Even in typically covalent molecules like benzene all the cdq C-C xc interaction energies reproduce very well the exact values. As we can see in Figure 8 most of the 1-2 interactions have relative cdq errors 10%. This improves for the 1-3 and 1-4 interactions. A summary of our results for the saturated hydrocarbons C n H 2n12 (n5225) is presented in graphical form in Figure 9. We find the surprising result that all 1-2 C-C cdq interactions are predicted with errors 2% while the cdq energies between the more distant 1-3 C-C pairs have errors about 5-6%. Nevertheless, the interactions between even more distant C-C pairs turn again to be calculated quite accurately (errors <1%) in the cdq approximation. The lr approximation fails completely to predict the 1-2 C-C interactions, but yields negligible errors for the 1-3 xc interaction energies. With regard to the C-H interactions, the situation is the opposite of that found for the C-C pairs: 1-2 C-H cdq errors are about 9-10% (except in ethane where the error is unusually large (54%)) whereas all except two of the 1-3 C-H cdq errors are <1%. For these two exceptions the error is not too large ($1:3%). We observe in Figure 9 that the cdq approximation improves considerably the cd results, giving 1-3 C-H interaction energies almost as accurate as the lr ones. Another surprising result in these systems concerns the cdq 1-3 H-H interactions: Contrary to what happens almost systematically, the cdq results are worse than the cd ones, albeit the relative errors in both approximations are acceptable ($223%).
We have considered staggered ethane as a representative example to analyze the correlation effects on the V AB xc energies. Our results for five representative AB pairs of this molecule are collected in Table 4. Correlation decreases (increases) the magnitude of the C-C interaction (H-H interactions), changes very little V CH xc when C belongs to a CH 3 group and H to the other, and enhances the intra-group V CH xc energies. The discussion of the above paragraph for the RHF results is still approximately valid for the correlated calculation. With the exception of inter-group H-H interactions, the cdq xc energies are closer to the exact V AB xc values than their full MP (V AB xc ) lr . As in many of the HF 1-2 interactions, the lr approximation fails to predict even a reasonable value for the C 1 2H 5 xc energy.
The different behavior of the lr and cdq approximations can be further illustrated with the case of the phenol dimer (Fig.  6). For this system, the relative errors versus the interatomic distance R A2B in these two approximations are plotted in Figure 10, both for intramolecular and intermolecular atomic pairs. Only two points, associated to intramolecular interactions, have a relative error (absolute value) ! 20% in the cdq calculation, while the error for all the lr points with R A2B < 2:64 (most of them associated to intramolecular pairs) is larger than 20%. However, for R A2B > 5:0 the lr approximation gives quite accurate xc interaction energies for all the pairs, whereas cdq errors are still important.
However, there is a general problem of the lr approximation that deserves to be commented: eq. (15) does not necessary converges to the exact xc interaction for large l 1 m 1 and l 1 m 1 values. This fact is illustrated in Figure 11, where the xc energies for some of the atomic pairs of the molecules in Figure 4 are represented versus l max 1 1l max 2 . The lr approximation suffers a systematic error in the C 1 2F 6 interaction of cis-CH 2 CF 2 and trans-CH 2 CF 2 molecules, regardless the value of l max 1 1l max 2 . In the case of trans-CH 2 CF 2 , the C 1 2F 6 xc interaction energy shows an oscillating behavior around an (erroneous) mean value. This pattern has been also observed in other cases. Contrarily, as we have repeatedly said in this section, catastrophic lr interactions (see, for instance Fig. 11) are still reasonable provided that the sum l 1 1l 2 is interrupted at a value approximately in the interval 2 l 1 1l 2 6.
This simple analysis shows that conclusions on the convergence of the mp expansion drawn exclusively from limited L5 l max 1 1l max 2 11 data cannot be trusted. Figure 12, where the ðV AB xc Þ lr energies for some 122 interactions of the H 2 O-H 2 O, HF-HF, and NH 3 2NH 3 dimers are shown, illustrates this fact in a crystal clear way. Cutting the mp expansion at L 8 the oscillatory behavior of Figure 11 would have been found indeed, but no catastrophe for larger L values would have been predicted. However, for L > 8 the mp expansion progressively deteriorates and for high L values the ðV AB xc Þ lr energies diverge. Again, the cdq approximation works fairly well, with all the 1-2 interactions in Figure 12 predicted with relative errors smaller than 3% except F 1 2H 2 (3.2%) and F 3 2H 4 (5.2%).
Summarizing, we have shown that the multipole series for the interatomic xc energies is conditionally convergent, and that the computational burden of the quasi-exact calculation of V xc when general QTAIM domains are used may be ameliorated by retaining up to quadrupole-quadrupole terms. With this approximation, reasonable errors are obtained in the medium-to long-distance range, sometimes even for directly bonded interactions. An overall image of the improvement of the cdq approximation over the cc (monopole-monopole) one can be grasped from Figure 13 that condenses all our calculations that span a six orders of magnitude range for V xc .
We have explored numerically the degree of fulfillment of the linear relation between V AB xc and d AB =ð2RÞ. Assuming V AB xc ' 2a½d AB =ð2RÞ1b, we have determined a and b for each molecule by fitting the computed values of V AB xc and d AB for every AB pair of this molecule to the above expression, verifying that b is always very small and a takes values relatively close to (but smaller than) 1.0. As representative examples, (a, b) for the H 2 O-H 2 O, NH 3 -NH 3 , and HF-HF dimers are ða; bÞ5ð0:8988; 20:0005Þ; ð0:8673; 0:0000Þ, and ð0:9150; 2 0:0010Þ, respectively. This result opens another possible route to the approximate but much cheaper computation of the xc interaction in those cases where the exact calculation is prohibitive or very expensive. An extensive analysis of this correlation in anion-p interactions which corroborates the above statement has been carried out by Foroutan-Nejad et al. [37][38][39]

Conclusions
We have shown that the interatomic exchange-correlation energies used in real space theories of chemical bonding, which measure the covalent contribution to a given interatomic interaction, can be approximated via a conventional multipole expansion. Rigorously, the series diverges when atoms are directly bonded, although it may be regarded asymptotically convergent. Truncation of the series up to l 1 1l 2 52 (including up to charge-quadrupole interactions) tends to provide results which are accurate to a few percent in 1 2 n, n > 2 interactions, and even to about 10% in many 1-2 directly bonded cases. In the n > 2 case the series converges in many cases, and including extra terms provides further accuracy. On the contrary, the consideration of larger l contributions in 1-2 interactions tends to seriously deteriorate the results. As the computational burden needed to calculate the multipole series is considerably smaller than that of the exact bipolar expansion, our results may be important to estimate covalent interactions in those cases where exact integrations are not feasible. They can also be used to ameliorate the computational cost in IQA decompositions of large systems, where many expensive, but small long-range xc terms can now be safely approximated without loss of precision. 1 m50; ffiffi ffi 2 p sin jmj/ m < 0; 8 > > < > > : (27) and P m l ðcos hÞ are the associated Legendre functions, defined for m ! 0 by P m l ðxÞ5 1 2 l l! ð12x 2 Þ m=2 d l1m dx l1m ðx 2 21Þ l : Finally, D l2m2 l1m1 ðr 1 ; r 2 ; RÞ in eq. (24) is defined as D l2m2 l1m1 ðr 1 ; r 2 ; RÞ54pð21Þ l1 X l11l2 l35jl12l2j V l1l2l3 ðr 1 ; r 2 ; RÞ T l3 l1m1l2m2 ðRÞ; (29) where the sum over l 3 runs in steps of 2, V l1l2l3 ðr 1 ; r 2 ; RÞ is a discriminant that takes different expressions in the four regions defined in Figure 15, and T l3 l1m1l2m2 ðRÞ is the angular factor T l3 l1m1l2m2 ðRÞ5 X 1l3 m352l3 d l3m3 l1m1l2m2 S l3m3 ðRÞ; where d l3m3 l1m1l2m2 is the Gaunt coefficient between the S lm ðh; /Þ's defined by d l3m3 l1m1l2m2 5 hS l3m3 jS l1m1 jS l2m2 i: Given that S lm is real, d l3m3 l1m1l2m2 is invariant against any permutation of the pair of indices (l i , m i ). These coefficients may be determined as described elsewhere. Using eq. (24) in eq. (11) one has A further simplification of K AB ij requires the explicit form of D l2m2 l1m1 ðr 1 ; r 2 ; RÞ. From the expression of V l1l2l3 ðr 1 ; r 2 ; RÞ [eqs. (B1)-(B9) of I] it follows that, as long as r 1 1r 2 R, this discriminant only takes a nonzero value in region III of the (r 1 , r 2 ) space (see Fig. 15). This condition will be exactly satisfied if R ! r max 1 1r max 2 , where r max 1 is the maximum value of the radial coordinate within X A , with an equivalent definition for r max 2 . In the present context, atoms A and B are said to be nonoverlapping if this condition is fulfilled, and overlapping otherwise. Although it may occur that the condition R ! r max 1 1r max 2 is not exactly satisfied, provided that the atomic basins X A and X B are well-separated in the space, we can expect that it is fulfilled in practical terms. The multipolar approach, intensively used to approximate the Coulomb repulsion in the modelization of biomolecules, is equivalent to the assumption that r 1 1r 2 R for any r 1 and r 2 . Thus, region III is identified with the complete first quadrant. In this region, D l2m2 l1m1 ðr 1 ; r 2 ; RÞ is given by D l2m2 l1m1 ðr 1 ; r 2 ; RÞ5ð21Þ l1 16p 2 D l1l2 r l1 1 r l2 2 R l11l211 T l11l2 l1m1l2m2 ðRÞ; where (33) D l1l2 5ð21Þ l11l2 ð2l 1 12l 2 Þ! l 1 ! l 2 ! ðl 1 1l 2 Þ! ð2l 1 11Þ! ð2l 2 11Þ! : Using eq. (33) in eq. (32), we get  (36) and the q X ij;lm have been defined in eq. (17). Finally, substituting eq. (35) in eq. (10) we obtain eq. (15), the MP for the exchange-correlation interaction, V AB xc À Á lr .
Keywords: chemical bonding Á atoms in molecules Á covalent energy Á molecular energy partitioning Á electron delocalization