Exact decompositions of the total KS-DFT exchange–correlation energy into one- and two-center terms Open Access

The accurate evaluation of the total energy of a molecular system is one of the most important challenges in quantum chemistry. However, the value of the energy itself provides little immediate chemical information. To extract chemical bonding information from the energy, one can make use of computational tools, which operate on the energetics of the bond/interaction formation, in particular the schemes that decompose the molecular or the formation energy into chemically meaningful terms. The so-called energy decomposition analysis (EDA) schemes are mainly developed to study a particular A–B interaction from the pair of fragments A and B. For that purpose, monomers or molecular fragments A and B must be defined and computed at the electronic state, which would better describe the bonding situation. The EDA terms that are obtained for the interaction energy are thus necessarily dependent on the nature of the reference’s isolated fragments.

Alternatively to that, in the case of the so-called interacting quantum atoms (IQA) approaches, the total energy is exactly decomposed (up to the numerical integration error) into one- and two-center contributions^1–3 solely using the wavefunction from the optimized AB system. The centers can be the individual atoms composing the system or groups of atoms, permitting the identification of energetic interactions between functional groups in the case of a molecular system or individual monomers in a complex. In IQA methods, no additional reference calculations are required, and the intra- and intermolecular interactions present in the complex are treated on equal footing.

Over the past two decades, this methodology has been widely employed to elucidate numerous chemical phenomena with a particular emphasis on the nature of chemical bonding. The descriptors derived from IQA analysis have been used to rationalize trends in the binding energies of first-row diatomics⁴ to investigate into the nature of different bonds, such as hydrogen bonds^5,6 and their cooperative effects,⁷ halogen bonds,⁸ and beryllium bonds,⁹ or, more recently, to propose a novel bonding type (i.e., collective interactions).¹⁰ In addition, efforts have been made to recover classical concepts from the IQA analysis, such as the steric repulsion in rotational barriers^11,12 and SN2 reactions¹³ or the gauche effect,¹⁴ to name a few. For a recent review of the IQA methodology, see Ref. 15.

IQA energy decomposition relies on real-space analysis and, hence, on the identification of the atom within the molecule (AIM), which is often chosen to be that from the quantum theory of atoms in molecules (QTAIM).¹⁶ It is important to stress that the real-space decomposition of properties such as the energy is not restricted to that particular atomic model. An atom A within the molecule may be more generally identified by continuous atomic weight functions w_A(r) ≤ 0 centered in the nucleus of the atom. Such atomic weight functions can be derived from a variety of approaches, including a number of Hirshfeld-type variants^17–19 or schemes that borrow elements of the QTAIM model, such as the topological fuzzy Voronoi cell (TFVC) approach.²⁰ 

Since the total energy can be expressed in terms of one- and two-electron density functions, it quite naturally decomposes into one- (atomic) and two-center (diatomic) contributions by applying the equations above (see in the following). However, contrary to electron distribution analyses (e.g., atomic populations, bond orders, or local spins), the formulation of the molecular energy decomposition scheme depends upon how the total energy is obtained for each particular electronic structure method. Popelier and Kosov first considered two-electron integrations over different QTAIM domains.²¹ Salvador et al. introduced the IQA scheme for the Hartree–Fock energy both for QTAIM and fuzzy-atom AIMs.^1,2 Later, Blanco et al. extended it to correlated wavefunctions³ [e.g., Configuration Interaction Singles and Doubles (CISD) or Complete Active Space Self-Consistent Field (CASSCF)] and introduced an efficient numerical quadrature algorithm for the two-center integrations.^22,23 More recently, proper formulations for Møller–Plesset second-order perturbation theory (MP2)^24,25 and Coupled-Cluster Singles and Doubles (CCSD)^26–28 energies have also been successfully introduced. Curiously enough, the extension of the method to the Kohn–Sham density functional theory (KS-DFT) has proven to be the most challenging, the origin of the problem being the local contribution of the exchange–correlation functional.

