Extension of the interacting quantum atoms ( IQA ) approach to B 3 LYP level density functional theory ( DFT ) †

An interaction between two atoms, bonded or non-bonded, consists of interatomic contributions: electrostatic energy, exchange energy and electronic correlation energy. Together with the intra-atomic energy of an atom, these contributions are the basic components of the Interacting Quantum Atom (IQA) energy decomposition scheme. Here, we investigate IQA’s proper use in conjunction with an explicit implementation of the B3LYP functional. The recovery of the total molecular energy from the IQA components is emphasised, for the first time. A systematic study of three model systems of biological relevance, N-methylacetamide (NMA), the doubly capped tripeptide GlyGlyGly and an alloxan dimer, shows the stabilization effect of B3LYP on most of the interatomic exchange energies (V X ) compared to their Hartree–Fock values. Diagrams of exchange energies versus interatomic distance show the clustering of interactions, one cluster for each 1,n (n = 1 to 6 where the atoms are separated by n 1 bonds). The positioning of some V X values outside their expected cluster marks interesting interactions.


Introduction
The spatial localization and partitioning of properties of matter is central to Chemistry.Fortunately, matter (at ambient conditions, and excluding metallic systems) allows itself to be described in terms of local features and behaviour, as well as recurring motifs.Indeed, molecular assemblies consist of molecules, which in turn consist of functional groups and largely transferable atoms (e.g.ref. 1-3).This locality and transferability enables Chemistry to classify, systematize and eventually explain and predict the phenomena it observes.
There are various alternatives to obtain this local information theoretically and computationally, most popularly by energy and electronic charge.Amongst the most used partitioning schemes are the ''energy decomposition analysis'' (EDA), 4,5 its variant the natural energy decomposition analysis (NEDA), 6 the natural bond orbital (NBO) analysis 7 or Quantum Chemical Topology (QCT), 8 which amongst its various segments incorporates the Quantum Theory of Atoms in Molecules (QTAIM) 9 and Interacting Quantum Atoms (IQA). 10The latter approach is conceptually the most minimal of all energy partitioning schemes.IQA is parameter-free although computationally the most expensive method.However, very recently, IQA has gained popularity through its implementation in the efficient computer program AIMAll. 11The IQA method is the subject of this paper, particularly, its use within the context of density functional theory (DFT).We also note again that IQA can operate on nonequilibrium geometries, and can thus go beyond the original QTAIM atomic energies, which are only valid at a stationary point (e.g. an equilibrium geometry).
QCT is a branch of theoretical chemistry that yields a wealth of calculated chemical information from the wavefunction of a molecule, molecular assembly, metallic system or ionic crystal.It uses the (mathematical) language of dynamical systems (e.g.attractor, critical point, separatrix, basin, gradient path) to obtain chemical insight and then ideally make predictions.It defines the so-called topological atom, which is a finite-volume object with a shape, which is specifically determined by the total system that this atom is part of topological atoms do not overlap nor leave gaps: they collectively exhaust space and their properties are additive.In other words, there is no need for corrections due to topological atoms intersecting, or leaving a gap between them (resulting in electron density not belonging to any one atom).
Each topological atom contributes its own intra-atomic energy to the total system's energy. 10In addition to this contribution, and more importantly for this study, IQA also provides all inter-atomic energy contributions.When summed, energy produced by an IQA partitioning of a B3LYP wavefunction, compared with that of a Hartree-Fock (HF) wavefunction, for the same molecule.
Note that IQA is not the only method that has identified the use of KS orbitals as a concern before.The same concern surfaces in the calculation and interpretation of delocalization indices (DIs), also denoted d(A,B).A delocalization index represents the number of electrons delocalized (i.e.exchanged or shared) between atoms A and B. Hence, it is not surprising that there is a relationship between d(A,B) and the interatomic exchange energy between A and B. Indeed, using a binomial Taylor expansion one can prove 34,35 that the interatomic exchange energy is approximately equal to minus d(A,B) divided by twice the internuclear distance.Thus, a large value for d(A,B) corresponds to a large absolute value for the interatomic exchange energy, and hence energetic stabilization.In 2002, Poater et al. 36 reported the significance of having to use DFT-KS orbitals within a HF-formalism in order to obtain a numerical value at DFT level.They produced d(A,B) values obtained from a HF treatment of a DFT wavefunction in order to gain useful chemical insight.However, as electron correlation was not fully considered, their values were consistently overestimated.Both d(A,B) and the interatomic exchange-correlation energies can provide valuable information of chemical significance for polyatomic molecules. 25,37,38e will present a direct comparison of the exchange energies calculated from a HF wavefunction and calculated from a HF treatment of a B3LYP wavefunction, both at exactly the same molecular geometry.In order to examine a range of both intramolecular and intermolecular interactions between atoms, the following biologically relevant molecules are used as case studies: (i) N-methylacetamide (NMA), which serves as a wellstudied prototype system for the ubiquitous peptide bond, (ii) triglycine (GlyGlyGly, with a peptide-bond cap at both termini) (TriGly), which is of relevance to the construction of our peptide force field QCTFF given the recent finding 39 that single amino acids are not sufficiently transferable towards proteins, and (iii) an alloxan dimer (Cambridge Structural Database (CSD) identifier ALOXAN), which is a remarkable crystal (for a reason explained just below) featuring upfront in an important review 40 on intermolecular interactions.

