A Quantum Chemical Topology Picture of Intermolecular Electrostatic Interactions and Charge Penetration Energy

Based on the Interacting Quantum Atoms approach, we present herein a conceptual and theoretical framework of short-range electrostatic interactions, whose accurate description is still a challenging problem in molecular modeling. For all the noncovalent complexes in the S66 database, the fragment-based and atomic decomposition of the electrostatic binding energies is performed using both the charge density of the dimers and the unrelaxed densities of the monomers. This energy decomposition together with dispersion corrections gives rise to a pairwise approximation to the total binding energy. It also provides energetic descriptors at varying distance that directly address the atomic and molecular electrostatic interactions as described by point-charge or multipole-based potentials. Additionally, we propose a consistent definition of the charge penetration energy within quantum chemical topology, which is mainly characterized in terms of the intramolecular electrostatic energy. Finally, we discuss some practical implications of our results for the design and validation of electrostatic potentials.


INTRODUCTION
Electrostatic interactions are central to molecular modeling because of their slow decay and strength. Especially when polar atoms or charged species are involved, they largely determine the stability and activity of biomolecules such as proteins, nucleic acids, or lipids, among others. 1,2 As such, a reliable description in molecular mechanics (MM) potentials is essential both in the short-and in the long-range.
Within the framework of MM methods, interactions comprising nonbonded atoms are usually represented by pairwise potentials such as the Lennard-Jones and the Coulomb ones. In the latter case, the use of point charges or higher order multipoles to avoid the integration of interacting charge densities has resulted in accurate electrostatics at longrange, with significant improvements to speed up and facilitate convergence such as the Ewald summation and its variants to perform, for example, molecular simulations in solution under periodic boundary conditions. 3−8 At short-range, however, the approximations taken for long distances become less accurate or invalid, 9 and a correct electrostatic description in this regime stills poses a challenge. Hence, there is a growing interest in improving short-range electrostatics (e.g., for troublesome hydrogen bonds), mainly focused on capturing the effects associated with the non-negligible interpenetration of densities, leading to the so-called charge penetration (CP) energy, which is typically defined as the difference between the electrostatic energy computed from continuous charge density distributions and that from multipolar approximations. 10 Thus, several investigations have been devoted in the last years to incorporate the charge penetration energy into the MM electrostatic potentials. 10−14 The separation of various energy terms as implemented in the MM potentials is somehow paralleled by the energy decomposition analysis (EDA) methods. 15 A major goal of any EDA approach is to ascertain the nature and type of the interactions among molecules as well as to rationalize their stabilizing or destabilizing roles, which may have implications for the design, parametrization, and validation of MM potentials such as the electrostatic ones. However, there is no unique recipe to decompose the energy, and thus many EDAs have been developed rooted in different approaches. Hence, symmetry-adapted perturbation theory (SAPT) makes use of a perturbative approach to differentiate the distinct nature of the intermolecular interactions, 16,17 while orbitalbased EDAs exploit a stepped scheme to calculate the different energies according to some reference electronic states, 18−20 and the interacting quantum atoms (IQA) method relies on a real space partition of the quantum mechanical (QM) density matrices, 21,22 being thus classified as a quantum chemical topology (QCT) method.
According to recent studies, in spite of their crude approximations, it may be feasible to improve the classical MM potentials by utilizing the information provided by EDAs. 23 More specifically, it has been shown that the SAPT energy components (electrostatics, induction, exchangerepulsion, and dispersion) can be modeled with relatively simple MM functions. 24,25 In particular, it has been demonstrated that the combination of empirical damping functions accounting for the CP energy with point multipoles results in electrostatic energies at short-range that are quite close to the SAPT ones. Actually, the SAPT electrostatic energy provides the required reference to parametrize and validate the CP-augmented potentials. However, different interpretations of short-range energetic effects involving the overlap of the electron densities of two or more fragments may be possible depending on the particular EDA of choice. 15 As such, other schemes such as the absolutely localized molecular orbital (ALMO) EDA, that relies on a different nonperturbative decomposition of energy terms, have also been proposed. 26 In this work, we reexamine the nature of electrostatic interactions under the prism of an orbitalinvariant, reference-free technique. The IQA approach fulfills these requirements as it is a QCT, real-space energy decomposition resorting to the partition of the reduced density matrices (RDMs). IQA distinguishes not only between electrostatic or exchange-correlation components of the interaction energy but also between intra-or interatomic (or fragment) contributions. Moreover, since IQA splits the total energy of a system and not only the interaction between selected fragments, it is capable of reconstructing (or dissecting) the energy ascribed to both covalent and noncovalent binding, allowing thus covalent bond energies to be characterized 27 as well as the accuracy of the energy components handled by QM fragment methods to be investigated. 28 Herein, we study in detail the electrostatic interactions involved in noncovalent complexes with a twofold goal. On the one hand, we aim to compare in a consistent and systematic manner the atomic and fragment contributions to the electrostatic energy as evaluated throughout a hierarchy of QM and MM approximations and at varying intermolecular distances. In this way, we seek to identify the best correspondence between the IQA and the MM electrostatic terms. On the other hand, we critically examine the CP concept and propose a novel definition relying on a joint orbital and real-space decomposition scheme, which can give new insight into the CP energy. To help fulfill these goals, the rest of the manuscript is structured as follows. First, we present and describe the theoretical scaffold that holds our work, paying particular attention to the IQAand its IQF variant energy decomposition, followed by subsections concerning the zeroth-order approximation, the electrostatic MM potentials, and a final assessment of the CP energy and the alternative definition proposed in this work. Subsequently, we describe some computational settings and the results of our test calculations, which were carried out on the S66 and S66x8 data sets. 29, 30 The various levels of description of the electrostatic interactions are then discussed based on their statistical correlation with benchmark data, their dependence with the intermolecular separation, etc. The QM and IQA calculations yield further information, not only about the magnitude of the CP energy, but more importantly, about its different role in the IQA descriptors. Finally, we conclude that the aim of improving the electrostatic description is essentially fulfilled at the expense of accounting for intramolecular effects.