The first term is the non-local exchange, which has the form of Eq. (9) but using the Kohn–Sham molecular orbitals (KS-MOs). It is modulated by the parameter a₀, ranging from zero (pure local exchange) to one (HF-like). The correlation part also can be expressed through non-local and local contributions. In double-hybrid functionals, the former borrows the form of the Møller–Plesset second-order energy correction formula again using the KS-MOs. In the case of range-separated functionals, the non-local and local parts of the exchange functional are modified according to the range-separation parameter, but the essence of the two contributions is kept.

According to Eq. (1), $Exclocal$ naturally decomposes into one-center (atomic) contributions. This fact has been recently used by some of us to develop an origin-independent decomposition of the electronic polarizability into atomic/fragment contributions.²⁹ However, for chemical bonding analysis, this situation is clearly unsatisfactory. First of all, since the HF-like non-local part of $ExcDFT$ does decompose into both one- and two-center terms, an IQA-type analysis would render a completely different picture of the atomic and diatomic interactions within the molecule when using a pure KS-DFT functional as compared to Hartree–Fock (with KS-DFT hybrids somewhat in between). Moreover, it has been repeatedly shown that the HF-like inter-atomic exchange contribution between bonded atoms is attractive and essentially responsible for the bonding. Ignoring this term would make most IQA inter-atomic energies of bonded atoms positive.

An obvious drawback of this strategy is that the total KS-DFT exchange–correlation energy is not recovered by the sum of all atomic and inter-atomic terms, i.e., the decomposition is not fully additive.

This strategy performed extremely well from a numerical point of view. Both the atomic and diatomic exchange energy components obtained with a local exchange functional (e.g., from a BLYP calculation) exhibited almost perfect correlation with the values obtained with the exact exchange formula (using the same KS-MOs and geometries).³² 

In 2016, Francisco et al. introduced an alternative strategy (henceforth, F-IQA) to decompose the KS-DFT exchange–correlation energy fulfilling Eq. (19).³³ The idea was again to use the exact exchange formula using the KS-MOs but incorporating properly defined atomic scaling factors to ensure additivity.

Hence, in the F-IQA method, the exact HF-exchange expression is always used to determine the atomic and inter-atomic exchange–correlation energies even in the case of purely local KS-DFT functionals.

The main purpose of this work is to assess the performance of the aforementioned KS-DFT IQA schemes. The atomic and inter-atomic exchange–correlation terms obtained for a molecular set using a local (GGA) functional, a hybrid (GGA) functional, and Hartree–Fock are compared. In addition, we also describe in detail a numerical procedure to improve the accuracy of the decomposition of the two-electron energy terms, which are the bottleneck of the IQA approaches. We show that in the context of overlapping atomic definitions, it is possible that the sum of all atomic and diatomic contributions to the two-electron energy (i.e., Coulomb and exact exchange if needed) exactly reproduces the analytical molecular value. We refer to this strategy as the (two-electron) zero-error scheme (ZES).

The two-electron contributions to the molecular energy, namely, Coulomb and Hartree–Fock exchange and correlation, in the case of wavefunction methods are both the bottleneck and the major source of the numerical error in the IQA energy decomposition schemes. They formally scale N⁶, N being the number of grid points, albeit efficient algorithms achieving N⁴ scaling have also been introduced in the QTAIM framework.³ In general, the numerical integration of one- or two-electron density functions over disjoint (e.g., QTAIM) atomic domains is more challenging than for fuzzy-atom (e.g., TFVC) domains. In the latter case, the atomic weight functions [Eq. (1)] tailor the molecular density functions. Then, the integrals are formally carried out over the entire space for which conventional and well-tested numerical integration schemes can be safely applied. Among them, the multicenter integration scheme introduced by Becke, which merges atom-centered spherical grids, is by far the most widely used.³⁴ 

In the first realization of the Hartree–Fock IQA decomposition in the general framework of fuzzy atoms, the numerical two-electron integrations required to decompose the Coulomb and exchange energy terms were carried out by using two sets of atom-centered grid points associated with the electron coordinates of electron 1 (r₁) and 2 (r₂).² It was shown that because of the $r12−1$ operator in the denominator, if exactly the same grid was used for both electron coordinates, all points where r₁ = r₂ (in fact, N grid points) had to be discarded to avoid the singularity. Thus, the overall accuracy of the integration is compromised, i.e., the sum of the one- and two-center terms compared to the corresponding molecular (analytical) value. In order to avoid this situation, the authors used two identical grids for both electrons but one of them rotated along the ϕ angle of the spherical coordinates to obtain the one-center terms. In this manner, sufficiently good accuracy was achieved using atom-centered grids consisting of 40 radial and 146 angular points. For the angular mesh, the grid associated with electron 2 was rotated 0.229 rad along ϕ.