Interacting quantum atoms (IQA)
Fig. 1 shows an example of a molecular system being partitioned into topological atoms.The picture shows a realm of atoms appearing as malleable boxes, or even ''bubbles'', negotiating energies and charges between them, and thereby forming a stable quantum system.As the nuclear positions change, so do these atoms' shapes and properties.It is clear that the topological atoms exhaust space: no overlap and no gaps.A remarkable property of crystalline alloxan (or pyrimidine-2,4,5,6-tetraone) is that hydrogen bonds do not appear in the crystal in spite of the two acidic NH entities for each molecule.
Instead, crystalline alloxan features 40 an array of short, almost orthogonal intermolecular CQOÁ Á ÁCQO contacts with distances of about 2.8 Å and CQOÁ Á ÁC angles in the range between 1551 and 1631.Two atomic interaction lines, each with a corresponding bond critical point, appear between the two monomeric units.
The IQA approach 10 follows the traditional QTAIM 9 formalism by partitioning a molecule into topological atoms.The IQA energies are calculated through the use of the first-order density matrix (non-diagonal) and the second-order (diagonal) density matrix, which allow the complete calculation of the Born-Oppenheimer energy of a molecule.To start, the intra-atomic and inter-atomic energies are given in eqn (1), where E IQA represents the molecular energy as obtained from an IQA partitioning, while E A IQA is the full atomic energy of atom A, E A intra is the intra-atomic energy of atom A (or ''self-energy''), while V AB inter represents the interatomic interaction energy between atom A and B. The intra-atomic contribution, E A intra , is further partitioned into: The intra-atomic energy comprises of the kinetic energy of the electrons (T A ), the electron-electron repulsive potential energy (V AA ee ) and the electron-nucleus attractive potential energy (V AA en ) within atom A.
The inter-atomic energy can be subdivided to give The inter-atomic energy comprises of the sum of four potential energies: the nucleus-nucleus repulsive energy (V AB nn ), electronnucleus attractive energy (V AB en ), nucleus-electron attractive energy (V AB ne ) and the electron-electron repulsive (V AB ee ) energy.The final energy can be divided further into a Coulombic (V AB Coul ) and an exchange-correlation (V AB XC ) contribution according to eqn (4), Here V AB XC is a combination of exchange (V AB X ) and correlation (V AB corr ) potential energies.
The calculation of the electron-electron potential V AB ee requires integration over atomic volumes of the second-order density matrix, which is computationally very expensive because of the evaluation of a six-dimensional integral, that is, three dimensions for each atom.The V AB ee energy can be expressed according to Go ¨rling-Levy theory 41 using perturbation theory as follows: The Coulombic term can be expressed as a functional of the electron density (first term of eqn ( 5)).However, the contributions from the exchange (second term) or correlation (third term) functions are more difficult as their analytical form remains unknown. 15Moreover, for DFT functionals such as B3LYP, the exchange term uses components from both the first-order and second-order density matrices. 15For these cases, the Kohn-Sham 42 second-order density matrix can be estimated from the first-order matrix, but will remain only approximate and differs from any analytical solution.Very recently, (dynamic electronic) correlation was incorporated 43 in IQA using closed shell coupled cluster theory.Compared to the exchange or Coulombic energies it is a minor quantity, at least for simple systems such as H 2 , N 2 , CO and H 2 O.
Returning to the conventional HF-compatible IQA, a rearrangement of the inter-atomic energies allows an electrostatic contribution and a covalent contribution to be separated and defined: where V AB elec (or V AB cl ) is the electrostatic contribution to the overall inter-atomic interaction energy between atoms A and B. The remaining exchange-correlation energy V AB XC , acts as a measurement of covalency between two atoms derived from the combination of inter-electron dependency and the Pauli Exclusion Principle.A final equation summarizes how interactions can then be typically analyzed through the combination of the electrostatic and covalent contributions to an interaction: At this point the brief review of the traditional IQA approach concludes.We will now address how both the traditional and the B3LYP-compatible IQA partitioning have been implemented into the computer program AIMAll. 11,44irst, a technical note is in place here.AIMAll slightly expands the IQA formalism explained above by additionally calculating so-called AA 0 inter-atomic energies.Here, A 0 represents every other atom in the system with the exclusion of atom A itself.Hence, A 0 is equivalent to a summation over all atoms B that are not equal to A, Note that recovering the total interaction energy within a molecule, which is formally P AB V AB inter , is more accurately calculated using V AA 0 inter as an intermediate quantity rather than where V AB inter is individually calculated for each atom pair.In theory, summing over B or working with the A 0 route should produce equal results but typically there are small differences in the energy values due to the nature of their separate calculation and associated errors (i.e. by grid quadrature or semi-analytical integration).In AIMAll, each of the interaction energy subcomponents (i.e.V AB en , V AB ne or V AB ee , see eqn (3)) can also be calculated in ''AA 0 mode''.AIMAll uses the AA 0 energies in the calculation of the total molecular energy instead of the summation over every interatomic AB interaction.Using AA 0 energies has two advantages: (1) a faster complete calculation of the IQA molecular energy partitioning albeit with reduced chemical insight, and (2) a more accurate recovery of the ab initio energy.Now we come to the actual incorporation of the explicit B3LYP functional into an IQA energy analysis.The strategy should be made clear upfront.The overall goal is to recover the total energy when using B3LYP, which has been glossed over in the literature so far.This goal is reached in this work by using the explicit functional only within a single atom, i.e. for the total atomic energy only.We will show that the functional cannot be used for an interatomic energy.Thus we have no choice but to adopt the Hartree-Fock-like expression for interatomic exchange energy but then using Kohn-Sham orbitals.
In the development of the B3LYP extension of IQA it is convenient to define energy contributions as total atomic energies, where ''comp'' can be XC, ne, en, nn or Coul.The quantity V A comp should not be confused with the previously defined intra-atomic energy (or ''self-energy'').This extra symbol is motivated by a future need (see below) to gather intra-and inter-atomic energy contribution, of a nature specified by ''comp''.
In previous versions of AIMAll, a B3LYP wavefunction would be treated identically to a HF wavefunction, which leads to a huge discrepancy between the total molecular energy and a summation of the various IQA energy components.In the older versions of AIMAll, both the intra-atomic and inter-atomic exchange energy is calculated from the HF exchange energy equation only: instead of from the (explicit) B3LYP exchange-correlation functional: 45 where a 0 = 0. is the Vosko-Wilk-Nusair local density approximation to the correlation functional.A ''pseudo-DFT/B3LYP'' compatible algorithm has quite recently been implemented in AIMAll in order to accurately recover a molecule's total ab initio energy.The qualifier ''pseudo'' refers to the ambiguity in the calculation of both interatomic and intra-atomic exchange.However, as explained below, a reasonable choice will be made to cope with this ambiguity.
Without extra thought or treatment, the introduction of a functional produces incorrect intra-atomic V AA XC and interatomic V AB XC energies (leading to incorrect V AA ee and V AB ee energies).As a result the total energy of the molecule is incorrectly calculated and one obtains a large energetic discrepancy.In order to obtain a correct result, the exchange-correlation functional needs to be explicitly incorporated.In a later AIMAll version (14.04.17 ref. 11), the total atomic exchange-correlation energy denoted V A XC (see eqn (10)) for the B3LYP functional was implemented.This implementation now allows the welldefined and correct calculation of the total V A XC component of V A ee for an atom A: In turn, the correct calculation of V A ee allows the complete recovery of an atomic energy E A IQA (see eqn (1)).When summing over all atoms, the ab initio energy is recovered from an IQA partitioning.It is important to keep in mind that, even with this implementation, standard functionals are constructed to include the interacting part of the kinetic energy, such that the atomic electron-electron repulsion in eqn ( 13) is inevitably contaminated by the residual kinetic energy term.
Interatomic exchange energies are chemically meaningful and therefore one must be able to calculate them, again in a DFT context.The division of the total atomic exchange-correlation into intra-(V AA XC ) and inter-(V AA 0 XC or V AB XC ) components is actually ambiguous.This ambiguity limits the correct extraction of interatomic terms but a reasonable choice or decision can be made.An ''amalgamated'' approach has been introduced to tackle this remaining limitation.This approach involves first calculating the inter-atomic exchange-correlation contribution (V AA 0 XC or V AB XC ) via This journal is © the Owner Societies 2016 the pure Hartree-Fock exchange equation only, but by inserting KS orbitals instead of HF orbitals, dr 1 dr 2 (14)   Note that the subscript X (rather than XC) emphasizes that, formally, we extract only the exchange part from the exchangecorrelation energy, using the HF framework.Secondly, this approach calculates the intra-atomic exchangecorrelation energy from the well-defined total atomic exchangecorrelation, V A XC,B3LYP , and the HF-defined inter-atomic exchange, V AA 0 X;HFðB3LYPÞ , or where V A XC,B3LYP is obtained by performing a three-dimensional integration over the volume of atom A, of E B3LYP XC defined in eqn (12).The calculation of V A XC,B3LYP is possible because the E B3LYP XC functional obviously only depends on the electron density.It is pivotal to appreciate the difference in treatment between inter-and intra-atomic exchange-correlation energies.Eqn ( 14) cannot be rewritten as a function of electron densities, and thus the B3LYP (which is a function of the electron density) cannot be introduced here.Thus, no interatomic exchange(-correlation) energy can be calculated from a 6D integral such as that in eqn ( 14) but then with the functional appearing in the integrand.
We also point out that a private communication with Dr Keith shows that V A XC,B3LYP is a transferable quantity as demonstrated from a series of tests he conducted on HÁ Á ÁH interactions in normal alkanes.
Finally, we note that in this proposed scheme, V AB X,HF(B3LYP) = V BA X,HF(B3LYP) , and that the total molecular exchange-correlation energy can be recovered as follows, The recovery of the total molecular energy from the intra-and inter-atomic contributions is made explicit in the following equation, where the quantity on the left hand side represents the molecular energy as recovered from the B3LYP-IQA analysis, which is to be compared with the original (unpartitioned) energy of the B3LYP wave function.
In this investigation we will analyse the effect of only using the HF-defined exchange equation (eqn (11)) in the calculation of interatomic exchange energies derived from KS orbitals.Note that, with the exception of the total atomic B3LYP exchangecorrelation functional energy, V A XC , all other intra-atomic and inter-atomic exchange-correlation energies are calculated incorporating only the HF exchange component without correlation.The appropriate notation will be applied throughout the analysis to maintain clarity.Subscript 'XC' is used for the total atomic exchange-correlation energy (V A XC ) along with the hybrid intra-atomic exchange-correlation energy (V AA XC ), and subscript 'X' is used for the inter-atomic exchange energy (V AB X or V AA 0 X ).

