Merging the Energy Decomposition Analysis with the Interacting Quantum Atoms Approach

Energy decomposition analysis (EDA) is a well-established approach to dissect the interaction energy into chemically sound components. Despite the inherent requirement of reference states has been a long-standing object of debate, the direct relation with the molecular orbital analysis helps in building up predictive models. The alternative molecular energy decomposition schemes that decompose the total energy into atomic and diatomic contributions, such as the interacting quantum atoms (IQA), has no external reference requirements and also the intra- and intermolecular interactions are treated on equal footing. However, a connection with heuristic chemical models are limited, bringing about a somewhat narrower predictive power. While efforts to reconcile the bonding picture obtained by both methodologies have been discussed in the past, a synergic combination of them has not been tackled yet. Herein, we present the use of IQA decomposition of the individual terms arising from the EDA in the context of intermolecular interactions, henceforth EDA–IQA. The method is applied to a molecular set covering a wide range of interaction types, including hydrogen bonding, charge–dipole, π–π and halogen interactions. We find that the electrostatic energy from EDA, entirely seen as intermolecular, leads to meaningful and non-negligible intra-fragment contributions upon IQA decomposition, originated from charge penetration. EDA–IQA also affords the decomposition of the Pauli repulsion term into intra- and inter-fragment contributions. The intra-fragment term is destabilizing, particularly for the moieties that are net acceptors of charge, while the inter-fragment Pauli term is actually stabilizing. In the case of the orbital interaction term, the sign and magnitude of the intra-fragment contribution at equilibrium geometries is largely driven by the amount of charge transfer, while the inter-fragment contribution is clearly stabilizing. EDA–IQA terms show a smooth behavior along the intermolecular dissociation path of selected systems. The new EDA–IQA methodology provides a richer energy decomposition scheme that aims at bridging the gap between the two main distinct real-space and Hilbert-space methodologies. Via this approach, the partitioning can be used directionally on all the EDA terms aiding in identifying the causal effects on geometries and/or reactivity.