This idea is somewhat reminiscent of the staggered mesh methods applied in solid-state calculations.^35–37 One could indeed use staggered atomic grids in the framework of Becke’s multicenter integration. This might increase the accuracy of the one-electron integrations, in particular. However, our present purpose is numerical two-electron integrations and, as we will see, the ZES scheme already affords the decomposition of the two-electron energy terms with full numerical accuracy.

We have observed that appropriate rotation angles depend mostly on the size of the angular grid. In addition, in a preliminary analysis, we explored the possibility of introducing a second rotation (with respect to the θ angle) to the grid for electron 2. No significant improvement was observed for a small set of molecular systems as compared to that obtained with a single rotation. Another observation was that the two-center two-electron terms are contributing to the numerical error several orders of magnitude less than the one-center ones while being rather unaffected by the rotation of the second grid.³⁸ 

However, the most relevant observation is depicted in Fig. 1. Here, the integration error in the two-electron B3LYP/cc-pVTZ energy of N₂ is shown with respect to the rotation angle of the second grid. The error is just defined as the difference between the sum of all one- and two-center contributions of the molecular two-electron energy and the exact (analytical) value, V_ee. One can see that if the atomic grid provided is sufficiently large, there is always a rotation angle for which the overall two-electron integration error vanishes. Thus, a 30 × 74 (radial and angular points, respectively) atomic grid is clearly insufficient, and with the optimum rotation angle, the error is still of ∼15 kcal/mol. However, using a 40 × 146 atomic grid, the integration error first crosses the zero-error at around 0.18 rad and then again at 0.217 rad, close to the default rotation angle used in Ref. 2. By using larger atomic grids, the zero-error line is crossed at smaller rotation angles. In addition, several poles are observed in the curves due to near singularities in $r12−1$⁠. Noteworthy, using sufficiently large atomic grids, such as 70 × 434 or 150 × 590, the shape of the curves is strikingly similar irrespective of the molecule. That is, the rotation angle that produces a zero-error in the two-electron energy lies within an extremely narrow range. This is illustrated in Fig. 2, where a 150 × 590 atomic grid for a set of ten small molecules was used. With this, we observed that the optimum rotation angle is within the 0.169–0.170 (in rad) range in all cases.

It can be readily seen that the one- and two-center terms thus defined exactly reproduce the total two-electron energy.

Ideally, the applied rotation angles in steps 1 and 2 should be previously optimized to provide small deviations with respect to the exact two-electron energy and, most importantly, of opposite sign. In that case, the γ value lies within the [0, 1] range and an actual damping (interpolation) is performed between the two estimates of V_ee. As such, a linear behavior of V_ee with the rotation angle of the second grid is implicitly assumed. As mentioned before, we have found that a combination of 150 × 590 atomic grids and a rotation of 0.169 rad typically overestimates V_ee by 1–5 kcal/mol, while using a rotation of 0.170 rad in the second step leads to a somewhat larger underestimation of V_ee. In Table I, we gathered the integration errors in V_ee for a set of 31 molecules obtained at the B3LYP/cc-pVTZ level of theory (for further details, see Sec. IV). It can be seen that the two-electron integration errors in step 1 are already rather small and negative with the only exception of H₂ for which the error is merely 0.1 kcal/mol. For most applications, these errors might be acceptable. The errors associated with the second step are somewhat larger but positive in all cases so that the γ values that afford the exact decomposition are within zero and one.

It is important to stress that with such (two-electron) zero-error scheme the two-center terms are evaluated only once in the first step. We have observed that their value is rather unaffected by a rotation of the second grid and exhibits integration errors 2–3 orders of magnitude smaller than those of the one-center terms.³⁸ This is because there are no near singularities caused by small $r12−1$ values and also because their contribution to the total two-electron energy is much smaller than that of the one-center terms, specially if heavy atoms are involved. Thus, it is only the much larger one-center terms that are slightly modulated to yield an overall exact decomposition of the two-electron energy.