Computational methods
The program GaussView generates sensible geometries for NMA and TriGly, which were subsequently geometry-optimized with the GAUSSIAN09 46 program, in order to obtain a local energy minimum for NMA at HF/6-311G(d,p), and for TriGly at HF/6-31+G(d,p) level.The latter level of theory has been recommended 47 before for its favourable balance between accuracy and computational demand.Note that no frequency calculations were carried out because this work does not depend on the confirmation of molecular geometries being true minima.Subsequent wavefunctions and ab initio energies were calculated at HF, B3LYP and MP2 levels of theory (with the respective basis set) without re-optimization of the geometry.Note that the MP2 wavefunction did not feature in any IQA energy partitioning but was only used to obtain net atomic charges, necessary for proper quantifying of polarity (see Discussion, Section 4).The latter served in a discussion of bond polarity quantified at HF versus post-HF level.We emphasize that the molecular geometries were those obtained at HF level, and that they were never re-optimized at B3LYP or MP2 level, in order to guarantee a direct comparison of energy components.The geometry of the alloxan dimer was obtained directly from the CSD, without geometry optimization, and the corresponding wavefunctions were again calculated at HF, B3LYP and MP2 levels (with 6-31+G(d,p)).
The extensive analysis of NMA and TriGly is preceded by a succinct analysis of H 2 and LiH, at HF, B3LYP and full CI level using 6-311G(d,p).The latter level of theory is computationally extremely expensive and hence only possible for the smallest of molecules.Yet, this analysis is useful in setting the scene while making contact with the most advanced non-DFT post-Hartree-Fock method.
The AIMAll 11,44 software package (both versions 13.10.19and 14.04.17) was used to perform the QTAIM analysis and the IQA partitioning calculations for the HF and B3LYP wavefunctions.AIMAll (v 13.10.19)was used for the IQA partitioning of the HF and B3LYP wavefunctions, via a HF-only IQA algorithm, to investigate and compare the V AB X values.The B3LYP wavefunctions were then analysed using AIMAll (v 14.04.17) to check if molecular energy is recovered when using the implementation of the B3LYP atomic exchange-correlation functional (V A XC ) in the IQA algorithm.Note that both versions produce the same inter-atomic V AB X,B3LYP values since both follow the HF exchangeonly algorithm.AIMAll was also used to calculate QTAIM atomic charges.All default parameters (e.g.''auto'' quadrature grids ''Boaq'' and ''Briaq'') were used within AIMAll, except the encomp parameter (short for ''Energy Components''), which was set to ''4'' in order to carry out a full IQA partitioning calculation.
AIMStudio, which is a component of the AIMAll package, and in-house software called IRIS were used for visualization. 48,49Within AIMStudio, all default settings were used.Some non-default settings for IRIS were employed, i.e. wireframe for the surface and altered transparency.
Throughout the analysis no interaction energy has been doubly counted.Within the AIMAll format, all interaction energies are split into V AB and V BA , referring to any subscript (e.g.elec, inter, XC, X), and V AB V BA .Simply adding V AB X and V BA X gives the total exchange energy between atoms A and B, or V AB X,tot .For clarity, and to condense the data, only the V AB X,tot interaction will be reported in the figures.In line with the IQA formalism, V AB X,tot will simply be reported as V AB X .For each system a table will be presented containing the quantitative analysis of the interaction energies.Within each table, the analysis is divided into sub-sections according to which interatomic exchange energy is more stable, V AB X,HF or V AB X,B3LYP , where the subscript identifies the nature of the wavefunction used.Within each sub-section, maximum and mean difference percentages between the HF and B3LYP interatomic exchange energies will be presented.The total number of interaction cases for each 1,n interaction type within each sub-section is also included.Percentages were chosen for the analysis to allow a transferable measure that can be analysed, irrespectively of the magnitude of the energy value in question.Across the analysis, the interatomic exchange energies span a range of 10 À6 kJ mol À1 to just over 10 3 kJ mol À1 .

Materials
Before the extensive analysis of NMA and TriGly, the scene is set with a brief analysis of H 2 and LiH.NMA and TriGly are studied to investigate intra-molecular interactions in small and medium-sized systems.The alloxan dimer system adds another dimension by including intermolecular interactions.Fig. 2 shows the three systems; note that TriGly has 17 H atoms and 22 non-H atoms, while the alloxan dimer has 4 H atoms and 20 non-H atoms.Through individually analyzing every interaction in each of the three systems we aim to rationalize the difference between interatomic HF exchange energies, V AB X,HF , and interatomic energies obtained from a HF-treatment of a B3LYP wavefunction, denoted V AB X,B3LYP .