■ INTRODUCTION
Understanding and accurately assessing intra-and intermolecular interactions is one of the main challenges in chemistry. In fact, the rational design of molecular systems consists of unravelling the physical origin of a particular chemical interaction/bond, often inaccessible directly from experiments. However, a common drawback is the absence of an exact quantum-mechanical operator that directly describes the chemical bond, giving raise to different approaches.
Among the number of developments, some methods focus on the analysis of the electron density in a system (AB) by comparing it with that from the composing fragments (A and B). The concept of deformation density 1 is commonly invoked in methods such as Voronoi deformation density charges 2 or charge displacement analysis, 3,4 for instance. The electron density of the AB system is a crucial component in the quantum theory of atoms in molecules 5 (QTAIM). By analyzing its topology, QTAIM provides a plethora of descriptors that can be utilized to classify various intra-and intermolecular interactions.
A better option is to focus directly on the energetics of bond formation and intermolecular interactions. Modern electronic structure methods are able to predict accurate formation energies, but the value itself bears little chemical significance. Energy decomposition schemes aim at decomposing the molecular (or formation) energy into physicochemical meaningful terms, to shed light into the nature of the chemical bonding. By comprehending the individual contributions to the overall energy, it becomes possible to rationally design molecular systems with desired properties, leading to a more predictive approach in molecular design.
One of the most widely used methodologies is the Ziegler− Rauk energy decomposition analysis (EDA), 6 derived from the pioneering work of Kitaura and Morokuma. 7 It considers the formal molecule formation AB (henceforth complex) from fragments A and B (atoms or molecular fragments). The overall stabilization energy (without basis set superposition error correction) reads as where E(X Y ) refers to the energy of the subsystem X at the optimized geometry of Y. Thus, ΔE stab is the energy of formation of the system AB from the isolated fragments in their ground states A and B. Within the EDA formalism, ΔE stab (eq 1) is decomposed as follows ) being ΔE int and ΔE prep the so-called interaction and preparation energy terms, respectively, which are defined as Here, E(A 0,AB ) represents the energy of fragment A computed at the optimized geometry of the complex (superindex AB) with a given electronic configuration (A 0 ), which may not correspond to that of the ground state for the isolated fragment. Defined as such, the preparation energy accounts for both the geometrical distortion of the fragments upon formation of the complex and the promotion energy from the electronic ground state to the chosen electronic configurations A 0 and B 0 . One often refers to strain energy when it only involves the geometrical deformation. Furthermore, it is necessarily positive (repulsive or destabilizing) because E(A A ) and E(B B ) are, by def inition, the ground state energies of the isolated fragments from both the electronic and geometric perspective. Instead, ΔE int will be negative (attractive or stabilizing) if the interaction between the fragments A and B while forming complex AB is favorable. Importantly, both the interaction and preparation energies depend on the choice for the states A 0 and B 0 as indicated in eqs 3 and 4, being crucial its appropriate selection (see below for further details). By introducing additional intermediate (pseudo)states built up at the optimized geometry of the complex, the interaction energy is further decomposed. Firstly, one considers a pseudostate of complex AB formed by the superposition of the undeformed (frozen) densities associated to the fragments in states A 0 and B 0 , namely (A 0 ∪ B 0 ), with its associated electronic energy E(A 0 ∪ B 0 ). We refer to (A 0 ∪ B 0 ) as a pseudostate because it does not have a well-defined antisymmetric wavefunction associated to it. 6,8 The energy difference with respect to the deformed fragments read as  ) where the AB superindex has been omitted for clarity. The term ΔE elec [A 0 , B 0 ] accounts for the electrostatic interaction of the frozen electron density of fragment A with the nuclei of fragment B and vice versa (attractive), the Coulombic repulsion of the frozen electron densities of A and B and the nuclear repulsion between A and B. i k j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z (8) where the molecular electrostatic potential (MEP) of fragment A in state A 0 , V A (r), reads as The potential originated from charge clouds is smaller than the one from point charges (nuclei) so that for neutral species the MEP of the fragments afford a favorable interaction that, at chemically relevant distances, overcomes the nuclear repulsion term. 9 In Kohn−Sham density functional theory (KS-DFT) there is an additional contribution from the exchange− correlation functional, 8,9 which is absent in wavefunction theory.
In a subsequent step, an intermediate state (A 0 B 0 ) is formed by Loẅdin orthogonalizing the occupied molecular orbitals (MOs) of the states (A 0 ) and (B 0 ), in order to build a proper antisymmetrized wavefunction. Orthogonalization is required as MOs belonging to different fragments are not orthogonal (in principle one could build a Slater determinant with nonorthogonal MOs but then the expectation value of the energy takes a much complicated form). The Loẅdin orthogonalization procedure does not induce charge transfer between the fragments, as the Hilbert-space based electron numbers of the interacting fragments are conserved. This will not be the case when applying a real-space analysis, as will be discussed later.
The energy difference between this intermediate step and that of the previous pseudostate reads as 10) which upon combination of eqs 6 and 10 leads to the so-called Pauli repulsion term (ΔE Pauli ) 8 MOs, and it is a well-defined quantity in the sense that involves properly antisymmetrized states 12) It is worth to note that E A B ( ) 0 0 0 is never considered, as is not evaluated explicitly. Once the electrostatic contribution is calculated using eq 7 or 8, the Pauli repulsion term is readily obtained from eq 12.
In the last step, the MOs of the complex are allowed to relax to the ground state of the complex. The energy lowering accompanying this process leads to the so-called orbital interaction term (ΔE orb ), that is necessarily negative (any intermediate state must be higher in energy than the ground state) All the steps along the EDA process are generally illustrated as in Scheme 1.
Localized orbitals can also be introduced in this type energy decomposition schemes. In the absolutely localized molecular orbitals EDA (ALMO-EDA), 10−13 the interaction energy is further decomposed into a frozen-density, polarization and charge-transfer terms by making use of (variationally optimized) block-localized orbitals, and explicitly avoiding any reference to intermediate pseudostates. Natural EDA 14 (NEDA) makes use of the well-known natural bond orbitals (NBO). 15−18 In the symmetry adapted perturbation theory (SAPT) schemes the interaction energy is perturbationally computed, thus avoiding the supermolecular approach, and decomposed into physically meaningful terms. 19 For further details about these EDA-like methodologies we guide the reader to refs 20 and 21.
Alternatively, the total energy of any molecular system can also be decomposed into intra-and inter-atomic contributions. Such decomposition does not require external references or predefined fragments, and treat intra-and inter-molecular interactions (covalent and non-covalent) on equal footing. Grouping the one-and two-center terms into intra-and interfragment contributions is only optional, but helpful in the case of dealing with intermolecular interactions between well-defined subsystems. The grouping of specific domains within the interacting fragments also allows the identification of main contributors and their mutual interactions. In Mayer's Chemical Hamiltonian Approach, atomic projector operators are used to decompose the Hamiltonian into one-and two-center terms. 22 Further developments in the Hilbert-space ultimately lead to the chemical energy component analysis. 23 Considering instead a decomposition of the real-space, the one-and two-electron contributions to the total energy readily afford one-and twocenter terms that exactly decompose (up to numerical accuracy) the molecular energy. Such methodologies rely on the identification of the atom within the molecule (AIM). Salvador and Mayer first decomposed the Hartree−Fock energy in the framework of QTAIM, 19 paving the way for the nowadays known as the interacting quantum atoms approaches. 24−31 In real-space analysis, a given quantity, F 1 , expressed in terms of a one-electron density function, f(r 1 ), is readily decomposed into one-center (atomic or fragment) contributions by integrating over the respective domains. Similarly, twoelectron quantities decompose into both one-and two-center components r r r r r r r r r r r r It is worth noting that real-space analysis is not restricted to non-overlapping disjoint domains such as those of QTAIM, where each atom is identified by its nucleus and its atomic basin. The AIM may be more generally represented by continuous atomic weight functions w A (r) ≥ 0 fulfilling Σ Α w A (r) = 1 so that the integration of molecular density functions over the atomic domains are effectively replaced by integrations over the whole real-space of atomic/diatomic effective density functions r r r r r r r r r r r r r r r r r r r r r Such atomic weight functions can be derived from a variety of Hirshfeld-type approaches 32 or even mathematical constructs borrowing elements of QTAIM theory. 33 Whether the AIM are allowed to overlap or not might be to some extent matter of taste. Using one or another AIM only has an effect on the actual numerical values obtained for the terms obtained by the IQA decomposition, but not on their definition and physical meaning.
Since the Born−Oppenheimer energy is entirely written in terms of one-and two-electron energy density functions, IQA naturally affords the decomposition of the molecular energy of a complex AB into atomic and diatomic contributions as where the ε i and ε ij terms account for the static net atomic and pairwise atomic interaction energies, respectively. The atomic and diatomic terms can be further grouped according to the composing fragments A and B, so that the total energy of the complex can be simply expressed as and Each of the intra-and inter-fragment energy term is built upon the physical components of the electronic energy, i.e. kinetic, electron nuclear attraction, nuclear repulsion and electronic repulsion. The latter may be further decomposed into the usual Coulomb, exchange and correlation contributions. One peculiarity of the IQA-type approaches is that the actual formulation depends upon the particular electronic structure method that is used to compute the molecular energies in the first place. Appropriate formulations have been developed for Hartree−Fock 24,26,34 and correlated methods (including CASSCF and CI, 27 MP2 28,35 and Coupled Cluster 29,30,36 ). Non-perturbative approaches explicitly require the second order reduced density matrix, which is not available in most electronic structure codes. Curiously, the KS-DFT case is the most problematic one because of the exchange−correlation energy nature. Within wavefunction theory framework, the exchange energy is expressed as a two-electron non-local contribution, that naturally decomposes into both one-and two-center terms. The later are essential to account for the stability of the diatomic bonding interactions. 21 In KS-DFT, the exchange−correlation energy is essentially written in terms of the exchange energy density as ) so that the straightforward real-space decomposition only affords one-center (atomic) terms. Different approaches have been introduced to (approximately) recover the chemically meaningful diatomic components from V xc . 25,31 Both EDA and IQA methodologies independently have been extensively used in the literature to gain deeper insight into the nature of the chemical bond and to characterize intra-and intermolecular interactions, allowing to understand and improve chemical reactivity, shedding light to the chemical-bonding picture of non-trivial systems and even most recently suggesting a new type of bond. 37 43,44 The authors decomposed the interaction energy between the two monomers A and B (E int AB ) into the sum of classical electrostatics (V cl AB ) and exchange−correlation (V xc AB ) contributions. They observed that the interaction was governed by the exchange−correlation, thus highlighting the importance of the covalent picture. On the other hand, the deformation energy of the proton acceptor moieties correlated well with the intermolecular charge transfer and classical electrostatic energy derived from IQA. Furthermore, by making use of the fragment's promolecular, polarized (by locating point charges) and fully relaxed densities, they observed that in weakly-bonded (almost non-overlapping) systems the quantities defined by other energy decomposition schemes, i.e. SAPT, KM, EDA and specially NEDA, can be obtained to a good approximation from the inter-fragment (AB) IQA terms. For instance, the electrostatic energy from NEDA was found to be roughly equivalent to the total inter-fragment interaction from IQA. 44 Pendaś et al. also critically analyzed the concept of steric repulsion from an IQA perspective. 45 The authors argued that Pauli repulsion is inherently dependent on the fragment's reference states in EDA. They applied IQA to decompose the Hartree−Fock interaction energy into fragment's deformation and inter-fragment interactions where the latter is further decomposed into its classical electrostatic and exchange−(correlation) contributions. The authors concluded that the Pauli repulsion is readily captured in the increase of the fragment's deformation energies of the intermediate (properly antisymmetrized) states. In the case of rotational barriers, the hyperconjugative effects are captured by the inter-fragment exchange contribution, enhanced due to electron delocalization. All in all, they show a certain degree of correspondence between EDA or NBO descriptors and those steaming from IQA. More recently, Racioppi et al. walked a reverse path. Instead of recovering EDA descriptors from IQA, they rearranged the EDA contributions to match those of IQA analysis. 46 In particular, in their pseudo-IQA energy decomposition the EDA contributions of Pauli repulsion, orbital interaction and electrostatic to the interaction energy are regrouped into overall variations of the kinetic, classical electrostatic and exchange−correlation contributions The same terms can be obtained by considering the usual reference-state IQA, which is based on decomposing the binding energy between two fragments A and B (ΔE bind IQA ) by subtracting the IQA terms from the fully relaxed complex's state from those obtained for the isolated fragments at the complex geometry. The authors showed excellent agreement between the like terms of both schemes in illustrative hydrogen bond and donor− acceptor interactions. 46 In this work, we pursue a different path, namely to enrich the conventional EDA approach by applying an IQA decomposition to each of the EDA terms of the interaction energy. Thus, in the EDA-IQA scheme we introduce herein, the electrostatic, Pauli repulsion and orbital interaction energy terms are decomposed into intra-and inter-atomic contributions, that can be further grouped into intra-and inter-fragment contributions.