THEORY AND METHODS
2.1. IQA Decomposition of QM Energies. The interacting quantum atoms method is a robust and physically sound approach to decompose the total QM energy of a system into chemically meaningful components. 21,22 It is based on partitioning the firstand second-order RDMs, as can be done with the real space partition proposed by Bader and coworkers within their Quantum Theory of Atoms in Molecules (QTAIM). 31 Thus, the three-dimensional space is decomposed into atomic regions (Ω I ) as the attraction basins of the gradient field of the electron density.
Given a global energy E of a system, IQA permits its decomposition into atomic components and pair interaction energies according to where E net I is called the net atomic energy and, under the Born−Oppenheimer approximation, represents the energy due to the kinetic energy of electrons plus all the interactions involved (i.e., electron−electron and electron−nucleus) inside the atomic basin of each atom I. Similarly, each E int IJ term comprises the interaction energy between the electrons (e) and nucleus (n) located in atom I with those ascribed to other atoms J, which can be separated into n−e, e−e, and n−n contributions.
In order to compute the potential energy, the pair density ρ 2 (r 1 , r 2 ) is required. This object can be split according to ρ 2 (r 1 , r 2 ) = ρ(r 1 )ρ(r 2 ) + ρ xc (r 1 , r 2 ) in two contributions. On the one hand, ρ(r 1 )ρ(r 2 ) represents a noncorrelated product of densities, whereas electron correlation is accounted for by the exchange-correlation (xc) density ρ xc (r 1 , r 2 ). Accordingly, the total interaction energy between two atoms can be decomposed into a Coulomb or electrostatic term E ele IJ and a quantum mechanical exchange-correlation one E xc IJ : the latter term comprising all the associated QM effects that other (e.g., perturbative) approaches identify separately as dispersion, charge-transfer, polarization, etc. However, such a decomposition of E int IJ into two terms is particularly relevant when assessing the nature of a given bond or interaction, since the electrostatic term is associated with ionicity and the exchange-correlation contribution with covalency. 22 IQA admits the grouping of atomic terms into fragment contributions (e.g., functional groups, molecules). Thus, a fragment decomposition similar to eq 1 of a molecular aggregate constituted by two moieties A and B involves where E net B can be calculated analogously. For practical purposes, we use the IQA acronym to refer to the atomic analysis, whereas for its fragment version the term interacting quantum fragments (IQF) is preferred.
In a previous work, 32 it was shown that IQF may be useful to dissect the formation energy of noncovalent complexes. Moreover, the IQA/IQF terms can be augmented with Grimme's D3 dispersion correction 33  ) energies between the two fragments, the latter being thus separated from the whole exchange-correlation one. Overall, the contribution of electrostatics and exchange-correlation to ΔE form is split between the intrafragment deformation and the interfragment interactions.
2.2. Electrostatic Energy from Continuous Charge Densities. The purely electrostatic energy for a given system with total charge density ρ(r) (ρ(r) ≡ ρ tot (r) = Z r R r ( ) ( ) , including both the electron density ρ e (r) and the nuclear charges Z I at positions R I ) is readily computed using the Coulomb law,   6) where, for the sake of simplicity, the electrostatic potential in this and the rest of the equations is expressed in atomic units. Interestingly, the QTAIM real space partition derived from ρ e (r) allows us to decompose the electrostatic energy at the atomic level, Similarly, the fragment-based decomposition can be readily accomplished in an analogous way, allowing thus the specific assessment of the electrostatic component of the formation energy ΔE ele of a two-fragment system AB as where ΔE ele is expressed in terms of two contributions, namely, the intrafragment variations of electrostatic energy in the formation process, E ele A Δ and E ele B Δ , and the interfragment electrostatic interaction, E ele AB . At this point, we note that although ΔE ele is commonly termed as a classical electrostatic interaction energy, we will refer to it as the electrostatic contribution to the formation energy of the A···B complex in order to help avoid confusions with the IQA/IQF interaction energy terms. When the charge density is constructed from the unrelaxed fragment densities as r r r ( ) ( ) ( ) , the electrostatic contribution to the formation energy, which is named here as the zeroth-order energy E ele 0 Δ , equals the Coulomb interaction between the unrelaxed densities: This energetic term corresponds to the so-called f irst-order polarization energy (or simply electrostatic energy) defined in SAPT, 16 which has been adopted as a benchmark electrostatic energy for the validation of recently developed short-range electrostatic potentials.
2.3. Electrostatic Potentials in Molecular Mechanics. To avoid the usage of continuous charge distributions, the MM methods typically invoke the multipolar expansion as detailed in the Supporting Information (SI), which approximates the zeroth-order energy defined in eq 9. Formally, the multipolar electrostatic energy E ele mp , 0 Δ is affected by two different error sources. On the one hand, the underlying expansion must be truncated at some order (l max = 0, 1, 2, ...), resulting thus in a certain truncation error. On the other, when r ( ) overlap to a significant extent, the rigorous application of the multipole expansion is impeded so that its usage at short distances implies some charge penetration error, which is normally assumed to be dominant. Nevertheless, the multipole-based potentials are still largely useful in many cases, and they enhance convergence by distributing multipoles throughout the molecule at the atomic sites and/or bond centers. 9,35,36 The MM electrostatic potentials can be classified into two groups. On one side, MM methods such as AMBER, 37 CHARMM, 38 GROMOS, 39 and OPLS 40 adopt simple electrostatic formulas with point charges (i.e., monopoles, with l max = 0) that are ultimately the result of a fitting procedure against the molecular electrostatic potential (ESP). On the other side, more sophisticated methods, such as NEMO, 41 AMOEBA, 42 or the QTAIM-based FFLUX, 43,44 include higher order multipoles (frequently up to the quadrupoles, l max = 2) in order to capture the anisotropy of the distribution of electrons in space. These potentials are generally built from the QM density matrix of the molecule of interest by means of the distributed multipole analysis (DMA) 36 or similar procedures. In addition, some methods (e.g., AMOEBA or NEMO) also refine the DMA multipoles to better reproduce the ESP values. In this way, the resulting charges/multipoles may include in an effective way both high-order multipolar contributions and CP effects. Actually, the performance of the MM potentials is examined statistically as a whole (i.e., using the full MM potential including all bonded and nonbonded terms) by various energetic and structural validation tests. A quite different approach is followed by the FFLUX force field. It makes use of QTAIM multipoles in contrast to the more widespread DMA methodologies and estimates them by means of a machine learning technique depending on each atom's environment.
In comparison with the atomic/multipolar methods that are massively employed in current simulation packages, the electrostatic MM potentials that go beyond the multipolar approximation are much less consolidated. In this category, we find different methodologies such as SIBFA, 45 EFP, 46 and AMOEBA+ 24 that complement the multipolar formulas with other (so-called damping) functions to capture very-shortrange electrostatics and to remove the CP error. In this way, these potentials (whose general form is shown in the SI) seek to reproduce E ele 0 Δ as evaluated by SAPT or similar methodologies. 24 2.4. Charge Penetration Energy. The CP energy E pen between two molecules has been defined 47 as the difference between the exact zeroth-order electrostatic energy E ele 0 Δ and its multipolar analogue E ele mp Conceptually, this straightforward definition of E pen is satisfactory. It also shows that E pen is not only an interfragment quantity but rather an energy that presents also intramolecular contributions according to the real space partitioning of the whole E ele 0 Δ . In this respect, the energetic definition suggests that the CP energy is not limited to the change in the electrostatic interaction between two atoms due to their electron cloud overlap and the associated loss of nuclear screening. 48 The rigorous evaluation of E pen for different systems at varying intermolecular separations would allow a deeper analysis of electrostatics and, eventually, the development of more accurate potentials. However, as noticed by Bojarowski et al., 47 dif ferent methods of obtaining multipole moments lead to different radii of (pseudo)convergence, dif ferent levels of multipole expansions at which (pseudo)convergence is achieved, and dif ferent values of penetration energy. Therefore, the value of the CP energy as evaluated with eq 10 may depend on the particular method used to derive the multipoles. Moreover, the usage of truncated expansions introduces some additional truncation error so that both truncation and penetration effects become somewhat mixed in the resulting E pen values. 49 An alternative to evaluate E pen has been proposed by Kairys and Jensen. 50 Having noticed the relationship between the CP energy and the magnitude of the orbital overlap, they attempt to recover such an effect from scratch, with a derivation of E pen independently from the multipolar model used to estimate electrostatics at first stage. However, the authors find that the inherent dependence on the set of molecular orbitals used may lead to different CP values.
2.4.1. A Novel IQF Definition of the Charge Penetration Energy. By combining both the Bader partitioning scheme ( , the following energy terms are obtained: In the notation used above the two interacting densities are encompassed by parentheses, while the basins they are integrated in are identified by the corresponding superscripts in the given order (only one if both are the same). Hence, for i n s t a n c e , t h e t e r m E ( , ) When the above double decomposition is applied to the electrostatic energies of the separate fragments, such as A, in the final complex, the electrostatic energy of the original species becomes Note that this partitioning is derived from the AB zerothorder (i.e., Hartree product) wave function and that geometry relaxation effects are not considered. By subtracting from eq 11 the previous fragment energies, the corresponding electrostatic contribution to the formation energy of the complex is obtained, Among the surviving terms in eq 13, E ( , ) at long distances, while the other three terms would present a similar behavior of increasing in magnitude when shortening the intermolecular distances R AB and canceling out in the opposite situation. Thus, those three terms can be directly related with the interpenetration of molecular densities and grouped in the IQF-like electrostatic charge penetration energy This term fulfills (and so its three Figure 1 represents the previous four terms between the partitioned A term between the total densities in each basin.

COMPUTATIONAL DETAILS
3.1. Molecular Geometries and Reference Interaction Energies. All the QM and classical electrostatic calculations were performed on the molecular geometries retrieved from the S66 database, 29 which contains a set of 66 complexes featuring the most common noncovalent interactions in biomolecules. These can be classified depending on the atoms involved into polar, nonpolar, and mixed. Analogously, the different complexes have been grouped into H-bond, Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article dispersion, and mixed according to the main interactions they experience (see Table S1). For representing both the atomic interactions and the subsets of complexes, a color code has been utilized: magenta for H-bond/polar, yellow for mixed, and blue for dispersion/nonpolar. In addition to the S66 set, a selection of 12 representative complexes from the S66x8 database, 30 which is an extension of the former to eight different fractions of the equilibrium intermolecular distances, were also considered. The benchmark CCSD(T)/CBS interaction energies collected in S66 were employed as the reference values for comparative purposes.
3.2. HF-D3 Calculations. HF/cc-pVTZ calculations were carried out on the S66 and the S66x8 geometries using the GAMESS-US package. 51 Grimme's D3 dispersion potential as implemented in the DFT-D3 code 52 was employed to incorporate the dispersion energy. Additionally, in order to correctly reproduce the asymptotic behavior of the dispersion energy at small distances, the Becke−Johnson damping function was chosen. 53 We selected HF because it lacks entirely dispersion energy and thereby yields a straight physical partitioning of energy in combination with the D3 potential. We also note in passing that HF-D3 has been shown to describe correctly and efficiently the structure and energetics of biomolecules 54 and that a variant of DFT-SAPT has been also developed in which the costly ab initio dispersion calculations are replaced by a reparametrized D3 potential. 55 In addition, the HF-D3/cc-pVTZ energies reproduce quite well the reference CCSD(T)/ CBS energies of the S66 structures (see Figure S1).
3.3. IQA Energy Decomposition Analysis. The decomposition of the QM and the electrostatic energies derived from continuous charge densities were performed with the PROMOLDEN code. 56 The integration settings comprised β-spheres with radii of 60% of the distance between each nucleus and its closest critical point. Within them, Lebedev angular grids with 974 points were used, along with Euler− McLaurin radial quadratures with 382 radial points. A bipolar expansion of r 12 1 − was selected with an l max of 6. On the other hand, the outer part of the basins (i.e., outside the β-spheres) employed the same angular and radial quadratures, albeit increasing their respective points up to 5810 and 512, with a maximum radius of 15 au. In this case, r 12 1 − was expanded by means of a Laplace expansion with l max = 10.
3.4. Point-Charge and Multipolar Calculations. Atomic charges were computed for the separate monomers in the S66 structures by means of the restrained electrostatic potential (RESP) method following the General Amber Force Field (GAFF) 57 prescriptions with a HF/6-31G* level of theory. In the case of the atomic multipoles, two different sets were employed. On the one hand, AMOEBA multipoles were derived up to the quadrupoles (l max = 2) following its corresponding parametrization protocol. 42,58 On the other, QTAIM multipoles were obtained by means of the PROMOLDEN program with an l max = 2. Both the AMOEBA and the QTAIM multipolar energies were obtained with the MPOLINT code. 59  Δ + Δ and the total kinetic plus exchange-correlation contributions ΔT + ΔE xc . The statistical analysis comprises the coefficient of determination R 2 , Spearman's rank correlation coefficient ρ, and the root-mean-square error RMS. Data corresponding to the whole set of complexes is depicted in black and that ascribed to the H-bond group is in magenta, while mixed and dispersion complexes are in yellow and blue, respectively. All the energies are in kcal mol −1 .

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article Additionally, a set of 12 S66x8 complexes was tested under the AMOEBA+ CP-corrected potentials. 24 For this, TINKER was used to calculate the respective CP energies as the difference between the CP-corrected multipoles and the multipolar energies previously derived. The parameters of the damping functions were directly taken from the literature. 24   Additionally, the zeroth-order electrostatic contribution to formation E ele 0 Δ is also included. The complexes are colored and displayed in columns according to the group they belong to, namely, H-bond, mixed, or dispersion, respectively. The energies (Y-axis) are given in kcal mol −1 , and the abscissas represent the intermolecular distances relative to the equilibrium ones (R AB /R eq ).

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article 4. RESULTS AND DISCUSSION 4.1. IQF-D3 Partitioning and Pairwise Approximation. The IQF-D3 decomposition of the HF/cc-pVTZ binding energies for the S66 complexes has been discussed at length in previous work. 32 Herein, we focus on the decomposition of the electrostatic descriptors into intra-and interfragment components. Interestingly, we found that the combination of the interfragment electrostatic interaction energy E ele AB with the D3 dispersion potential yields pairwise energies that are quite well correlated with the S66 benchmark values, the coefficient of determination being R 2 = 0.990 with RMS errors of 5.7 kcal mol −1 (see Figure 2 and Table S2). Thus, the IQF E ele AB descriptors in conjunction with the D3 potential capture the essential electrostatic and dispersion interactions that determine the relative stability of the noncovalent complexes. When addressing both terms independently ( Figure S2), we find that the pure electrostatic E ele AB term exhibits a satisfactory overall correlation (R 2 = 0.943) due to the fundamental role of electrostatics in H-bond complexes. On the other hand, the D3 descriptor has a null global correlation with the S66 reference energies, although it is reasonable (R 2 = 0.820) for the dispersion complexes as expected. However, the mixed complexes are not well-described by either the electrostatic or the dispersion energies separately, and their combination becomes critical.
In contrast to the ability of the E D3 ele AB + descriptors to capture the main features of noncovalent binding, the combination of ΔE ele , which includes both the intra-and the intermolecular electrostatic effects, with the D3 potential deteriorates the global correlation (R 2 = 0.888) and results in larger RMS errors (17.3 kcal mol −1 ). The full IQF decomposition (eq 8) explains this unbalanced description because the intrafragment electrostatic energies, which contribute to the deformation energies, tend to cancel out with the QM energy terms (electronic kinetic energy and exchange-correlation) that are not required in the simple electrostatic + dispersion picture (see Figure 2 (Table S3). This behavior is also satisfactory within the S66 subsets: R 2 = 0.989 and 0.988 for the polar H-bonded systems and the dispersion-dominated complexes, respectively, albeit the correlation is somewhat reduced in the case of the mixed complexes (R 2 = 0.755). Further support for the use of the zeroth-order energies comes from the atomic level, where a high degree of coincidence between the diatomic zeroth-order E ele IJ 0, and fully relaxed E ele IJ energies is also found at the equilibrium geometries (R 2 = 0.995, see the SI).
When addressing the distance dependence of the previous term (see Figure 3), both E ele AB and E ele AB 0, follow the same trends at varying intermolecular separations R AB (given as relative to the equilibrium distances R eq ). As expected, they start diverging at short distances due to the strengthening of charge polarization, charge-penetration, and charge-transfer effects that attenuate the pairwise electrostatic forces. The magnitude of these effects is clearly system-dependent, as well as the shape and slope of the E ele AB and E ele AB 0, curves, revealing thus further details about the role of electrostatics in these complexes. Thus, the electrostatic stabilization of the four H-bond complexes and others (e.g., the π-complex of the uracil dimer) is continuously reinforced upon shortening the monomer−monomer distance, reflecting the major electrostatic control of these systems. In contrast, the T-shaped benzene complexes with methanol or N-methylacetamide reach an electrostatic minimum at a distance longer than the equilibrium one while the small electrostatic energies of the dispersion dimers (i.e., +1, −1 kcal/mol) change very little along the curves (some small leaps are due to residual errors arising in the numerical integration over the atomic basins).
In Figure 3 the deviation between the global E ele 0 Δ energies and the interfragment E ele AB 0, anticipates the underlying CP effects associated with the density overlap. For the H-bond and some of the mixed complexes, the two curves decrease with lowering separation, but they split gradually for R AB /R eq < 1.6. The global E ele 0 Δ stabilization nearly doubles E ele AB 0, at R eq , showing thus the large impact of intramolecular electrostatics as defined in the IQF framework. For the π-complexes (benzene−dimer, benzene−methanol, ...) or the weakly interacting neopentane dimer, the inter-vs intramolecular balance is differently modulated because the deviation between the global and the interfragment electrostatics becomes significant only at very close distances (e.g., R AB /R eq < 1.1), which are indicative of mutual overlap. In these systems, E ele 0 Δ is thus reinforced by several kcal mol −1 , which are ascribed to the intrafragment electrostatic stabilization achieved by the fragment-overlap (i.e., CP) effects. Such effects have a minor influence on the small E ele AB 0, energies (<1−2 kcal mol −1 ), which tend to remain nearly constant or become slightly attenuated. As shown below (Section 4.5), the IQF analysis of the CP energy gives further insight about the behavior of E ele AB 0, and E ele 0 Δ with R AB /R eq .

Comparison between E ele
and Pairwise MM Energies. The pairwise approximation that emerges from the IQF-D3 decomposition and the validity of the zeroth-order approximation for the electrostatic interactions provide an insightful theoretical support for the construction of noncovalent MM potentials. In this scenario, E ele AB 0, can be seen as the most suitable IQF descriptor to assess the approximate electrostatic potentials. Hence, we calculated the interfragment electrostatic energies using the RESP atomic charges and the AMOEBA multipoles, as well as the QTAIM multipoles up to the quadrupoles.
According to the statistical data in Table 1, either the RESP atomic charges or the QTAIM/AMOEBA multipoles give interfragment electrostatic energies that correlate considerably well with E ele AB 0, (R 2 > 0.9 and RMS errors ∼ 1 kcal mol −1 ) for the full S66 set and also for the H-bond/dispersion subsets. These point-charge/multipolar electrostatic energies are less satisfactory for the less abundant mixed complexes, although the multipolar potentials yield a more accurate description (R 2 ≃ 0.6−0.8) than the RESP charges (R 2 ≃ 0.5). In addition to E ele 0 Δ and the fully relaxed and zeroth-order IQF pairwise Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article terms, Figure 3 also displays the distance dependence of the QTAIM/AMOEBA/RESP energies, that results quite close to that of the interfragment E ele AB 0, energies. Nevertheless, a closer inspection reveals that the QTAIM/AMOEBA/RESP energies tend to overestimate the stabilizing/destabilizing character of E ele AB 0, for the H-bond/dispersion dimers, respectively. The good agreement between the multipolar and the RESP energies in Table 1 and in Figure 3 suggests that the RESP fitting procedure may incorporate in an effective way higher order effects even at short distances. In addition, our results point out that the pure QTAIM multipoles can be employed in the construction of accurate electrostatic potentials, free from the inclusion of other effects that may be present when the DMA multipoles are fitted against the molecular ESP. In fact, the QTAIM multipoles, which are already considered in the FFLUX force field, readily reproduce the ESP without the need of any constraint. 63 4.4. Comparing Diatomic Electrostatic Interactions. IQA permits an unambiguous decomposition of the continuous-density intermolecular interaction energy into a sum of atomic and diatomic terms that enables a thorough analysis of the global molecular properties based on their atomic origins, and a close comparison with the various MM descriptions at this atomic level.
As expected, the IQA diatomic terms correlate almost perfectly with the QTAIM multipolar ones E ele mp IJ , 0, (see Figure  4). On the contrary, the AMOEBA and RESP energies are significantly less correlated (R 2 of 0.7 and 0.4, respectively) and have large RMS errors. For example, the largest discrepancies between E ele IJ 0, and the QTAIM-multipolar E ele mp IJ , 0, in the acetic acid dimer (about 6 kcal mol −1 ) arise from the atoms involved in the OH·O H-bonds, the rest of pair interactions having much lower differences (<0.5 kcal mol −1 ; see Tables S7−S9) ), the largest discrepancies amount to hundreds of kcal mol −1 and involve not only short polar contacts but methyl C atoms too (see Tables S10−S15). values was not entirely unexpected given that the RESP charges are derived from the molecular ESP and the AMOEBA multipoles are obtained by the DMA protocol. In fact, a difference of 1 order of magnitude between the atom−atom electrostatic interactions from IQA and MM potentials has also been noticed previously. 64 The present results show in further detail the actual discrepancies between the various atomic representations and suggest that, although the diverse atomic multipoles employed in classical potentials yield similar molecular electrostatic energies, the atomic decomposition is more questionable, which, in turn, can negatively affect the interpretation of localized electrostatic interactions and/or result in artifacts while dealing with QM and MM short-range electrostatics in hybrid QM/MM methodologies. ρ ρ between the tails of each molecular density that penetrate into the opposite basin, as described in eq 14. This constitutes an ef fective penetration energy in the sense that the molecular identity between two overlapping fragments becomes necessarily blurred so that fragment properties are dependent upon the scheme followed to dissect the global charge density into its constituents. Nevertheless, the topological analysis of ρ 0 yields a consistent identification of molecular fragments so that we believe that the associated charge-penetration analysis can give useful insight into the electrostatics of noncovalent complexes.
The application of eq 14 to E ele 0 Δ results in the energy contributions shown in Figure 5. On the one hand, the interfragment energy E ( , )  ∈ Ω ). In light of these results, e−n attraction greatly overcomes e−e repulsion between different zeroth-order densities inside the same basin and gives rise to the significant stabilizing energies observed. There is also a minor repulsive contribution owing to the purely electronic repulsion between the penetrating A Further insight can be gained by analyzing the distance dependence of the various energy terms as shown in Figure 6. The plots confirm that the three components of E pen IQF tend to zero when R AB /R eq > 1.5 and further highlight the role of the intrafragment terms. Interestingly, the E ( , ) ρ ρ is damped out or inverted at the shortest distances. This is not entirely unexpected given that, as two   to the reference E ele 0 Δ , which is equivalent to the SAPT electrostatic energy. In effect, Figure 7 shows that E ele AMOEBA Δ + nearly matches E ele 0 Δ . Concerning the CP energies, it is important to note again that the AMOEBA+ reference for measuring the CP energy is different from that provided by the IQF-QTAIM approach. Nevertheless, the two penetration energies exhibit a similar behavior with R AB , particularly for the more stable H-bond complexes, which resemble also the variations experienced by the intramolecular CP terms, ρ ρ . Therefore, we conclude that the AMOEBA+ CP and similar corrections account mainly for intramolecular electrostatics.

CONCLUDING REMARKS
In this work we have analyzed the short-range electrostatic interactions in the S66 and S66x8 data sets through a hierarchy of approximations at both the molecular and the atomic levels.
We have shown first that the IQA/IQF decomposition augmented with the D3 dispersion terms gives support to  values, which are obtained from the unrelaxed densities of the isolated monomers (i.e., the zeroth-order approximation).
According to our results, the intermolecular E ele AB 0, energy turns out to be the most appropriate IQF descriptor to analyze and/or compare with electrostatic MM potentials. In particular, we have considered two widely used potentials relying on the RESP atomic charges or the AMOEBA distributed atomic multipoles, respectively, as well as the multipolar potential up to the quadrupoles derived directly from the QTAIM basins. The three MM pairwise approximations correlate satisfactorily with the zeroth-order IQF term Figure 7. Comparison between the AMOEBA+ model and the zeroth-order IQF energies for our model S66x8 complexes. The complexes have been displayed according to the group they belong to (either H-bond, mixed, or dispersion). Distances (X-axis) are relative to the equilibrium ones (R AB /R eq ) and energies (Y-axis) are in kcal mol −1 . values are further decomposed into diatomic contributions, large discrepancies between the RESP or the AMOEBA atom−atom interactions and their zeroth-order IQA counterparts are unveiled. Although this is understandable in terms of the specific details of the RESP/AMOEBA charge/multipole derivations, it contrasts sharply with the nearly perfect match between the QTAIM atomic multipolar energies and the IQA reference values. Hence, MM potentials based on the QTAIM multipolessuch as the QCT-based FFLUXmay provide a more consistent description of electrostatic interactions at both the molecular and the atomic levels.

Journal of Chemical Theory and Computation
Besides forging links between the IQF/IQA quantities and the MM electrostatic potentials, we have studied the charge penetration effects that arise from the mutual interpenetration of the zeroth-order molecular densities in their opposite QTAIM basins as built from the final ρ 0 of the complex. This QTAIM perspective allows us to dissect the CP energy into different contributions that emphasize its intramolecular character, which, in turn, is dominated by the attraction between the nuclei of fragment A (B) and the penetrating tail of density B (A). In this way we may clarify some practical issues related with the CP corrections for MM potentials. For example, adding CP corrections to MM potentials like RESP/ AMOEBA, which target the zeroth-order interfragment electrostatic energy, results, necessarily, in an unbalanced description. This aspect, which has been overlooked in previous works, 10,48,65 implies also that the electrostatic energy employed in popular MM force fields (AMBER, CHARMM, ...) cannot be compared with the global E ele 0 Δ energy derived from continuous charge distributions, but with its interfragment component. On the other hand, CP corrections have been derived to improve the description of the QM−MM electrostatic interactions in hybrid QM/MM methodologies. 12 In this case, such corrections should mitigate short-range electrostatic artifacts, particularly those associated with the QM−MM covalent linkages. However, considering the highly dissimilar interatomic electrostatic energies produced by the QM densities and the RESP/AMOEBA potentials, the usage of electrostatic parameters more akin to the QM densities at the atomic level may have a larger impact in improving the QM− MM electrostatics.
Finally, concerning the novel MM potentials inspired by the QM SAPT methodology, it is clear that the multipolar electrostatics (interfragment) must be augmented by the CP potentials (intrafragment) if one seeks to reproduce the global electrostatics E ele 0 Δ . Nevertheless, the IQF/IQA approach (and other EDAs) points out that the intramolecular electrostatic energy is closely related with other energy changes induced by fragment overlap (e.g., deformation and interfragment exchange-correlation energy), suggesting thus that the separate treatment of these effects by means of independent potential terms might be inefficient and hamper parameter development and transferability.