IQA molecular energy (E IQA ) vs. ab initio (E WFN )
Recent versions of the AIMAll program (version 14.04.17 and later) include the implementation of the explicit B3LYP exchange functional into the IQA partitioning.As a result, the IQA energies accurately recover the ab initio molecular energy (typically within o1 kJ mol À1 ).Table 1 shows the IQA energy recovery of each system, when using the incorrect HF exchange (AIMAll version 13.10.19)and then the correct B3LYP exchange (14.04.17).
Clearly, the ab initio energy of each molecular system is not recovered for B3LYP wavefunctions in AIMAll version 13.10.19because the values for DE run in the thousands of kJ mol À1 .The DE values are also seen to grow with system size.This growing error is due to the increasing number of V A ee energies being incorrectly calculated, which we expect to be more dramatic for non-hydrogen atoms.Alloxan has 3 non-hydrogen atoms more than TriGly (in spite of having fewer atoms in total), which probably explains why its energy error is larger than that of TriGly.
When the explicit B3LYP V A XC functional is used in the IQA partitioning, the ab initio molecular energy is accurately recovered.For TriGly and the alloxan dimer, the B3LYP recovery error is much smaller than that of HF.An pitiable IQA energy recovery of B3.8 kJ mol À1 was observed for the alloxan dimer, which is caused by poor integration errors (L(O) values) for 7 of the 8 carbon atoms in the dimer (see Table S1, ESI †).Extra integrations with non-default quadrature grids show that V AB X values are numerically stable, even when L(O) is quite large.Hence we deduce that another energy component, most likely the intraatomic energy, is sensitive to the magnitude of L(O).Next, the breakdown of the interatomic exchange energies, V AB X , will be investigated for each system.X energies for each interatomic interaction, together with the V AB X,HF = V AB X,B3LYP bisector.The figures are grouped into a triplet for each molecular system.Within the grouped figures, the first figure always incorporates all interatomic interactions in the system (66, 528 and 276 for NMA, TriGly and alloxan, respectively).Typically, the first figures of each triplet are dominated by the strongly covalent backbone interactions within the molecule, corresponding to V AB X energies in the region of hundreds of kJ mol À1 .The second figure of each triplet only identifies V AB X energies within the range of À60 to À1 kJ mol À1 ('medium-range').The third figure of each triplet only identifies V AB X energies within the range of À1 to 0 kJ mol À1 ('long-range').These two energy ranges can be loosely translated to strong and weak interatomic interactions.Finally, Fig. 4 and 7 show a direct overlay of |V AB X | values obtained at both HF and B3LYP level for the interactions in TriGly and the alloxan dimer, respectively.Fig. 6 is identical to Fig. S9 of the ESI, † and is shown in the main text for convenience and also due to its importance.
3.3.1 Small molecules and full CI.Table 2 compares three wave functions, in terms of interatomic exchange energy V XC (in kJ mol À1 ) and bond order d, at Hartree-Fock (HF), B3LYP and full CI (FCI) level of theory, using the same basis set of 6-311G(d,p) for all (which is also used for NMA below).
In general, the results are in line with previously obtained insight.V XC scales proportionally with d, as expected.In covalent systems (such as H 2 ), the Coulomb correlation localizes the electrons in the atomic basins, which means that d decreases and that the HF wavefunction is too delocalized.This is a general phenomenon that will be clearly and repeatedly recovered in NMA and TriGly.This exaggerated delocalization by HF is why the |V XC | value of H 2 is considerably smaller at the FCI level.In ionic molecules, such as LiH, the situation is reversed (although now the effect is much less pronounced).Here the mean field solution (i.e.HF) returns very well localized ions, while electron correlation transforms them somewhat back into neutral systems.This effect makes both d and |V XC | increase.The effect of electron correlation on ions is to change their electron density just a little.B3LYP captures this effect nicely, and the B3LYP and FCI results are in very good agreement.
In particular, d HH must be equal to 1, exactly (1.000000), both at the HF and B3LYP levels.This value does only depend on the one-determinant structure of the (pseudo)-wavefunction.This also explains why B3LYP fails to provide a good V XC value.Note that this means that other energy components are greatly affected by this inadequacy (kinetic, electron-nucleus, etc.) because the total B3LYP energy is close to the FCI one.The change in V XC from HF to B3LYP is small since the real change caused by electron correlation in this system (as opposed to that found in LiH) is not in the electron density, but in the second-order density matrix.The HF, B3LYP, and FCI densities in H 2 are quite similar but it is the pair behavior that is different.This can only be taken into account using a multi-determinant expansion, which provides smaller d and |V XC | values.
3.3.2NMA.NMA contains 12 atoms and hence (12 À 1) Â 12/2 = 66 unique interatomic interactions.NMA features regularly in biochemistry as a small and computationally inexpensive system that contains a peptide bond.NMA's interactions span a range of interaction types from 1,2 to 1,6.To elaborate on this notation, a 1,2 interaction consists of two neighbouring atoms separated by a single 'typical' covalent chemical bond; a 1,3 interaction is separated by two bonds, a 1,4 interaction is separated by three bonds, while an 1,n interaction is separated by n À 1 bonds.Thus, NMA serves as a good example for a small variety of interaction ranges encompassing the so-called bonded (1,2; 1,3 and 1,4 interactions) and non-bonded interactions (1,5 and higher).
Fig. S1-S3 (ESI †) compare the HF and B3LYP V AB X values in NMA, along with a bisector.The V AB X energies are distributed both above and below the bisector, indicating that increased stability may occur at either HF or B3LYP level.Fig. S1 (ESI †) shows all 66 interactions, and at this grand scale, a clear stabilization is perceived for only 3 interactions, which are all 1,2.To see more detail, Fig. S2 and S3 (ESI †) zoom in on smaller energy windows, which are À60 to À1 kJ mol À1 and À1 to 0 kJ mol À1 , respectively.Fig. S2 (ESI †) shows the vast majority of 1,3 and 1,4 interactions.Energies are scattered around the bisector, with almost twice as many data points in the |V AB X,B3LYP | 4 |V AB X,HF | half (which is hard to see, but clear from Table 3

below).
In Fig. S3 (ESI †), which shows only three 1,4 interactions and all 1,5 and 1,6 interactions, most data points now appear below the bisector, which means that B3LYP stabilizes these interactions compared to HF.This stabilization trend can also be seen as a 'shift' in the data points, in going from HF to B3LYP.This analysis benefits from more quantitative information about the distribution of the V AB X values at either side of the bisector.This information is provided in Table 3, which groups Table 2 A comparison between Hartree-Fock (HF), B3LYP and full CI (FCI) wavefunctions in terms of interatomic exchange energy V XC (in kJ mol À1 ) and bond order d.The molecules were geometry optimized only at HF level using at 6-311G(d,p) and the geometry kept constant for all three levels of theory.Note that V XC is actually V X for HF, and V X,HF(B3LYP) for B3LYP   3 shows only percentages, which characterize the degree to which the introduction of B3LYP changes the HF V AB X values.To understand how these percentages (maximum and mean) were calculated, it is instructive to focus on the 4 example interactions (N1-H8, N1-C9, N1-C2 and C2-O3) marked in red in Table S2 (ESI †).These interactions also appear in Fig. S1 (ESI †) where they lie in the |V AB X,B3LYP | 4 |V AB X,HF | region, that is, in the lower half, below the bisector.Table S2 (ESI †) makes clear that the largest relative energy shift, going from HF to B3LYP, is À111.8 kJ mol À1 (for N1-C2).Note that negative differences correspond to a stabilization caused by B3LYP.All B3LYP energies in this paper are always referred to HF, so it makes sense to work with the ratio of this energy difference to the original HF.Converting this ratio to a percentage, we define 100 Â |V AB X,B3LYP À V AB X,HF |/V AB X,HF as the ''absolute difference percentage''.The absolute difference percentage associated with C2-O3 then appears as 15.2% in Table 3 under the 1,2 entry referring to 4 interactions with |V AB X,B3LYP | 4 |V AB X,HF |.Table 3 also lists the mean difference percentage (i.e.9.6%), where the mean represents an average over these 4 interactions.
Table 3 shows that the interactions predominately stabilized by B3LYP are the 1,4 interactions and higher.The 1,2 and 1,3 interactions are mostly destabilized by B3LYP, in a ratio of about 2 to 1. Secondly, when |V AB X,B3LYP | 4 |V AB X,HF |, the mean difference percentage increases monotonically with n in 1,n over the whole range of n.However, when |V AB X,HF | 4 |V AB X,B3LYP |, there is no apparent trend linking difference percentage to the interaction type.Fig. 3 shows the shift in V AB X energies in going from HF to B3LYP for all 66 interactions in NMA.The logarithmic energy scale brings out the shift in terms of energy ratios.Fig. 3 confirms, at a glance, that the larger n in an 1,n interaction, the more pronounced the relative energy stabilization caused by B3LYP.This effect is consistent with the behaviour of the mean difference percentage in Table 3, for the stabilizing half on the right of this table.This effect may also be related to the fact that B3LYP is not an asymptotically correct functional, so it behaves pathologically at long-range.
It is pleasing to see that V AB X energies broadly fall apart into non-overlapping clusters, each cluster corresponding to an 1,n interaction, where n = 2, 3, 4, 5 or 6.The gap between the 1,2 cluster and the 1,3 cluster is particularly pronounced, while the 1,3 and 1,4 clusters are contiguous.The overall behaviour of V AB X energies as a function of internuclear distance is such is that the lowest values occur at the largest distances.The descending line in black captures this behaviour coarsely, with a correlation coefficient r 2 of 0.85, referring to the B3LYP data only.This broadly respectable correlation does not attempt to spot any non-linearity in the data, nor does it detect a fine structure emerging by grouping V AB X values per interaction pair AB.
Closer inspection of the clusters reveals some energetic anomalies, that is, interactions appearing in unexpected clusters.For example, four 1,5 interactions (green) turn up far outside their cluster.Fig. 3 marks these as O3-H12, O3-H11, O3-H10 and H5-H8.The first two interactions of this set appear near the overall (B3LYP) correlation line (in black).Hence they are not anomalous in terms of how a V AB X energy varies with AÁ Á ÁB distance.On the other hand, O3-H10 is anomalous, not just by connectivity ''distance'' (i.e.1,5) but by actual (geometrical) distance.Indeed, this interaction is curiously strong for its relatively large internuclear distance.The reason for this observation is not clear.Inspired by earlier work 25 that showed that V H-H X values are anomalously large when part of a planar H-C-C-H arrangement in alkanes, it is tempting to invoke the near-planarity of the H10-C9-N1-C2-O3 fragment (where the H10-C9-N1-C2 dihedral angle is À1631 (off by only 171 from 1801) and the C9-N1-C2-O3 dihedral angle is 41).Finally, note that V N1-O3 X is the strongest interaction in the 1,3 cluster, quite far to the right of the black line.This enhanced delocalization across the C-N peptide bond to the oxygen is reminiscent of the resonance canonical imparting double bond character to the CN interaction.