■ THEORY
Let us consider again the formation of the complex AB from fragments A and B. The application of eqs 14−16 to the complex's final ground state (AB) readily affords the real-space decomposition of the interaction energy into intra-and interfragment terms, namely where For clarity, in the previous equation and henceforth we omit the explicit dependence of the EDA term on the reference states (A 0 and B 0 ). The Δε int,A and Δε int,B account for the energy gain/ loss by the fragments when going from their isolated reference state to their ef fective state within the final complex. It is worth to note that in the context of real-space analysis, these contributions do not only originate from changes in the MOs upon complex formation, but also by the fact that the fragments share the physical space once the complex is formed (in intermolecular interactions the second effect should be dominant). In ref 44 the authors refer to these terms as fragment's electronic deformation energies. We will adopt here this nomenclature, so that Δε int,A ≡ Δε def.el,A and Δε int,B ≡ Δε def.el,B .
On the other hand, the term Δε int,AB describes the energy gain upon complex formation that can be purely ascribed to interfragment interactions. The net interaction energy is thus seen as a balance between the prize the fragments must pay to share the physical space and be electronically prepared, and the gain originating from the new interactions that were absent before the complex's formation.
In a similar fashion, by applying again eqs 14−16 to the complex's intermediate state (A 0 B 0 ) one can also obtain an analogous decomposition of the orbital interaction EDA term, namely The intra-fragment terms account for the net energy gain/loss upon relaxing the wavefunction from the intermediate state to the ground state of the AB complex. This relaxation comes with a change in the electron density. If the underlying AIM definition depends upon this scalar (e.g. QTAIM, TFVC or iterative Hirshfeld approaches), these terms contain also a contribution from the change on the boundaries of physical space going from AB to A 0 B 0 . The latter could be removed by using the same AIM definition for states AB and A 0 B 0 . In the QTAIM context that means integrating the density functions of state A 0 B 0 on the atomic basins obtained from the AB state. In the case of overlapping AIM schemes, it implies using the same atomic weight functions throughout. Such strategies have been already used in the context of QTAIM and fuzzy atoms in similar contexts. 44,47,48 In the present case, since it is actually impossible by construction to use the same AIM definition for the complex and the isolated fragments, we opt for using the AIM definition derived from each state.
The IQA decomposition of state (A 0 B 0 ) readily affords an analogous decomposition of ΔE Pauli + ΔE elec , by taking the isolated fragment states A 0 and B 0 as reference. On the other hand, since each term in ΔE elec involves the electron density and/or potential from different fragments (see eq 7), one may argue that this term is entirely of intermolecular nature. In that case, ΔE elec would contribute solely to the inter-fragment term, and consequently one would have the following decomposition for ΔE Pauli However, such scheme is not satisfactory neither numerically nor conceptually. The main concern is that Δε Pauli,AB thus defined mixes up real-space and Hilbert-space quantities, while in this case they behave quite differently. Indeed, as mentioned above, there is no net charge-transfer between fragments A and B when building the intermediate state A 0 B 0 according to Hilbertspace analysis (e.g. Mulliken and Loẅdin populations add up to the number of electrons of each fragment). This is not the case when performing a real-space analysis (using disjoint or fuzzy domains), again because the fragments within the complex share the physical space Hence, the frozen density of isolated fragment A when brought to the complex's geometry does not entirely belong to fragment A, and similarly for fragment B. This influences the numerical values obtained using eq 27 and, for consistency, this effect should be also taken into account when applying the realspace analysis to the other EDA terms, and in particular to ΔE elec . One should essentially ignore the original allegiance of the fragment's frozen densities and treat the integrand in the exactly same manner as one does it with the electron-nuclear and the Coulombic contributions to the energy in the conventional IQA scheme, namely Here, we introduce the fragment's net electrostatic potentials V A net (r) and V B net (r). They are different from the electrostatic potentials V A (r) and V B (r) of eq 9, because in the electronic term the integration is carried out within the fragment's domain. Since only part of the fragment density is used,

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article Figure 1 depicts the topology of the net electrostatic potentials and fragment's densities along the inter-fragment distance for NaCl and water dimer. In the vicinity of the nuclei, the electrostatic potential is large and positive, as the electronic part cannot compensate for the nuclear term. At longer distances from the nucleus, the net potential slowly decays to zero. However, in the case of a strong acceptor or an anionic fragment, since there is an excess of electron charge compared to the nuclear one, the net potential becomes negative, and tends to zero from below (see light orange curves in Figure 1). This effect is much more pronounced when the donor of charge is anionic (Cl − vs H 2 O). On the other hand, there is a fraction of electron density of B that penetrates into A (i.e. w A (r)ρ B 0 (r)) and vice versa. It corresponds to the dark orange and dark blue curves in Figure 1. As expected, the density of the charge donor penetrates more and deeper into the acceptor region than the other way around. The interaction of that density from B with the net potential of the acceptor A results in the electrostatic contribution assigned to A, Δε elec,A . It corresponds to the integration of the grey curve in Figure 1b. This term is negative for the acceptor (notice the negative sign on the r.h.s. of eq 30) and can be significant if the density of the donor is able to penetrate deep into the acceptor's domain. However, in the case of the donor, the net potential can be negative in the region where it interacts with the density penetrating from the acceptor    The numerical value of this contribution will account to which extent the net potential of fragment A penetrates into fragment B and vice versa, compensated by the point charge nuclear repulsions. As we will see, this term can be positive of negative. In any case, one can readily see that It is fair to note that Jimeńez-Gravalos and Suaŕez recently achieved a similar decomposition of the electrostatic interaction in the QTAIM framework for a different purpose. 49 They did not explicitly considered fragment's electrostatic potentials, but it can be seen that their E ele ) terms correspond to our Δε elec,A and Δε elec,B , respectively. Jimeńez-Gravalos and Suaŕez further decompose the inter-fragment electrostatic contribution into a dominant term E ele AB (ρ A 0 , ρ B 0 ) that tends to the overall ΔE elec at long distances, and a residual one ) which, together with the intra-fragment contributions, accounts for the charge-penetration energy. We shall see that Δε elec,AB from eq 31 also converges smoothly to ΔE elec at long inter-fragment distances, so for the present purpose we do not consider such additional decomposition.
Subtracting the contributions of eq 32 from those originating from the IQA decomposition of ΔE Pauli + ΔE elec finally yield the appropriate real-space decomposition of the Pauli repulsion term The final EDA picture is completed by the inclusion of the preparation energies from eq 4 and, if required, a dispersion correction. In the case of the semiempirical dipole−dipole model of Grimme, the dispersion correction is added to the interaction energy and has no influence in the intermediate steps, being trivially decomposed by construction as This will not be the case if one uses more sophisticated density-dependent dispersion corrections such as VV10. 50 Finally, the basis set superposition error (BSSE) correction can be estimated a posteriori via the counterpoise formula, 51,52 which is also additive Summarizing, the present approach affords a real-space fully additive decomposition into intra-(A or B) and inter-fragment (AB) contributions of all terms occurring in the EDA scheme To construct all EDA states, the wavefunctions of the dimer and the isolated fragments at the optimized and dimer geometries were evaluated with Gaussian16. The pseudostate (A 0 ∪ B 0 ) electronic structure was constructed using the MOs of the isolated fragments at the dimer geometry. This step was performed with the local program APOST-3D, 59 providing its electronic structure information in a formatted checkpoint (.fchk). Transformation of the formatted into unformatted (.chk) checkpoint file was realized with the unfchk tool from Gaussian. Finally, its associated total energy was extracted using the created chk file as starting guess and forcing to skip the SCF procedure (e.g. SCF = (MAXCYCLE = −1) keyword in Gaussian16). In these calculations the symmetry use was fully disabled to prevent any atomic basis set position difference.
Energy decomposition calculations were also performed with the APOST-3D code using the TFVC atomic definition. 33 For the production results, one-electron numerical integrals were realized using 150 radial (Gauss−Legendre quadrature 60 ) and 974 angular Lebedev−Laikov 61 grid points per atom, while twoelectron numerical integrals have been performed using 150 radial and 590 angular grid mesh. Numerical error minimization of the one-center two-electron contributions was performed using the zero-error strategy. 62