Of course, the zero-error scheme could also be applied independently to each of the contributions to V_ee, namely, Coulomb and exact-exchange (and correlation in case of correlated wavefunction methods), provided their exact value is known beforehand. However, we find that the behavior of V_ee in the vicinity of the optimum rotation angles is very similar for HF, B3LYP, and BP86, as illustrated by Fig. 3. Since the overall errors are very small anyway, it is more efficient to apply the zero-error scheme only once to reproduce the exact V_ee.

In this vein, Table S1 of the supplementary material gathers the γ values for the molecular test set at the HF, B3LYP, and BP86 levels of theory. The integration errors have been omitted for clarity, all of them being lower than 0.5 kcal/mol (in absolute value), as proceeding from the one-electron integrations. For the later, the same 150 × 590 atomic grid was used, the main source of the numerical error being the atomic kinetic energy contributions due to the oscillatory character of the kinetic energy density near the nuclei. The C₆H₆ molecule is the worst case with errors of up to 0.4 kcal/mol. Of course, the one-electron grid can still be further improved to decrease the residual integration error if necessary. More importantly, using the aforementioned 0.169 and 0.170 rad rotation angles for the two steps of the two-electron zero-error scheme resulted in the desired interpolation in almost all cases, disregarding whether V_ee contains any exact-exchange contribution.

Summarizing this section, the two-electron zero-error scheme smoothly worked for the molecular systems tested, obtaining most of the γ values within the [0, 1] range using a system-independent fixed setup. In very few cases, the corrected two-electron energy terms were obtained by extrapolation but with γ values still close to one. Alternatively, in the case the γ value would be far off the [0, 1] range, one could consider performing second grid rotations iteratively until the desired numerical conditions are fulfilled. It is worth to mention that the robustness of the zero-error scheme had already been put into stringent test. Some of us recently showed that the elements of the α tensor can be expressed through the second derivative of a zero-th order field-dependent energy so that an energy decomposition of the latter readily affords the decomposition of α.²⁹ Using very large atomic grids for the one-electron energy terms and the zero-error scheme for the two-electron part was absolutely critical to obtain numerically converged atomic contributions to α. Overall, the zero-error scheme appears to be a promising strategy for obtaining accurate IQA energy decomposition terms for large molecular systems while keeping an affordable computational cost.

The atomic and inter-atomic exchange–correlation energy terms obtained using the aforementioned F-IQA and SM-IQA schemes for the molecular set (see Sec. IV) are compiled in Tables II–IV. We have considered the systems described at the HF, B3LYP, and BP86 levels of theory at their own geometry optimized structures. We focus exclusively on the exchange–correlation (xc) terms as the remaining ones are exactly the same with both KS-DFT IQA approaches. Regarding the one-center xc terms, we discuss separately the results for the hydrogen atoms as in our previous study we observed significant differences between HF and pure KS-DFT due to the absence of core electrons.³² In addition, the atomic xc energies for H are of the same order or magnitude as the inter-atomic xc terms of bonded atoms, while the corresponding values for heavier atoms are one or two orders of magnitude larger.

Let us start with the analysis of the inter-atomic xc components with up to 43 values compiled in Table II. The first observation is that all xc (exchange-only in case of HF) inter-atomic contributions of chemically bonded atoms are negative as it is well known for IQA decompositions. Their magnitude (in absolute value) is deeply connected with the covalent bond order so that bonds with higher multiplicity tend to exhibit larger (more negative) inter-atomic xc contributions. Thus, the larger values are obtained for N₂, NO⁺, and the C–C triple bond in C₂H₂.