TriGly (GlyGlyGly).
There are 528 interatomic interactions in this oligopeptide (Fig. 2) ranging up to 1,15 interaction  types.Ignoring a hydrogen atom each time, NMA occurs four times in GlyGlyGly, which therefore marks a major upscaling in the size of our analysis, while also establishing links between the NMA and GlyGlyGly.Table 4 shows that the interactions that B3LYP predominately stabilizes, are the 1,4 interactions and higher, as before for NMA, but this effect tails off at 1,13 and even reverts at 1,14.When lumped together, the 1,2 and interactions are neither stabilized nor destabilized by B3LYP.In other words, 46 interactions are less stable under B3LYP while 40 are more so.As in the case of NMA, the mean difference percentage increases monotonically with n (in 1,n) up to 1,8 interactions, again when |V AB X,B3LYP | 4 |V AB X,HF |.Thus, in TriGly, the range of the monotonic trend observed up to the 1,6 interactions in NMA, is extended to 1,8.However, at this point, the mean difference percentage decreases (but not monotonically) until the 1,15 interactions are reached.When |V AB X,HF | 4 |V AB X,B3LYP |, there is still (as for NMA) no apparent trend linking difference percentage to the interaction type, throughout the whole range of interactions.
On a technical note, such large difference percentages for long-range 1,n interactions could be a cause for concern in the analysis above.However, the reliability of the values for V X should not be called into question, as we now explain.Table S3 (ESI †) shows how a variation in the outer angular integration grid (within AIMAll) affects the V X energies for the very high 1,n interactions for two example interactions in TriGly: C22Á Á ÁH27 (1,14) and H24Á Á ÁH28 (1,15).A change in the outer angular quadrature integration grid changes the integration error L(O) calculated for each atom.A computationally less expensive grid, for example, 'Low (Lebedev)' will typically 50 result in a higher absolute L(O) value for a given atom, compared to that of 'Sky High (Lebedev)'.Table S3 (ESI †) shows that the variation of the integration grid does cause a significant change in L(O) but importantly V X remains very stable throughout.Such small changes in the values of V X , compared to that of L(O), allow us to trust the V X values, even for very high n in 1,n.This stability of V X suggests that a large L(O) error affects another IQA energy component more, presumably E A intra .Fig. 4 shows the effect on V AB X energies of introducing B3LYP, compared to the original HF values, for all 528 interactions in TriGly computed at the HF/6-31+G(d,p) and B3LYP/6-31+G(d,p) level of theory.As in Fig. 3, it again makes sense to express the overall behaviour of V AB X energies as a function of internuclear distance in a coarse, linear way.However, here we draw two lines, one for HF and one for B3LYP, instead of the single black line in Fig. 3.The descending line in black captures this behaviour, admittedly coarsely, with a correlation coefficient r 2 of 0.91 for both HF and B3LYP.It is clear that the cloud of B3LYP energies (in red) occurs right of the cloud of HF energies (in blue), which means that most interatomic interactions are indeed more stable at B3LYP level.It is interesting to again analyze anomalies in the same manner as in Fig. 3.It is best to truncate the analysis at 1,6 interactions because energies of about 1 J mol À1 or less are not really significant (unless hundreds of them are considered as a group and their energies simply added, as allowed by the additive nature of V X ).
Fig. 5 shows again the broad behaviour seen before in NMA (Fig. 3): the further two atoms are apart, the less they exchange (i.e.reduced delocalization).In other words, well separated  atoms have a small V AB X energy.This behaviour is typical for saturated systems, such as this tripeptide, which lacks aromatic amino acids or other extensive electronically delocalized fragments.Secondly, as seen in NMA, the various 1,n interactions do cluster as well, albeit with more overlap compared to NMA.The most substantial overlap takes place between the 1,5 (green) and 1,4 (red) cluster.Such overlap is not surprising because TriGly is a more complex system, which contains a wider variety of interactions than found in NMA.In particular, because TriGly is curled (see inset) there are many instances where the actual (geometrical) distance is smaller than the distance based on bond connectivity (i.e.labelled by 1,n).This mismatch between geometrical and connectivity-based distance shows itself in 1,6 interactions where, for example, the H10Á Á ÁH19 distance is only 2.68 Å.The corresponding V AB X (B3LYP) energy is a respectable 0.4 kJ mol À1 , whereas the mean V AB X energy of all 1,6 interactions is at least an order of magnitude smaller.
The next set of interactions of interest is the same as in NMA: chemically meaningful quantities such as V AB X are anticipated to be transferable.We expect the NÁ Á ÁO interactions associated with the peptide bond to be very similar in energy value.QCT is known for providing highly transferable properties.Indeed, the V NÁ Á ÁO X energy averaged across all four peptide bonds in TriGly and the one peptide bond in NMA, is À113.4 kJ mol À1 with a standard deviation of only 1.0 kJ mol À1 at B3LYP level.The standard deviation at HF level is only 0.2 kJ mol À1 , which is a remarkable testimony of transferability.As expected, the V NÁ Á ÁO X energies again (see NMA, Fig. 3) appear at the strong end of the 1,3 cluster (Fig. 5, at the extreme right of blue data points, marked as ''OQC-N group'').Finally, we mention that the 1,5 interaction H5Á Á ÁO8 also appears at the strong extreme of the 1,5 cluster, but energetically this interaction is less anomalous than any OÁ Á ÁN interaction because V H5Á Á ÁO8 X lies very near the black line.In contrast, the 1,6 interaction O8Á Á ÁO14 lies quite a bit off the black line, and marks an unusually strong stabilization of about 6 kJ mol À1 by ''through-space'' oxygenÁ Á Áoxygen delocalization.This observation sets the scene for exploring another type of interaction that received a great deal of attention in the last decade, called the np* interaction.
TriGly is an ideal system to make contact with the research theme of Raines and co-workers on so-called np* interactions in proteins. 51The interaction occurs between protein backbone amides, and is characterized as arising from the delocalization of a lone pair of electrons (n, in fact, the p-rich lone pair) of an amide oxygen, to an antibonding orbital (p*) of the subsequent carbonyl group.This research group has elevated this type of interaction to one that may be as important as hydrogen bonding, which they characterize as an n ss* interaction.Here, there is an s-rich lone pair of an amide oxygen delocalizing with the antibonding s* orbital of the N-H bond of an amide bond, but now four amino residues further down the protein backbone.Although the np* interaction can continuously get stronger, Raines et al. propose that the interaction should unmistakably take place if the distance d between O B and C C in C A QO B Á Á ÁC C QO D is smaller than 3.2 Å, and if the angle y, defined as O B Á Á ÁC C QO D , lies between 991 and 1191.Note that the quantum mechanical interpretation of Raines et al. is couched in the language of Natural Bond Analysis (NBO), 7 which suffers from practical issues of computational stability and lack of transparency.However, IQA is alternative candidate to express delocalization effects but now in an orbital invariant way, as well as computationally stable and transparently.In TriGly, the C A QO B Á Á ÁC C QO D system that comes nearest to obeying those geometrical criteria is C 6 QO 8 Á Á ÁC 2 QO 14 , where the distance d is 3.4 Å and y is 871.It turns out that that the V C2Á Á ÁO8 X energy is 2.8 kJ mol À1 .This energy is close to a value quoted by Raines et al., still in kcal mol À1 , and amounting to 0.5 in that energy unit, obtained after an NBO analysis of an a helix.Having said this, the 1,5 interaction C2Á Á ÁO8 (marked in Fig. 5), occupies a rather unremarkable place in the green cluster or with respect to the black line.A likely reason is that the C 6 QO 8 Á Á ÁC 2 QO 14 fragment falls out of the range proposed by Raines et al., strictly speaking.
3.3.4Alloxan dimer.This system is the first example of a van der Waals complex, which is of interest to the current analysis because of its surprising total lack of hydrogen bond as discovered by Bolton 52 in its stable crystalline structure.However, Dunitz and Schweizer ask if alloxan is really a ''problematic'' structure 53 although it has been regarded as such for about 40 years.They write that ''It (''alloxan'' red.) is held together mainly by whatever factors contribute to the cohesive energies. . . the CQ Q QOÁ Á ÁCQ Q QO type of interaction would seem to play an important part''.Unpublished work from our lab shows that the electrostatics between topological atoms plays an important role in the stability of the complex.Here we only discuss the exchange energy component.The investigation of the intermolecular interactions allows us to compare with the intra-molecular interactions, and present a complete analysis of the V AB X energies.The alloxan monomer is a heterocyclic planar system.Within each alloxan monomer, intramolecular interactions only reach to the maximum value of 1,6.When added to the intermolecular interactions, a total interaction count of 24 Â (24 À 1) = 276 is obtained, in the knowledge that each alloxan monomer consists of 12 atoms.In the current work, only one of the two known alloxan dimer crystal orientations 53 was investigated (see Fig. 2), the other being a dimer held by bifurcated hydrogen bonds.For the purpose of this analysis, a single dimer arrangement is sufficient.The monomers are roughly T-shaped (with a tilt), with two carbonyl groups of one monomer pointing towards the plane of the other monomer, approximately bisecting this plane.
Fig. S7-S9 (ESI †) again compare the HF and B3LYP V AB X values, but now in alloxan.Fig. S7 (ESI †) shows all 276 interactions, and at this grand scale, clear stabilization materializes for only a dozen or so interactions, which are all 1,2.The introduction of B3LYP renders the strongest interactions stronger.Fig. S8 (ESI †) (À60 to À1 kJ mol À1 ) shows the vast majority of 1,3 and all but one of the 1,4 interactions, along with very few 1,5 interactions, as well as up to a quarter of inter-alloxan interactions.Energies are scattered around the bisector, with about 10 times as many data points in the |V AB X,B3LYP | 4 |V AB X,HF | half.Fig. S9 (ESI †) repeatedly shows a larger shift towards stabilization by B3LYP.The interactions from 1,4 to 1,8 follow trends identical to those observed for NMA, with greater stability favoured in the B3LYP wavefunction.Fig. 6 (which is identical to Fig. S9, ESI †) covers the range 0 to À1 kJ mol À1 , which corresponds to 1,5; 1,6 and the inter-alloxan interactions.Here, 83% interactions are stabilized by B3LYP.The predominately stabilizing trend caused by B3LYP also appears in purely intermolecular (i.e.alloxanÁ Á Áalloxan) interactions.In total, there are 144 of these interactions, which amounts to 52% of the complete total of 276 interactions.Table 5 shows these intermolecular interactions as a separate entry.There are now 27 cases for which |V AB X,HF | 4 |V AB X,B3LYP | and 117 for which |V AB X,B3LYP | 4 |V AB X,HF |.The distribution of cases over these two stabilization regimes can be rationalized as follows.From the TriGly data we learnt that interactions are more stable at HF level if they are (very) long range (or 1,2 and 1,3 but this very short range is not relevant here).Hence, we expect the average internuclear distance for the 27 cases to be large.Indeed, it turns out to be 6.5 Å as shown in Table S6 (ESI †).The value of 6.5 Å falls in between 6.3 Å and 6.8 Å, which are the respective average internuclear distances for 1,8 and 1,9 interactions in TriGly.Fig. 7 demonstrates the energetic stability that B3LYP brings to the vast majority of interactions in the alloxan dimer: the red cloud (B3LYP) is indeed predominately shifted to the right compared to the blue cloud (HF).This figure is the counterpart of Fig. 4, and again shows two lines, one for HF and one for B3LYP, with almost identical correlation coefficients r 2 of 0.85 and 0.84, respectively.Note that the x-axis is logarithmic hence |V AB X | must be used.It is clear that the shift becomes more pronounced as the interatomic distance increases and the exchange energy decreases.The shift is more pronounced here than in TriGly (see Fig. 4).
Finally, Fig. 8 shows again the broad behaviour seen before in NMA and TriGly (Fig. 3 and 5) but now for the alloxan dimer.The black line represents the overall relationship between V AB X and internuclear distance.The corresponding correlation coefficient r 2 is 0.84.This time the clusters represents the exchange energy patterns of the alloxan monomers only, and show again minimal overlap.Four intramolecular OÁ Á ÁO interactions stand out by their unusual strength: O7Á Á ÁO8, O7Á Á ÁO12 in one alloxan, and O19Á Á ÁO2O, O19Á Á ÁO24 in the other.The average  V OO X B3LYP energy is À30.8 kJ mol À1 with a pleasingly small standard deviation of 0.9 kJ mol À1 (original data in Table S5, ESI †).These four OÁ Á ÁO interactions appear outside of their 1,4 cluster.This observation is reminiscent of that made in TriGly where the 1,6 interaction O8Á Á ÁO14 appeared far outside its cluster.
Fig. 8 also marks inter-alloxan interactions (in grey-brown) having values of |V X,B3LYP | 4 0.01 kJ mol À1 .These intermolecular interactions mainly populate the weaker end of the energy spectrum: they start in the utmost left region, which is weaker than the intramolecular 1,6 interactions (flat disk), and extend over the 1,5 interactions (green disk) well into the 1,4 interactions (red disk).In fact, two intermolecular interactions even outdo the 1,4 interactions in terms of strength: O7Á Á ÁN17 and O12Á Á ÁC14 (marked in Fig. 8, and ignoring the 4 intramolecular OÁ Á ÁO interactions outside the 1,4 disk).They are the strongest intermolecular interactions, and remarkably correspond to the two interatomic interaction lines found between the two monomers.The O12Á Á ÁC14 interaction materialises between two CQO groups, and can hence be linked to Dunitz and Schweizer's comment on the important role of CQOÁ Á ÁCQO type interactions played in the stability of the crystal.However, a convincing argument on the stability of the alloxan dimer can only be made when considering all IQA energy contributions.