■ RESULTS AND DISCUSSION
We have considered the set of intermolecular complexes from ref 20, that essentially includes hydrogen bonded species, cationdipole, cation−π, halogen−π and π−π interactions between dimers. In addition, we have also considered several anion-π complexes from Quinõnero et al. 63 Except for the π−π interactions, one can identify electron donor and acceptor moieties, which entails certain charge-transfer upon complex formation. We will henceforth refer to fragment A as the net acceptor of charge and fragment B as the donor of charge.
Let us start by analyzing the electrostatic contribution of EDA. Table 1 gathers the IQA decomposition of the electrostatic contributions for the whole set of systems at equilibrium. Note how cation−π and anion−π interactions result in similar values of ΔE elec , but their IQA decomposition reveals a completely different mechanism. In the former, the cation is largely stabilized (large and negative Δε elec,A values) upon complex formation because of its interaction with the frozen density of the donor. This is largely compensated by a positive interfragment electrostatic term Δε elec,AB , that becomes more repulsive as the equilibrium inter-fragment distance shortens from K + to Li + . In the case of anion−π interactions, the Δε elec,B contribution is positive, in line with the situation of Cl − in NaCl previously depicted in Figure 1, and the inter-fragment term is positive. In dispersion bound systems, both the overall electrostatic and their IQA components are very small, and in most cases within numerical accuracy. In the hydrogen bonded and cation-dipole systems, one cannot see a clear trend neither for the intra-fragment contributions (in almost all cases negative) nor for the Δε elec,AB values.
To shed light into the origin of these numerical values, we have also considered the evolution of ΔE elec as well as its IQAdecomposed terms along the dissociation pathway of representative intermolecular complexes. The results are depicted in Figure 2. As it is well-known, when the frozen densities of the two fragments are brought at the complex's optimized geometry, ΔE elec is favorable 9 and the shorter the inter-fragment distance, the more negative the total ΔE elec contribution. The real-space decomposition of ΔE elec yields further insight on this interaction. As previously discussed, the intra-fragment contributions originate from the net electrostatic potential of one fragment interacting with the density of the other fragment that is able to penetrate into its domain. These terms are strongly attractive, particularly in the case of the acceptor A (blue curves in Figure 2), as the more ρ B 0 (r) is able to penetrate into the Ω A domain, the more negative the Δε elec,A contribution becomes. Thus, Δε elec,A is enhanced as the interacting fragments come closer in all cases. Furthermore, this contribution is much larger for cationic acceptor species than for neutral ones (notice the different scales in the examples of Figure 2).
In the case of Δε elec,B (donor of charge), the trend is similar but the magnitude is much smaller, as the amount of density from the acceptor A able to penetrate into the donor is much reduced. In the case of a donor interacting with a hard cation like Li + , this term is essentially zero at all interatomic distances (see Li + ···H 2 O and Li + ···C 6 H 6 curves in Figure 2). However, when the donor is anionic, the trend for Δε elec,B is completely reversed. Since N B 0 > Z B its net electrostatic potential on Ω B can be negative, and thus any ρ A 0 (r) able to penetrate into Ω B leads to positive Δε elec,B values. This effect is clearly seen in the C 6 F 6 ···F − case of Figure 2.
The usefulness of the IQA decomposition of ΔE elec is most clearly seen in the case of Li + ···C 6 H 6 . As shown in Figure 2, the yellow curve is surprisingly flat, and even becomes less attractive at very short distances, totally at odds with the expected behavior. Yet, the overall picture of the intra-and inter-fragment contributions for this system is strikingly similar to that of Li + ··· H 2 O or Na + ···C 6 H 6 . Close inspection to Figure 2 indicates that the behavior of ΔE elec in Li + ···C 6 H 6 can be explained by an insufficient enhancement of the intra-fragment contribution of Li + at short distances.
It is worth to note that both Δε elec,A and Δε elec,B tend asymptotically to zero as the inter-fragment distance increases. This is the expected behavior since at large distances the fragments are essentially in their reference state. Consequently, the inter-fragment Δε elec,AB contribution tends to the overall ΔE elec value. As the distance decreases, however, Δε elec,AB becomes less favorable and even repulsive at very short distances. Thus, the Δε elec,AB value for a given complex at equilibrium geometry may be slightly positive (e.g. Li + ···C 6 H 6 ) or negative (Li + ···H 2 O), but the behavior of the components is analogous in both cases.
Still, the Δε elec,AB contribution at equilibrium distance is very sensitive to the nature of fragments A and B. When both A and B are neutral, the electron-nuclear attraction compensate the nuclear−nuclear repulsion and the Δε elec,AB values are very small (ca. ± 2 kcal/mol). However, when the donor B is anionic the picture is reversed and at equilibrium Δε elec,AB is negative. The case of C 6 F 6 ···F − behaves opposite to the other systems (i.e. it becomes more negative as the inter-fragment distance decreases). The second term on the r.h.s. of eq 31 is key to explain this behavior. Since B is an anion, ρ B 0 (r) holds an excess of electrons with respect to Z B . In addition, the net potential of A in Ω B is governed by the nuclear contribution (hence positive). In that scenario, the closer the fragments, the larger the potential and consequently, even if part of ρ B 0 (r) smears into Ω A , the more dominant the negative term becomes.
Let us proceed by analyzing the Pauli repulsion EDA term, ΔE Pauli , whose contribution originates from the intermediate state A 0 B 0 . Bickelhaupt and Baerends showed that the antisymmetrization of the frozen fragment densities to build A 0 B 0 induces an electron density flow from the intermolecular region to the atomic regions. 9 By decomposing ΔE Pauli into kinetic (ΔT 0 ) and potential (ΔV Pauli ) terms, they showed that the contraction effect translates into an increase of the kinetic energy, and a concomitant decrease (more negative) of the potential energy. The latter is due to the fact that more density is accumulated at regions (e.g. close to nuclei) where the Coulombic potential is larger. The IQA decomposition of ΔE Pauli recovers this picture from a real-space perspective. By definition, kinetic energy contributions only have intra-fragment character upon IQA decomposition and, consequently, they are captured by the Δε Pauli,A and Δε Pauli,B terms. In other words, the Δε Pauli,AB term solely contains potential energy contributions. The kinetic energy increase is so dominant that these terms are expected to be positive and increase along the shortening of the inter-fragment distance. This is exactly the behavior depicted in Figure 3. Furthermore, the corresponding values at equilibrium distance ( Table 2) indicate positive contributions for both the donor and the acceptor fragments with very few exceptions.
On the other hand, the Δε Pauli,AB contributions are large and negative in all cases, and also become more favorable at shorter distances. The origin of this behavior is that, according to eq 33, this term does not explicitly contain energy differences between the intermediate and isolated fragment's states, as there is no inter-fragment term associate to the latter. Deeper analysis indicates that the classical part of the potential energy differences cancels (particularly in the neutral complexes), so the inter-fragment exchange−correlation contribution becomes the dominant term.
Note that the aforementioned contraction effect also increases (becomes more negative) the overall exchange−correlation energy of A 0 B 0 with respect to that of A 0 and B 0 . The dominant exchange contribution is governed by the density close to the nuclei, by virtue of its ρ(r) 4/3 dependence. It might appear counterintuitive that a charge depletion in the inter-atomic region leads, nevertheless, to a negative inter-fragment exchange−correlation. As pointed out by Salvador and Mayer, neither the bond order nor the Hartree−Fock exchange energy components are directly related to overlap populations, but to part of the density localized on the atoms that leads to a correlation between the fluctuations of the atomic populations, even in the absence of overlap. 25 Indeed, inter-atomic exchange energy contribution in the Salvador−Mayer KS-DFT IQA formulation originates on the bond order density between a pair of atoms, which is actually large in the vicinity of the nuclei. 25 An even simpler explanation is that part of the exchange− correlation energy of the A 0 B 0 state is assigned to inter-fragment character by the IQA decomposition, while, once again, there is no inter-fragment contribution from the isolated fragments to compensate for it, as shown in eq 33.
The net result (see Table 2) is that the Δε Pauli,AB contributions are systematically large and negative. On the contrary, the intrafragment Δε Pauli,A and Δε Pauli,B terms are positive, mimicking the behavior of ΔT 0 , but bearing not just kinetic but also intrafragment electrostatic and exchange−correlation contributions. Beyond the overall trend, it is not easy to compare the values of the inter-and intra-fragment contributions from one system to another, especially among different interaction types. Again, even though the behavior of the IQA components is analogous for all interaction types, the actual numerical values are largely dictated by the respective equilibrium distances.
The orbital interaction from EDA, ΔE orb , originates from the relaxation of the MOs of the complex's intermediate state A 0 B 0 to the final complex's AB ground state. It is by definition a negative contribution (if the final state of AB is the ground state), that compensates for the repulsive Pauli term. At short intermolecular distances the intermediate state A 0 B 0 is higher in energy, so that the relaxation energy to the state AB becomes more negative. This behavior can be observed in Figure 4 (yellow curve) for all systems. The orbital relaxation induces an increase of electron density in the inter-atomic (and thus intermolecular) region, making the inter-fragment exchange− correlation contributions of the AB ground state larger (in absolute value) compared to the ones from the intermediate state A 0 B 0 . This is captured by the Δε orb,AB term (grey curve in Figure 4), that closely follows the trend of the global ΔE orb value, with the exception of the C 6 F 6 ···F − system but for reasons that will be disclosed later. The trends observed for the intrafragment terms (blue and orange curves) vary according to the nature of the donor and acceptor moieties. The intra-fragment contribution of the electron donor, Δε orb,B , vanishes at long distances but as the fragments approach it becomes destabilizing. At distances much shorter than the equilibrium the term becomes less repulsive and can even be stabilizing in the case of the water dimer. On the contrary, the intra-fragment contribution for the acceptor, Δε orb,A , is very small (particularly at equilibrium distances) but usually stabilizing along the dissociation profile.
The decomposition of ΔE orb at the equilibrium geometries can be found in Table 3. It is well-known that the ΔE orb contribution accounts for both polarization and charge-transfer effects from the intermediate to the final state. It is precisely the amount of charge-transfer that largely dominates these intrafragment contributions to ΔE orb . The more charge is transferred to the acceptor A going from the intermediate A 0 B 0 state to the final state, the more stabilizing the Δε orb,A contribution, as shown in Figure S1 of the Supporting Information. In the case of the donor moieties the correlation is not as good, but the contributions follow the same trend: the more charge is transferred to the acceptor, the more destabilizing the Δε orb,B values are.
Remarkably, the anion−π systems exhibit an opposite trend. The anion donates charge upon interaction, yet the Δε orb,B contribution is stabilizing. This holds along the whole dissociation profile, as shown in Figure 4. At the same time, the acceptor gains charge but its Δε orb,A contribution is destabilizing. One can also note in Figure 4, the wrong asymptotics of the intra-and inter-fragment contributions for C 6 F 6 ···F − at long distances. This is in fact a clear fingerprint of delocalization error in the KS density, coming from the BP86 functional. First, the dissociation profile could not be further extended at longer distances due to severe SCF convergence problems but, most importantly, the partial charge on F − actually increases from a value of −0.841 at 3.79 Å distance to −0.865 at equilibrium distance, which might explain the aforementioned opposite trend of these systems. It is beyond the scope of the present work to examine the dependence of the decomposed terms on the underlying density functional approximation, but it appears the chosen level of theory is not particularly appropriate to describe these anion-π interactions.
Of course, since the present EDA−IQA decomposition is fully additive, one can obtain the intra-and inter-fragment decomposition of the total interaction energy, ΔE int , by adding the corresponding electrostatic, Pauli repulsion and orbital interaction terms (and dispersion, if included). Numerically, this is not necessary as one can simply perform a conventional IQA decomposition of the final AB state of the complex and subtract the isolated fragment's energies of A 0 and B 0 to obtain the intrafragment or deformation contributions.
For completeness, the IQA decomposition of ΔE int along the dissociation profile of the representative systems is shown in Figure 5, while the corresponding values at equilibrium geometries are gathered on Table 4. Similarly to the orbital interaction contribution, the electronic deformation energies (Δε def.el,A and Δε def.el,B ) at equilibrium are governed by the amount of charge transfer, in this case between the final state and the that of the isolated free fragments. Note that this charge-  Table 4 are larger (in absolute value) as well. The correlation between the electronic deformation energies of both the donor and acceptor moieties and the respective amount of charge-transfer is excellent (r 2 = 0.95, see Figure S2 of the Supporting Information). However, the correlation curve does not cross the (0,0) point but slightly above. That is, even though the acceptor A can eventually gain a small amount of charge (e.g. 0.05e for Na + in Na + ···C 6 H 6 ), the corresponding electronic deformation energy is still slightly positive (+2.7 kcal/mol), due to the accompanying polarization of the fragment's density within the complex. Finally, as usual in the conventional IQA analysis, the Δε int,AB contributions are largely stabilizing along the dissociation profile and also at equilibrium, even for the dispersion−bound complexes (notice that the interaction energies in Table 4 do not contain the dispersion correction). There is also a decent correlation (r 2 = 0.82, see Figure S3 of the Supporting Information) between the Δε int,AB and ΔE int values at equilibrium, even considering the unreliable anion-π complexes. So far, we have presented cases where the EDA-IQA terms are conveniently grouped to match the fragments selected in the EDA step straightly. However, the pairwise nature of the IQA terms can also be used to identify the directionality of each specific interaction. As an example, we have chosen the series of lithium carbanions: LiCF 3 , LiCHF 2 , LiCH 2 F, and LiCH 3 . For  Table 5 summarizes the EDA results for the Li−C bond using the ionic reference fragments Li + ( 1 S) and CR 3 − ( 1 A 1 ). The dissociation energy values show that the more hydrogen atoms are in the molecule, the stronger the Li-CR 3 bond is. The interaction energy follows the same trend, as the preparation energy represents small energy penalties upon deformation. The nature of the chemical bond from the reference fragments is mainly ionic, since the electrostatic interaction represents between 85 and 90% of the total stabilizing interactions. All Pauli repulsion, electrostatic and orbital interaction terms increase in absolute value going from LiCF 3 to LiCH 3 . Applying EDA-IQA to the series brings further insight into the reason behind these trends. The results are gathered in Table 5. The IQA decomposition of ΔE elec shows a large stabilization of the Li + moiety, that increases along the series. As discussed before, this stabilization originates from the CR 3 − density penetration into the Li + domain. The presence of highly electronegative F atoms in the CF 3 − moiety reduces the electron density at the carbon atom, which is closer in space, compared to the less electronegative H atoms in CH 3 − . With the frozen densities of the CR 3 − and Li + fragments interact, the more H atoms the more charge penetration into the Li + domain, resulting in an enhanced Δε elec,Li contribution. The inter-fragment Δε elec,Li−CR3 contribution is also stabilizing, an also increases along the series. Further IQA decomposition into atomic and diatomic contributions helps to rationalize the trend. The contribution of each Li···F interaction to the electrostatic term (ca. 110−115 kcal/mol) is significantly larger than the contribution of the less ionic Li···H counterparts (ca. 20−30 kcal/mol). However, it is the direct Li− C interaction that is most affected by the nature of the substituent R. While this term is largely destabilizing in LiCF 3 , cancelling out to some extent the Li···F interactions, it is even slightly stabilizing in LiCH 3 . The EDA-IQA decomposition provides actual quantification for the qualitative argument that an electron deficient C atom (like in CF 3 − ) will exhibit electrostatic repulsion with the cationic Li + moiety.
The orbital relaxation represents only ca. 10% of the total attractive interactions. Its IQA decomposition reveals that the leading term is the inter-fragment Δε orb,Li-CR3 , but the trend is explained by the contribution of the Li + fragment. The orbital interaction term gathers both charge-transfer and polarization effects. When using charged reference fragments, like in this case, charge transfer should dominate, because the ionicity of the final state is reduced upon bond formation. Going from LiCF 3 to LiCH 3 species the bond ionicity is reduced, hence more charge is transferred to the Li + fragment, enhancing its contribution from −0.6 to −6 kcal/mol, respectively. The inter-fragment contribution to the ΔE orb is still the leading term, but it remains essentially constant along the series, because the orbital relaxation effects are essentially the same for the Li···F and Li···H contacts.
Regarding to the Pauli repulsion EDA term, its total value monotonically increases from Li + ···CF 3 − (35.2 kcal/mol) to Li + ···CH 3 − (51.4 kcal/mol). As previously discussed, the highly repulsive kinetic energy contribution is gathered by the atomic EDA-IQA terms, while the stabilizing potential energy contribution spreads more importantly over the diatomic terms. In this case, the trend along the series is better captured by the highly repulsive contribution of Li + . The less electron rich C atom of CF 3 − leads to a decreased electron reorganization on the Li + fragment upon orthogonalization and antisymetrization to build the intermediate state. On the other hand, it is interesting to note the difference between the contributions of the interatomic Li···C and Li···H/F terms. The former is strongly stabilizing, because the orthogonalization and antisymmetrization affects to a larger extent the inter-fragment region (i.e. the Li−C σ-bond). On the contrary, the same electron reorganization weakens the Li···H and Li···F diatomic terms, and their contribution to the Pauli term is repulsive. The balance of these interatomic contributions makes the interfragment contribution to the Pauli repulsion indeed favorable but almost along the series.
Finally, it is worth to note that the atomic or group contributions of the EDA terms are numerically affected by the particular shape of the atomic weight functions. It is known, albeit not much discussed in the literature, that in the case of Hirshfeld-type approaches the values of the atomic weight functions of light atoms on top of the nucleus is not exactly one, as it is the case for QTAIM or TFVC schemes. For instance, the weight of C atom on top of each H nucleus in CH 4 using conventional Hirshfeld scheme is ca. 0.1. 26 Hence, the charge penetration predicted by Hirshfeld-type schemes can be much larger. Charge penetration is not too significant in terms of electron population/charge, as it is very small compared to the overall atomic density. However, its effect on the energetics of the electrostatic term (i.e. eqns. 30 and 31) can be much more relevant, as the nuclear position is precisely where the nuclear potential is larger.  However, charge penetration effects are even more dramatic in the case of Pauli repulsion, again particularly for the charged species. In this case, it originates from the fact that the individual terms are obtained by differences between IQA energies of the isolated fragments and those of the molecular complex, where the fragments share the physical space. Such interprenetration is much more significant using HI and as a consequence, Δ Paulil,Li for Li + ···H 2 O changes from +20.6 kcal/mol using TFVC to up to 72.1 kcal/mol using HI. On the other hand, the decomposition of the orbital interaction contribution leads to very similar results for both TFVC and HI, precisely because the terms are obtained by comparing IQA energies of the molecular complex at the final and intermediate state (i.e. in both cases the fragments share the physical space). Adding up all EDA terms obtained with HI leads to a final EDA-IQA decomposition of the interaction energy with much larger intra-and inter-fragment contributions of different sign that compensate each other. For this reason, we do not recommend the use of Hirshfeld-type schemes for the present EDA-IQA scheme.