The inter-atomic xc terms obtained with B3LYP and BP86 functionals correlate extremely well with the HF values using either the F-IQA and SM-IQA formulations, as shown in Fig. S1 of the supplementary material (worst case exhibits R² = 0.98). It is more interesting to focus on how the contributions differ from each other in each case. With the F-IQA formulation, the inter-atomic xc components for B3LYP are systematically more negative than the HF ones with only two exceptions (F₂ and H₂S). In both cases, this discrepancy can be attributed to significant differences in the wavefunction itself (F₂ is poorly described at the HF level, and the shape of the atomic boundaries in H₂S was already found to change significantly from one method to another²⁰). This observation is a direct consequence of the scaling factors used in Eq. (27). The total xc energy in B3LYP contains both exchange and correlation contributions, whereas in HF it only contains exchange. Since the total xc value is more negative, the scaling factors are greater than one, and hence, all terms (both atomic and inter-atomic) become more negative. Going to BP86, the picture is very similar. In fact, the mean unsigned deviation (MUD) between the BP86 and B3LYP xc components is merely 2.2 kcal/mol, whereas the respective MUD values with respect to HF are 14.5 and 14.2 kcal/mol, respectively.

With the SM-IQA formulation, the trends are different. Most of the B3LYP inter-atomic xc contributions (32 out of 43) are less negative than the HF ones. The differences are significant in the C–C and C–H bonds of the alkane series. However, the opposite trend is observed mainly for the systems with large (in absolute value) inter-atomic xc contributions, namely, triple bonds, such as in acetylene, CO, N₂, or NO⁺ (see Fig. S1). Similarly to F-IQA, the MUD between BP86 and B3LYP SM-IQA values is very small (3.5 kcal/mol).

It is worth to note that the general effect of electron correlation in wavefunction theory is precisely the weakening of the covalent bonds. This is readily observed in equilibrium geometries, where bond lengths tend to increase upon inclusion of electron correlation. It is also well-known that electronic bonding indicators, such as bond orders, also decrease. For instance, the bond order in N₂ goes from 3.04 at the HF level to 2.83 for a CISD WF using the Ángyán–Mayer formulation, which only involves the first-order density matrix (i.e., the exchange density). Including the contribution from the cumulant of the second-order reduced density matrix (the so-called delocalization index³⁹) further decreases the value to 2.20.⁴⁰ In the seminal work of Blanco et al. about the IQA approach for correlated wavefunctions, the authors showed that the IQA decomposition of the cumulant of the RDM2 (electron correlation contribution) results in positive inter-atomic contributions for systems such as H₂, N₂, or the O–H bond in H₂O. Popelier et al. performed a more systematic study of the role of electron correlation in the IQA-MP2 decomposition.⁴¹ The authors found that the inter-atomic correlation energy contributions are usually positive for covalent bonds. In addition, the stronger the bond the larger the inter-atomic correlation contribution. Exception to the rule were hydrides and weak interactions, exhibiting much smaller and negative contributions. Of course, the total electron correlation energy contribution is negative; therefore, the atomic correlation terms are large and negative to compensate.⁴² In the F-IQA formulation, this effect should be captured by the scaling factors. As shown on Table III, the atomic xc energies from both B3LYP and BP86 are systematically more negative than the HF ones, which means that the net effect of correlation is as anticipated. Table IV gathers the atomic xc energies for hydrogen atoms, which are much smaller than for second and third period atoms (in fact, they are of the same order of magnitude as the inter-atomic terms for bonded atoms). The differences between the HF and F-IQA values are very small. Except in the aforementioned case of H₂S, where the shape of the atomic domains changes significantly going from HF to a DFT density, the MUD between the HF and F-IQA atomic contributions for H atoms is merely 7–8 kcal/mol. More importantly, the atomic xc values with BP86 are less negative than the HF values in several cases. In contrast, with the SM-IQA formulation, the atomic xc contributions are systematically more negative than the HF ones; in the case of the H atoms by ∼20 kcal/mol on average. This means that the observed trends of the atomic correlation contributions in IQA-MP2⁴² are better captured by the SM-IQA formulation at least for the set of studied systems.