Discussion
The interatomic exchange energy results provide several points to consider more deeply.These points will be addressed in turn, starting with a brief comment on the comparison between delocalization indices and V AB X .A systematic increase in energetic stability is seen between pairwise 1,4 bonded interactions and higher non-bonded interactions.This trend was observed unequivocally until B1,14 interactions.However, interactions of 1,6 and higher, in general, had interatomic exchange energies that fell below the integration error of the atoms involved, eventually reaching values as little as B1 Â 10 À6 kJ mol À1 in magnitude.Such interactions would be considered negligible for any application.Hence, up to B1,6 interactions, served as an appropriate cut-off for interaction types.
A comparison can be made between the V AB X results presented and the trend previously observed in measuring pairwise electron delocalization, via the delocalization index (DI).The consistent overestimation of DI values, which has been observed before, matches the current results for non-bonded interactions (1,5 and higher).This overestimation provides evidence that some Coulomb electron-electron correlation (which is included in the one-electron density) is partly included in the pairwise exchange results, despite the HF-formalism not directly addressing it.This is a well-known phenomenon, which allows a simple explanation in terms of the different behaviour of Kohn-Sham and Hartree-Fock orbitals.KS orbitals are more localized compared to HF orbitals, which is common in largely covalent interactions.This increased localization leads to smaller DIs.Because exchange energies are proportional to DIs, they will be smaller as well (in absolute value).The opposite will be true when KS orbitals become more strongly delocalized than their HF counterparts.In any case, the exact behaviour will never fully be included without the electron-pair density being defined within the HF-formalism.
The 1,2 and 1,3 interactions are neither stabilized nor destabilized by B3LYP.To rationalize this fact, it helps to consider work of Poater et al. 36    A more detailed analysis is also possible that allows for a physicochemical interpretation of the systematic differences found in our data beyond the general discussion already presented.Tables S2, S4 and S5 (ESI †) gather all the (1,n) n = 2,3,4 pairs for NMA, TriGly and the alloxan dimer, respectively, listing their V AB X,HF and V AB X,B3LYP values, their difference, and absolute difference percentages.In order to shed some light on them we take into account the major factors that influence V AB X .We will use orbital arguments, which are easier to grasp, although fully orbital invariant insights may also be given.At the HF or KS levels, So the first, and probably most important factor in determining V AB X , is the degree of delocalization of the orbitals involved in a given AB interaction.It is easy to understand from eqn (18) that only orbitals delocalized between atoms A and B will contribute non-negligibly to V AB X .Inclusion of electron correlation effects, even at DFT level, tend to localize orbitals involved in direct covalent links, and to delocalize their long-range tails.For a given AB pair, other factors involve the size of the topological atoms, and the amount of electron charge contained in them.For instance, it is clear that in 1,2 interactions, a shift in the position of the inter-atomic surface (separating A and B) towards the ideal mid-partition will increase V AB X .In summary, homopolarity works in favour of covalency.Conversely, heteropolarity works against covalency.Similarly, V AB X is usually dominated by the diagonal i = j term in eqn (18), the so-called self-interaction contributions, which behave like products of densities.An increase in the electron population of an atom will necessarily increase these self-interaction terms, adding to a larger V AB X .There are other effects that the above expression (only valid for single-determinant expansions) will never capture.There are other effects that the above expression (only valid for singledeterminant expansions) will never capture.At any standard HF or KS DFT level, a hydrogen molecule will be described by one Slater (pseudo)determinant consisting of a doubly occupied totally symmetric s g orbital, independently of the internuclear distance.Its DI will be exactly equal to 1, and its V AB X value will be very large, not vanishing at dissociation.This is the well-known single determinant dissociation problem of H 2 rephrased in terms of DI's or V AB X 's.A proper wavefunction correlated method, as simple as a CAS[2,2] calculation, will mix at least two determinants, leading to a destructive interference that will make V AB X vanish at the large distance limit.We will now show how a combination of these effects rationalizes our findings.
It is useful to introduce net atomic charges (denoted q) into this discussion.Tables S7-S9 (ESI †) provide the net charges of all atoms in NMA, TriGly and the alloxan dimer, respectively.These data show that B3LYP produces net accurate charges that are very close to those generated by MP2, which models correlation effects independently from DFT.
As noticed above, our data favour a clearly biased distribution for 1,4 and higher interactions.When |V AB X,B3LYP | 4 |V AB X,HF |, the sum of the net charge on A and on B at the B3LYP level (i.e.q(A) + q(B)) is lower than the sum at the HF level across all 1,n interaction types.In other words, the total electron population contained in the B3LYP (topological) atoms is larger.This is in agreement with the comments on the role of self-interaction posed in the above paragraph, and larger atomic electron populations will be, in general, accompanied by larger |V AB X,B3LYP | values, particularly when the implied extra electrons are delocalized between A and B, such that the product c i 2 (r 1 )c i 2 (r 2 ) contributes significantly to V X (see eqn (18) and recall that we integrate electron 1 over A and electron 2 over B).Almost all the exceptions to this rule may also be rationalized.For instance, in the C2-C10 interaction in the alloxan dimer, q(A) + q(B) at B3LYP is lower, because these carbon atoms participate in very polar CQO interactions at the HF level.However, |V AB X,B3LYP | o |V AB X,HF |.In this case, the smaller net charge of the B3LYP carbon atoms is due to electrons provided by the neighbouring carbonylic oxygens, which do not affect the electrons participating in the C2-C10 interaction.The latter are still much more delocalized at the HF than at the B3LYP level, giving rise to a larger V X in the former case, as the following paragraph explains.
Similarly, across all 3 systems, all the 1,2 bonded interactions that showed increased stability at the HF level were C-H or C-C interactions.In the case of C-H, the larger HF exchange stems from a spurious polarization giving rise to negatively charged hydrogen atoms, caused by a too delocalized C-H bonding orbital.However, B3LYP recovers the correct polarity by shifting the interatomic surface (or equivalently localizing the bonding orbital) and thereby decreasing the exchange.In the case of C-C, we face the well-known fact that too delocalized HF orbitals lead to exaggerated covalency, counteracted again by localization in the DFT description.We should add that it is a general phenomenon that inclusion of electron correlation considerably decreases both the DI and the exchange-correlation energy in 1,2 homoatomic bonds.
An analysis of the 1,3 interactions along the previous lines also explains the stability of the interactions at each level of theory.The 1,3 interactions that are more stable at HF, consisted of a wider variety of elements (or types of topological atoms, e.g.CÁ Á ÁN, CÁ Á ÁO, NÁ Á ÁH) across the three case studies.Typically, at HF level these interactions typically exhibit a lower sum (not absolute value) of the net charge on A and on B. In some interactions (e.g.NÁ Á ÁH), the polarity of the H atoms involved, changed according to whether the HF or B3LYP wavefunction was used.At HF, a number of H atoms had exaggerated charges and were often anionic.Some of these atoms exhibited a change in polarity, becoming cationic at B3LYP.The previous anion-anion behaviour would lead to an increase in the |V AB X | at HF level.The small number of 1,4 cases that exhibited greater stability at HF were solely due to a decreased sum, q(A) + q(B), at this level.Typically, these cases were limited to N-H, N-N and H-H.
We should add, in passing, that obtaining intermolecular exchange contributions from KS orbitals is not new, because it is actually a standard procedure in the SAPT-DFT method originally proposed by Williams and Chabalowski, 54 and later improved by Hesselmann and Jansen. 55This method has become an accurate, yet cheaper alternative to the very expensive symmetry-adapted perturbation-theory calculations.
Overall, this analysis validates the current expansion of IQA to B3LYP level of theory without reason for concern, in spite of the absence of an explicit second order reduced density matrix.