■ CONCLUSIONS
In this work, we have presented the implementation of the IQA decomposition of the individual terms arising from the Kitaura-Morokuma (KM) EDA methodology, namely electrostatic, Pauli repulsion, and orbital interaction. The EDA-IQA approach has been illustrated for a set of complexes, covering different types of intermolecular interactions. In this context, the atomic and diatomic contributions obtained for each EDA term have been conveniently grouped into intra-fragment and inter- Figure 5. Energy evolution (in kcal/mol, y-axis) of ΔE int (yellow) and its IQA-decomposed terms, i.e. Δε def.el,A (blue) and Δε def.el,A (orange) and Δε int,AB (grey), along the dissociation pathway (in Å, x-axis) of H 2 O···H 2 O, Li + ···H 2 O, NH 4 + ···H 2 O, Li + ···C 6 H 6 , Na + ···C 6 H 6 and C 6 F 6 ···F − molecular systems. Equilibrium distance marked with a vertical line. fragment terms. Although the EDA-IQA terms can be grouped to match the selected fragments, the methodology affords atomic and pairwise interactions between all atoms of the system in all EDA terms, helping the precise identification of ruling effects. This in-depth analysis affords a better rationalization of the trends of the bonding along the LiCF 3 to LiCH 3 series. Through the lens of real-space analysis such as IQA, the electrostatic interaction from EDA can no longer be seen as intermolecular in nature, but also results in meaningful and nonnegligible intra-fragment contributions, because the interacting fragments share the physical space once the complex is formed. The EDA-IQA decomposition of the Pauli repulsion shows destabilizing intra-fragment contributions, particularly in the case of fragments that are net charge acceptors. On the contrary, the inter-fragment Pauli contribution is strongly stabilizing. The intra-and inter-fragment ΔE Pauli contributions closely mimic the behavior of the classical decomposition of Pauli repulsion into kinetic and potential terms, respectively. In the case of the orbital interaction term, the sign and magnitude of the intra-fragment contribution at equilibrium geometries is largely driven by the amount of charge transfer: the net acceptors of charge stabilize and the donor moieties destabilize. The proper asymptotics profile of all EDA-IQA terms is also confirmed along the intermolecular dissociation path. Finally, while this work focuses on the particular implementation of (KM) EDA-IQA analysis for intermolecular interactions, it can be readily applied to other EDA schemes relying in intermediate states such as ALMO-EDA.
The Supporting Information contains the correlations between the intra-fragment orbital interaction and interaction energy components (IQA), with the fragment charge differences associated to their associated states ( Figures S1 and S2, respectively), the correlation between the total interaction energy (EDA) and its inter-fragment (IQA) component ( Figure S3), fragment charges (TFVC) for each of the states of the EDA process (Table S1), fragment charge differences (TFVC) associated to each EDA energy term (Table S2) and a comparison of TFVC and Hirshfeld-Iterative EDA-IQA results for selected systems (