It is fair to note that the numbers discussed so far have been obtained at the corresponding stationary points on the potential energy surface of the given level of theory. Thus, several factors influence the observed differences between HF and BP86 or B3LYP one- and two-center xc contributions, namely, the way the total exchange (correlation) energy is expressed, the shape of the MOs, and the geometry. To explore the latter effect, we have studied in more detail F₂ for which a deviation from the general trend was observed (i.e., the F-IQA inter-atomic xc values are less negative, and the SM-IQA values are more negative than the HF ones). In Fig. 4, we show the evolution of the inter-atomic xc values with the inter-atomic distance for HF, BP86, and B3LYP wavefunctions. Notice that the HF optimized geometry (1.328 Å) is significantly shorter than the experimental value (1.412 Å), while BP86 slightly overestimates the distance. The inter-atomic xc energies for KS-DFT using either F-IQA or SM-IQA formulations are more negative than the HF ones for all inter-atomic distances. An important observation is that the differences between all schemes are roughly constant along the (rather short) inter-atomic distance profile, while the shorter the inter-atomic distance, the more negative the inter-atomic xc value. Consequently, the IQA-DFT xc values at their equilibrium distances may lie above the HF values at the compressed HF equilibrium distance even if the IQA-DFT xc values are systematically more negative than the HF along the profile. What makes the bond in F₂ different from the other cases is that the SM-IQA values lie well below the F-IQA ones. A plausible explanation might be the unexpectedly large value of the bond order around 1.4 for the considered distances.

As a final illustrative example, we applied the IQA decomposition along the energy profile [i.e., intrinsic reaction coordinate (IRC) obtained at each level of theory] of the Diels–Alder reaction between 1,3-butadiene and ethylene. As shown in Fig. S2, Hartree–Fock largely overestimates the barrier height (∼45 kcal/mol), whereas the pure GGA BP86 functional underestimates it (∼15 kcal/mol). Hybrid functionals, such as B3LYP, do a particularly good job in this case, predicting a barrier height (∼25 kcal/mol) close to the experimental estimate of 27.5 kcal/mol.⁴³ In any case, this reaction profile poses a good example where the three methods considered in this work exhibit quantitative different behaviors. The IQA decomposition along the IRC in each case can provide insight onto the origin of such differences. On the one hand, one could analyze each atomic and diatomic contribution separately. Because of the large number of energy contributions (IQA descriptors) obtained, this usually requires the assistance of some sort of statistical analysis, such as the Popelier’s REG scheme.⁴⁴ An alternative method involves examining the change in total energy compared to the energy of the isolated reactants (diene and dienophile) and organizing the IQA terms according to the corresponding fragments. This results in fragment deformation energies (the IQA energy of the fragment relative to the energy of the free isolated fragment) and inter-fragment contributions. This also permits to establish proper analogy with the so-called activation strain model (ASM),^45,46 where the energy along a reaction profile is decomposed into strain energy (computed as the difference between the energies of the isolated reacting fragments at the constrained geometry of each point of the profile and that of the isolated fragments) and interaction energy. The former is destabilizing by definition, while the latter ultimately accounts for the formation of the products. Casals-Sainz and co-workers⁴⁷ recently showed how the endo/exo preferences of Diels–Alder reactions can be explained in terms of IQA descriptors without recurring to additional reference states, such as in the ASM model.