Conclusion
To the best of our knowledge, this study is the first to calculate and study interatomic exchange energies V AB X by using an explicit density functional, in this case that of B3LYP.As such, the IQA energy decomposition has been satisfactorily extended to DFT.Indeed, we recover the original (unpartitioned) energy of the system from the IQA energy contributions.The correct version (14.04.17 and higher), of the computer program AIMAll contains a pseudo-DFT IQA algorithm that makes this recovery possible within B0.7 kJ mol À1 of the original ab initio B3LYP energy.
We studied the influence of B3LYP on V AB X , compared with the HF values, for three systems of biochemical interest: N-methylacetamide (NMA), doubly capped tripeptide GlyGlyGly (TriGly) and an alloxan dimer.Over the three systems, a grand total of 870 V AB X energies for interactions ranging up to 1,15 interactions were carefully analyzed.The 1,6 and higher interactions, with V AB X energies corresponding to B0.01 kJ mol À1 , were disregarded because they are negligible.The introduction of B3LYP invariably stabilized 1,4 and higher interactions, showing a consistent energy shift towards increased exchange energy, in terms of absolute value.Such a shift in stability is almost always due to an decrease in the sum of net atomic charges (q(A) + q(B)) caused by the use of Kohn-Sham orbitals.
However, for 1,2 and 1,3 bonded interactions, the V AB X energies were not always shifted in the same direction, i.e. sometimes the HF energies were more stable than the B3LYP energies.The exchange stability for each 1,2 interactions can be understood by considering a combination of the (de)localization of the KS or HF orbitals, the position of the interatomic surface and the sum of net atomic charges (inferring polarity) of the involved atoms.Through the orbital description, DFT methods have some correlation incorporated, resulting in more localized orbitals than HF.Only C-C and C-H interactions are the 1,2 interactions with greater V AB X energies using HF orbitals (over KS orbitals).Here, HF over-delocalized their atomic orbitals leading to greater exchange in the interaction.Contrastingly, the increased heteropolarity, experienced by some of the 1,2 interactions (e.g.N-H, N-C, C-O) through a HF orbital description, counteracts the lower sum of net atomic charges, leading to less interatomic exchange at HF but more so at B3LYP.All remaining 1,2 interactions had a reduced sum of net charges at B3LYP, resulting in greater interatomic exchange stability.The 1,3 interactions can be explained similarly to those of the 1,2 interactions.
Plots of V AB X versus internuclear distance shows that the vast majority of interactions emerge in well-defined clusters, each associated with a bond-connectivity 1,n.Similar AÁ Á ÁB pairs yield very similar V AB X values, a hallmark of high transferability.Some anomalously strong interactions occur, shedding new light on patterns of exchange stability beyond that of the traditional covalent bonds.
In keeping with previous views and work on delocalization indices, HF-like B3LYP interatomic exchange energies (V AB X ) also provide chemically relevant values, but now energies.This realization occurs despite only first-order (and not electronpair) electron density correlation being accounted for.Overall, the extension of IQA to B3LYP can thus be justified.

3. 3
Fig.S1-S9 (ESI †) contrast the HF and B3LYP HF-calculated V AB X energies for each interatomic interaction, together with the V AB X,HF = V AB X,B3LYP bisector.The figures are grouped into a triplet for each molecular system.Within the grouped figures, the first figure always incorporates all interatomic interactions in the system (66, 528 and 276 for NMA, TriGly and alloxan, respectively).Typically, the first figures of each triplet are dominated by the strongly covalent backbone interactions within the molecule, corresponding to V AB the V AB X values according to which side of the bisector they lie, i.e. either |V AB X,HF | 4 |V AB X,B3LYP | (upper half in Fig. S1-S3, ESI †) or |V AB X,B3LYP | 4 |V AB X,HF | (lower half).Table

Fig. 3
Fig. 3 Logarithmic plot of |V AB X | versus interatomic distance for all 66 NMA interactions at HF/6-311G(d,p) and B3LYP/6-311G(d,p) level of theory.Interaction types ''1,n'' have been encircled where possible, showing the clustering as well as the outliers.
Fig. S4-S6 (ESI †) again compare the HF and B3LYP V AB X values, but now in TriGly (counterpart of Fig. S1-S3, ESI †).Fig. S4 (ESI †) shows all 528 interactions, and at this grand scale, clear stabilization materializes for only a dozen or so interactions, which are all 1,2.Fig. S5 (ESI †) (À60 to À1 kJ mol À1 ) shows the vast majority of 1,3 and 1,4 interactions, along with about a third of the 1,5 interactions.Energies are scattered around the bisector, with 2.6 times as many data points in the |V AB X,B3LYP | 4 |V AB X,HF | half.Fig. S6 (ESI †) repeatedly shows a larger shift towards stabilization by B3LYP.The interactions from 1,4 to 1,8 follow trends identical to those observed for NMA, with greater stability favoured in the B3LYP wavefunction.

Fig. 5
Fig. 5 Logarithmic plot of |V AB X | versus interatomic distance for TriGly interactions up to 1,6, at HF/6-31+G(d,p) and B3LYP/6-31+G(d,p) level of theory.Key outlying V X values have been labelled in both the plot and the inset molecule image.The black line shows the overall correlation of the B3LYP energies, with a correlation coefficient r 2 of 0.91.

Fig. 6 V
Fig.6V AB X energies for long-range interactions (0 to À1 kJ mol À1 ) in the alloxan dimer, from both HF and B3LYP wavefunctions.Note that this figure is repeated in the ESI † as Fig.S9because there it is part of a triplet of plots on alloxan.
regarding a similar analysis on d(A,B) rather than on our V AB X values.From their study of ethane and diborane, one concludes that d(A,B) values of the 1,2 and 1,3 interactions involved can either increase or decrease.Given the proportionality between |V AB X | and d(A,B), mentioned in the Introduction, we can conclude that the results of Poater et al. support our own findings, stated at the beginning of this paragraph.From this parallel behaviour of d(A,B) and V AB X results, we conclude that some first-order density correlation must also being included in V AB X .

Fig. 7
Fig. 7 Logarithmic plot of V AB X versus interatomic distance for all 276 interactions in the alloxan dimer at HF/6-31+G(d,p) and B3LYP/ 6-31+G(d,p) level of theory.

Fig. 8
Fig. 8 Logarithmic plot of |V AB X | versus interatomic distance for alloxan interactions with |V X,B3LYP | 4 0.01 kJ mol À1 , at HF/6-31+G(d,p) and B3LYP/ 6-31+G(d,p) level of theory.Interaction types ''1,n'' have been encircled where possible, showing the clustering as well as the outliers.Bracketed interaction labels indicate the interaction also forms an AIL.The black line shows the overall correlation of the B3LYP energies, with a correlation coefficient r 2 of 0.84.

Table 1
IQA molecular energy (E IQA ), ab initio energy (E WFN ) and IQA residual error (DE) of NMA, TriGly and the alloxan dimer, at HF and B3LYP theory levels, using AIMAll versions 13.10.19and14.04.17.All energies are in kJ mol À1 .Note that, in the case of B3LYP, E IQA (from AIMAll version 14.04.17) is the same energy as the left hand side of eqn(17)or E molec This journal is © the Owner Societies 2016

Table 3
Stabilization and destabilization of V AB X energies in going from HF to B3LYP for all 66 interatomic interactions in NMA.Maximum and mean absolute difference percentage (between and HF and B3LYP) sorted by interaction type ''1,n'' and by whether |V AB X,HF | 4 |V AB X,B3LYP | or |V AB X,B3LYP | 4 |V AB

Table 4
Stabilization and destabilization of V AB X energies in going from HF to B3LYP for all 528 interatomic interactions in TriGly.Maximum and mean absolute difference percentage (between and HF and B3LYP) sorted by interaction type ''1,n'' and by whether |V AB X,HF | 4 |V AB X,B3LYP | or |V AB X,B3LYP | 4 |V AB X,HF |

Table 5
Stabilization and destabilization of V AB X energies in going from HF to B3LYP for all 276 interatomic interactions in alloxan.Maximum and mean absolute difference percentage (between and HF and B3LYP) sorted by interaction type ''1,n'' and by whether |V AB X,HF | 4 |V AB X,B3LYP | or |V AB X,B3LYP | 4 |V AB X,HF |.Just over half of the interactions (144 in total) are of intermolecular nature This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 20986--21000 | 20997 This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 20986--21000 | 20999