In the present work, we focus on the origin of the above-mentioned differences in the activation barriers obtained with the three electronic structure methods. Figure 5 illustrates to what extent these differences are translated into the intra-fragment (namely, the dienophile (F₁) and the diene (F₂) and inter-fragment F₁–F₂ interaction terms from the respective IQA analyses. In line with the picture obtained with the ASM model, the fragment deformation energies dominate the energy profile starting from the reactant complex. Both F₁ and F₂ deformation energies monotonically increase as the reaction proceeds until well passed the transition state (TS). The deformation of the diene is slightly larger than that of the dienophile before reaching the TS. Surprisingly, according to the F-IQA results, the fragment deformation is initially smaller for the HF profile during the early stages of the reaction but then abruptly increases and surpasses the KS-DFT values before reaching the TS. As the reaction progress, the differences between the three methods become minor yet follow the same trend as the total energy along the profile (i.e., HF > B3LYP > BP86). The inter-fragment interaction is always negative and steadily increases (in absolute values) until the formation of products (where the fragments are covalently bound). Although the differences in this term for the three methods are minimal, they increase after the TS. The analysis suggests that the reason behind the overestimation of the reaction barrier by HF is the increased fragment deformation, which is not counterbalanced by a significantly higher inter-fragment stabilization.

When comparing the behavior of the two IQA-DFT formulations, it appears that the SM-IQA results exhibit more noticeable discrepancies among the three levels of theory (perhaps somewhat exaggerated taking into account the relatively small overall energy differences between the methods in Fig. S2). Specifically, the fragment deformation of both the diene and dienophile is considerably greater for HF throughout the entire reaction profile. Additionally, the inter-fragment stabilization is also significantly larger for HF than for the DFT methods. Figure 5 illustrates that the trend of HF > B3LYP > BP86 is maintained for all terms throughout the entire reaction profile.

A new scheme to eliminate the numerical error of the sum of two-electron energy contributions (i.e., Coulomb and exact exchange) has been introduced. It is readily applicable in the framework of overlapping atoms such as TFVC or Hirshfeld approaches. Such a two-electron zero-error scheme provided robust numerical integrations for both equilibrium structures and along the potential energy surface.

A critical comparison of the performance of two fully additive approaches for the IQA decomposition of the KS-DFT energy (F-IQA and SM-IQA) has also been performed for a wide molecular test set and along the reaction coordinate of a Diels–Alder reaction. The atomic and diatomic exchange–correlation energy components obtained with both approaches are in very good agreement and also exhibit excellent correlation with the Hartree–Fock results. As a general trend, the SM-IQA diatomic xc components tend to be less negative than the Hartree–Fock ones in accordance with the known effect of electron correlation on covalent bonds.

All ab initio calculations have been performed using the Gaussian09 package⁴⁸ at the Hartree–Fock (HF) and at the BP86^49,50 and B3LYP^51,52 Kohn–Sham DFT levels of theory coupled with the cc-pVTZ full electron basis set.⁵³ Molecular (equilibrium) structures have been obtained for each level of theory and confirmed by vibrational frequency analysis (no negative frequencies).

Real-space energy decompositions have been performed with the APOST-3D program⁵⁴ using the topological fuzzy Voronoi cell atomic definition²⁰ for both the F-³³ and SM-IQA³² approximations of the inter-atomic exchange–correlation energy term. For the production results, both one- and two-electron numerical integrals were evaluated using 150 radial and 590 angular points distributed according to Lebedev–Laikov.⁵⁵ The zero-error scheme was applied in all cases for the two-electron energy contributions.

The test set of 31 molecules consists of carbon and sulfur oxides, particularly SO₂, SO₃, CO, and CO₂; the hydrocarbon series of C₂H₆, C₆H₆, C₂H₄, and C₂H₂, a set of second and third row hydrides with general formula XH_n, where X goes from Li to Cl and other neutral and charged inter-atomic molecules exhibiting different bond multiplicities, e.g., CN⁻, N₂, or F₂.

Geometries and wavefunctions of the prototypical Diels–Alder reaction between 1,3-butadiene and ethylene were obtained at the HF, BP86, and B3LYP levels by means of IRC calculations starting from the optimized transition states. IRC step numbers were increased to ensure convergence to both reactants and products at each level of theory. The Hessian matrix was recomputed at each step of the IRC to ensure smooth potential energy curves.

Damping factors [see Eq. (31)] for the dataset (Tables S1 and S2), correlations between the HF and KS-DFT atomic and diatomic xc components (Fig. S1), and energies along the intrinsic reaction coordinate for the Diels–Alder reaction for HF, B3LYP, and BP86 methods (Fig. S2).

Financial support has been furnished by the Ministerio de Ciencia, Innovación y Universidades (MCIU), Grant No. PGC2018-098212-B-C22. M.G. acknowledges the Generalitat de Catalunya and Fons Social Europeu for the predoctoral fellowship (Grant No. 2018 FI_B 01120).

The authors have no conflicts to disclose.

Martí Gimferrer: Data curation (equal); Formal analysis (equal); Software (equal); Writing – original draft (equal). Pedro Salvador: Conceptualization (equal); Funding acquisition (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material. Extra data are available from the corresponding author upon reasonable request.