Energetics of Electron Pairs in Electrophilic Aromatic Substitutions

Abstract

The interacting quantum atoms approach (IQA) as applied to the electron-pair exhaustive partition of real space induced by the electron localization function (ELF) is used to examine candidate energetic descriptors to rationalize substituent effects in simple electrophilic aromatic substitutions. It is first shown that inductive and mesomeric effects can be recognized from the decay mode of the aromatic valence bond basin populations with the distance to the substituent, and that the fluctuation of the population of adjacent bonds holds also regioselectivity information. With this, the kinetic energy of the electrons in these aromatic basins, as well as their mutual exchange-correlation energies are proposed as suitable energetic indices containing relevant information about substituent effects. We suggest that these descriptors could be used to build future reactive force fields.

1. Introduction

Computational and theoretical chemistry (CTC) has undoubtedly come of age. Advances in both electronic structure methods and computer technology now allow researchers to perform simulations of realistic systems with outstanding accuracy within affordable times, so that CTC is already playing a crucial role in fields like big pharma industry or materials design [1,2,3,4,5,6]. Machine learning (ML), on the other hand, has fully revolutionized the field [7,8], and problems which were once thought to be out of reach in a foreseeable future, like protein folding, are close to being solved [9]. Given this favorable scenario, it should never be forgotten that the simulation of complex systems, or the training of deep neural networks, rely at one point or another on human understanding of an underlying energy model, be it a well parameterized force field (FF) in the first case, or a specific neural network (NN) in the second. Many successful FFs or NNs proposed in the CTC context share a common structure rooted in some of the most cherished chemical concepts: atoms that form molecules through well-defined spring-like chemical bonds, which can be formed or broken [10]. Although the formation and rupture of these entities that chemists call ‘bonds’ lies at the core of chemical intuition, the formal construction of atomistic models has usually followed, perhaps as a result of historical precedence, the physicist’s viewpoint, based on a many body expansion of the energy landscape. Early atomistic models [11] already introduced the idea that pairwise additive pair potentials were essentially enough for simulation purposes, with smaller three body corrections that could be added ad hoc, if necessary. This mantra was inherited in the first wide purpose force fields and has remained practically untouched to date.

The mismatch between the bond-centric chemical and the atomistic points of view has led to several aberrations, which may well be exemplified by the concept of atom types. When atoms are used as the FF building blocks, specific parameterizations for different coordination or bonding environments are needed. Up to ten different types of nitrogen atoms exist in the generalized Amber FF (GAFF) [12], for instance. The situation in ML-based methods is not better, since NNs normally inherit the atom types logic, although the actual number of types is not fixed, but inherited from the information provided for the training of the network. It is clear that much of this zoo could be avoided if electron pairs, not atoms, were selected as the energy model building blocks, as it is electrons (i.e., bonds) that rearrange in the course of a reaction. Moreover, the direct use of electron-pair units (which behave like bosons) may decrease the complexity of the parameterization, much influenced by the Fermi statistics of electrons.

We very much think that such an approach is worth pursuing, and we have proposed in this regard to use real space partitioning techniques that are able to isolate electron-pair regions and to couple them to real space energy partitioning methods [13,14]. This way, we aim to find how well these techniques are able to recover known chemical insight and to move slowly toward bond-centric energy landscapes that can offer accurate results at a much smaller and less arbitrary parameterization effort.

In a previous work [14] we have shown how the abovementioned coupling recovers and justifies the bond charge model (BCM) proposed by Parr and Borkman [15,16]. To do that we exploited the topology induced by the Electron Localization Function (ELF) proposed by Becke [17], given its ability to define electron pair regions in real space [18], and its good performance in the understanding of chemical properties and structure [19,20], as well as in reactivity [21,22]. We coupled the ELF with the interacting quantum atoms (IQA) approach [23,24], a powerful real-space quantum-chemical topology tool that allows partitioning the total energy of a system into chemically meaningful contributions, and has also been successfully applied in the analysis of a number of chemical phenomena from an intuitive chemical bond perspective [25,26,27]. With this methodology, it was shown that all kinetic, electrostatic, and exchange-correlation terms for the electron pairs follow straightforward scaling laws that were proven to be valid in widely different bonding regimes, ranging from homopolar to highly polar. It was clearly shown that the original BCM model, as well as its updated version, in which it is coupled with the ELF [28], works well in as much the quantum mechanical correction represented by the exchange-correlation term is considerably smaller than all the other contributions entering into the BCM energy. In this sense, it is remarkable to note that this model has recently been applied to the description of a variety of C–C bonds [29].

Kinetic, electrostatic and exchange-correlation terms have been shown to follow distance dependencies that are akin to straightforward scaling laws. Electron pairs, as described by the ELF, can be modelled from the homogeneous electron gas point of view. The proposed scaling behaviors have proven to work reasonably well for different polarities (from homonuclear to highly polar bonds, including non-conventional bonds). Whenever the electron pair fails to describe bonding, like in metallic bonds, or when highly delocalized bonding regimes set in, the model is expected to fail. We have also shown how simple local measures can be introduced to identify these situations [14].

Application of the IQA-ELF methodology (the combination of IQA energy decomposition scheme with the ELF topology in the BCM framework) is still in its infancy, although the ease with which it deals with otherwise difficult-to-model entities, like lone pairs, encourages us to progress in its development. As a further step in such an enterprise, we have decided to invest our efforts in checking to what extent it may be able to recover more subtle electronic effects, particularly those which are usually explained from an atomic-based rather than from a bond-based point of view. Out of the many possible examples out there, we have decided to study the effect of substituents in electrophilic aromatic substitutions (EAS) [30]. These are thoroughly studied reactions, characterized by strong orientation/activation effects triggered by the nature of the substituents [31,32]. As every fresh organic chemistry student learns, both inductive as well as resonance or mesomeric components are usually invoked to rationalize the experimentally observed behavior. Rules to predict the ortho-/para- or meta-orienting ability of substituents are an active topic in chemical research and have been recovered from almost every imaginable point of view [33,34,35,36,37,38,39,40,41]; and even the EAS reaction mechanism is still under debate [30,42,43,44,45,46,47]. Among the various previous approaches for EAS study, we are especially interested in the ELF due to its direct relation with the BCM model and the IQA-ELF coupling. Notice that the ELF application to this kind of processes was first addressed by Silvi and coworkers [48,49,50], who were able to predict the most reactive position in terms of an ELF-based reactivity index in mono-substituted benzene molecules, and then extended it to polyaromatic molecules and to aromatic heterocycles.

Since standard wisdom is atom-centric, here, we want to examine whether the bond-centric IQA-ELF point of view is able to provide further insights into this problem. In this first incursion we will just examine IQA-ELF quantities in isolated substituted benzene rings, much in the light of previous treatments performed within the conceptual density functional theory framework [48,51].

An interesting point regards the absence of σ-π channels in the orbital invariant treatments provided by real space techniques, and, for this distinction, it becomes essential to isolate inductive (I) from mesomeric (M) effects [52]. In this sense, it is relevant to recognize that, in the last decade, a number of works have clearly shown that we can relate these chemical concepts to the spatial decay rate of several electronic reduced density matrices. The exponentially decaying nature of delocalization measures in insulators becomes linked to Kohn’s shortsightedness of the one-particle density and to short-range inductive effects, while the generally oscillating and geometrical decay rate of this quantities in metallic-like systems is connected to long-range mesomerism [53,54]. It will be shown that we can recover the classical insights about I and M effects from the decay of IQA-ELF quantities with the distance to the substituent.

2. Results and Discussion

Table 1 contains a succinct list of X substituents classified according to their activating/deactivating nature as well as to their orienting ability [31]. Its purpose is illustrative, and we will refer to it to contextualize some of our findings. We start by considering briefly the topology of the ELF gradient field, which is shown in Figure 1 for some representative examples. Notice that the results are shown in a graphical manner for the sake of clarity, while the numerical data and some alternative figures are provided in the Supplementary Materials.

2.1. ELF Topology

As originally put forward by Fuster et al. [48], monosubstitution of benzene does not alter the basic ELF topology of the aromatic ring. The phenyl subspace displays the same number and type of basins as benzene, with standard core regions for all carbon atoms, C(C_i) i = 1,6, as well as V(C_i,H) and V(C_i,C_i±1) disynaptic valence domains (note that the list i = 1.6 is circular, so C(C₁) is linked to V(C₁,C₂) and V(C₁,C₆)). In this sense, the ELF field is structurally stable against electrophilic aromatic substitution, easing the comparison of Ph–H (Figure 1a) with Ph–X derivatives (see Figure 1b–f for 5 substituted systems). Substitution does alter the phenyl subspace bifurcation tree, and two kinds of diagrams were found according to the nature of the X substituent that could resolve ortho/para from meta orienting moieties [48]. Although this analysis showed how the perturbative approach that examines the electronic structure of the reactant alone can be correlated to the observed experimental behavior, it did not offer any clear physical rationalization of it. In order to offer such a link, it is advisable to turn to basin observables, be they based on electron populations and their fluctuations or, as we will show, on IQA-ELF quantities.

2.2. Recognizing Induction and Mesomerism

It is interesting to check how inductive and mesomeric effects can be detected in a bond-centric decomposition such as that provided by the ELF field. A number of basin expectation values can be used in this regard, although the simplest one is, no doubt, the basin electron population. Since we are mainly interested in changes in the aromatic ring, we will only consider V(C_i,C_i±1) populations in the following. Let C₁ be the substituted carbon atom, and let us label the V(C_i,C_i+1) = V(C_i,C_j) bond basin and its electron population as bij and nij, respectively, as depicted in Scheme 1. Figure 2 gathers some representative basin populations across the ring, while a full list can be found in Supplementary Table S1. We will base our arguments in the following paragraphs on the results found in this figure, although they are general to all the systems we have studied. Notice that in less symmetrical substituents (like Ph–SH), nearly equivalent basin populations like n23 and n56 were found in some cases to be slightly dissimilar and thus averaged. The same procedure was applied to kinetic and exchange-correlation energies and delocalization indices.

Most of the chemical knowledge regarding the electron density accumulations or depletions that the X substituent induces in the aromatic ring should be found in the nij populations. Both the carbon cores and the V(C,H) basins do surely show variations with X, but these are not directly related to ring electron delocalization. Our data are mainly in line with standard wisdom, although several interesting anomalies stand out.

In the first place, it is rather clear that, a bit unexpectedly, the total number of electrons in the aromatic ring, which can be measured by the average nij population (red spot in Figure 2), is generally larger than that found in the unsubstituted benzene (red dashed line), see Supplementary Figure S1 for nij relative to benzene (Δnij). It is only smaller with strong EWGs, like BH₂ or CN, but not in the case of strong deactivating groups like NO₂ or in weaker ones like the halides. There is indeed a tendency toward larger/smaller total populations in the case of activating/deactivating groups, but a loose one. This warns about the blind use of arrow pushing techniques, which dominate textbook explanations.

A salient feature in Figure 2 is the very different behavior of the population of the bond basins with X. On the one hand, we find systems like toluene, Ph–CF₃, or Ph–NH₃⁺, for which the substituent effect decays with the distance, n12 > n23 > n34. This extinction of the density perturbation induced by X as the distance to the substituent increases agrees with an insulator-like behavior of the one-particle density matrix, which we here associate to (dominant) inductive effects. Notice, however, that the methyl group, usually classified as a weak EDG, leads to a distribution not so far from that found for CF₃, which is however a traditionally strong EWG. It is also clear that large inductive effects lead to stronger perturbations, and that n34 falls well below the reference Ph–H value in Ph–NH₃⁺, for instance. Activating groups tend to display smaller n12 values (in the case of O⁻, n12 = 2.49 e), while deactivating groups usually show larger n12 populations, peaking at 2.92 e for Ph–NH₃⁺.

A very different kind of distance behavior is also obvious. In many of our systems, the nij values oscillate, and n23 is maximum. A comparison of Table 1 and Figure 2 (or Supplementary Table S1 for the full set of molecules) leads to identifying these systems as those where mesomeric effects are important. Interestingly, molecules displaying +M and −M classical resonance give rise to the same kind of oscillatory evolution, with large n23 and smaller n12 and n34 populations. In the +M case, n23 tends to be clearly larger than in the −M one, peaking at n23 = 3.29 e when X = O⁻. For both +/−M systems a clear tendency to a n12 < n34 configuration is also evident. This behavior is relatively consistent with the naïve resonance structures that include a C₁ = X double bond, thus reinforcing the quinoid-like C₂–C₃ links in most resonance forms, at the expense of the C₁–C₂ and C₃–C₄ ones (see Scheme 2).

Halogens display a peculiar behavior. On the one hand, their ring bond populations decay with the distance, pointing to dominating inductive effects. On the other, their average number of electrons is quite large, in line with other activating systems in the upper part of Table 1. Fluorine, in particular, is close to showing a mesomeric-like maximum at n23 that points toward its standard classification as a π-donor.

The nitro group also deserves some specific comments. Despite being usually classified as a strong −I/−M deactivating group, in the same category as CN, our data does not find any (dominant) mesomeric effect, but a behavior much more in line with other −I deactivating groups such as CF₃ or NH₃⁺. Cyanide, on the contrary, does clearly show an oscillating M pattern. Given that the −M classification of NO₂ is based on the purported positive π hole of the N atom in the NO₂ Lewis structure, something which is not supported by real space chemical bonding analyses, we think that the present results might be better representing the true nature of the nitro group deactivating features. When the π-Lewis acid nature of a group is out of question, as in the BH₂ case, the oscillations are clearly present. Finally, Ph–Li is also interesting on its own. The Li atom has a very low electronegativity and no π density channels, basically leaving a Ph⁻ anion characterized by a very small aromatic population that responds mesomerically to the perturbation. Actually, its average number of aromatic bond population is the second lowest after BH₂.

As already commented, these observations are compatible with knowledge about how density perturbations decay with the distance in molecular systems. As proven by Walter Kohn in a seminal paper [55], insulators are characterized by an exponential drop of electron-electron correlations with the distance, which ends up in the non-diagonal elements of the one-particle density matrix (1DM) also decaying exponentially,ρr;r+s≈^e−ηs. On the contrary, the 1DM falls algebraically in metals,ρr;r+s≈a^s−d , where the exponent d depends on the dimensionality of the system. This difference supports the textbook chemical picture of localized/insulating, delocalized/metallic electrons. We have also shown [53,54] that the algebraic decay is found in unsaturated alternant hydrocarbons, where an overall inverse quadratic envelope is enriched by an oscillatory pattern of the real space delocalization indices (DI), which measure the fluctuation of the electron populations in two different spatial regions, or the number of shared electrons between them. It was found that when atomic regions are used, the DI(A,B) values with A and B separated by an even number of bonds (grossly corresponding to meta pairs) vanish in the tight-binding approximation, and that they turn out to be much smaller than DIs between pairs separated by an odd number of bonds in actual calculations.

2.3. Bond Delocalization and Mesomerism

The findings in the previous subsection have prompted us to examine the ability of the different aromatic valence basins to delocalize their electrons among themselves as a way to rationalize the substituent effects in electrophilic attacks. To do so, we examine vicinal DIs, i.e., delocalization indices between adjacent V(C_i,C_i+1) and V(C_i+1,C_i+2). These DIs can be identified by the common C_i+1 atom, which may correspond to either the substituted position, ipso, or to the ortho, meta, or para ones, as shown in Scheme 3. We will thus refer to them as the DI(a) descriptors, with a = i,o,m,p. A summary of our data for the same set of selected systems as for Figure 2 is gathered in Figure 3. Values are collected relative to the unsubstituted benzene molecule in order to ease the analysis in terms of substitution induced enhancements. The full set of results can be found in Supplementary Table S2.

Before performing a detailed analysis of the data, it is relevant to point out that the DI(a) descriptors do have a clear physical meaning. They exactly measure how correlated the populations of vicinal aromatic bonds are, since DI(A,B) = −2·cov(n_A,n_B), with A,B being two spatial domains, n_A,n_B their respective electron populations, and cov(n_A,n_B) the covariance of these populations. In other words, the DI(a) descriptor tells us about how easily electrons traverse the C_A atom, flowing through it. Since in the electrophilic attack of the C_A moiety the arenium ion strongly modifies its vicinal links, we expect ΔDI(a) = DI(a,X) − DI(a,X = H) to correlate with some of the experimental results.

A quick look at Figure 3 shows that this is indeed the case. First, we mention that these descriptors are not completely independent from those in Figure 2. For instance, ΔDI(i) is proportional to n12, and whenever the latter population is larger than the one in benzene, ΔDI(i) is positive and vice versa. This is to be expected, since the number of shared pairs of electrons depends on the basin populations. Again, the importance of mesomeric effects for a substituent can be grasped from the data through a negative ΔDI(i). When inductive effects dominate, ΔDI(i) is positive, and this again includes the nitro group.

Let us now focus on the changes in ΔDI(i,o,m,p). First, classical activating groups such as NH₂, SH, or CH₃ display a +,−,+ pattern for the o,m,p descriptors, respectively (see Table 1). This corresponds to a greater delocalizing ability than plain benzene at the ortho and para sites, while a smaller one at the meta position. Notice that the first two groups are classified as −I/+M, while the methyl substituent is usually understood as a pure +I substituent. Although the magnitude of ΔDI(i) is related to the strength of the substituent effect, care has to be taken due to the influence of the underlying basin populations. For instance, fluorine is found to be ortho/para directing quite as much as NH_2, but its larger basin populations (Figure 2) mask its strength. A suitably normalized descriptor might solve this issue, but has not been pursued in this contribution. We do not skip the fact that the oxido group displays a large negative ΔDI(p) value. We do not have a fully consistent explanation for this fact, but as in the case of the nitro group behavior, we tend to think that the traditional classification may hide other interesting effects that should be further investigated.

Similarly, deactivating groups display negative ΔDI(o,p) values, with meta contributions either positive or negative. The change from a +,−,+ to a −,+,− pattern as we change X is very interesting, and demonstrates that bond-centric descriptors do also pick up a coherent image of the electronic effects under study. X = CF₃ (and also CHF₂, or CH₂Cl, see Supplementary Table S2) displays triply negative values, BH₂ a considerable enhancement of the meta index, and the NH₃⁺ group a small activation of the ortho position, which might be masked in a true electrophilic attack due to electrostatic repulsion with the attacking moiety.

The experimental substituent effects have been usually rationalized by means of the Hammett equation, in which ratios of reaction rates are interpreted in terms of σ and ρ constants which consider substituent and reaction conditions effects, respectively. Figure 4 shows how σ_m − σ_p correlates reasonably with DI(m) − DI(p) for a number of systems. The Hammett constants have been taken from Ref. [56]. Similar correlations were found through local ELF bifurcation values [48], demonstrating that our physically based descriptors are able to capture, at least grossly, the subtleties of standard organic chemistry knowledge.

2.4. IQA-ELF Descriptors

We have shown how basin electron populations and fluctuations can be used to offer a consistent picture of inductive and resonance effects without an explicit σ − π separation. We now turn to the less studied basin energetic expectation values. We will start with the kinetic energy. As in previous works, we use the positive definite kinetic energy density, see Equation (1).

t=12∇·^∇′ρ(r;^r′)

For only in the case of QTAIM basins we can use any Laplacian based kinetic energy density. This is a common practice, and avoids possible negative basin expectation values.

Kinetic energies of the aromatic valence basins are collected in Figure 5 (and Supplementary Table S3 for the complete set of systems), with a notation equivalent to that of Figure 2 (provided in Scheme 1), and the kinetic energy with relative to benzene (Δtij) is shown in Supplementary Figure S2. As expected, the basin kinetic energies are proportional to the basin populations, so that as a first approximation, Figure 2 and Figure 5 show similar trends.

More interesting is the consideration of the domain kinetic energy per electron, which can be constructed immediately from the data. This is 1.166 au in benzene, and deviations from this result are quite small in all cases. A detailed analysis shows that t23/n23 and t34/n34 are almost unaffected by substitution, and that the largest variations are found in t12/n12 which is larger than in benzene. This indicates that the perturbation of the nodal structure of the wavefunction, or of the curvature of the 1DM, falls quite fast with the distance to the substituent. For instance,Δ_t12n12is equal to 0.033, 0.028, and 0.027 au for the oxido, amino, and nitro X moieties, respectively.

Finally, the only remaining IQA-ELF energetic parameter that we expect to be related to orientation effects in electrophilic aromatic substitutions is the inter-basin exchange-correlation energy,_VxcAB . Electron-nuclear attractions are density dependent quantities that are not expected to respond clearly to X effects, so they will be skipped in this first contribution. As repeatedly shown [25], V_xc is a measure of the covalent energy associated to the interaction between the electrons contained in domains A and B. For the sake of clarity and given its direct connection with the delocalization index, we represent V_xc for the pairs of basins under study with respect to benzene (ΔV_xc) basins in Figure 6 (see Supplementary Table S4 for the numeric values of the complete set of molecules). We again label the interaction by the common carbon atom in the adjacent bond basins, as indicated in Scheme 3 for delocalization indices.

As it can be immediately checked, exchange-correlation energies mimic population fluctuations [57] and thus inter-basin electron delocalization. It is particularly interesting that the energy descriptors we have investigated in this contribution follow rather closely their population analogs. On the one hand, kinetic energies are found to evolve as populations. On the other, exchange-correlation energies can be sensed by electron delocalization indices. This is of considerable importance since the computation of energetic basin observables, particularly in the case of two-electron ones, is extremely computationally intensive. In the pseudo one-determinant approximation that we are using here, for instance, the calculation of V_xc requires 6D numerical grids, that can be computed with O(N⁴) effort if bipolar expansions are used [58]. Meanwhile, the calculation of delocalization indices is a simple byproduct of the construction of the domain overlap matrices for the set of occupied orbitals, just needing 3D grids.

Overall, we think that building substituent-responsive force fields based on electron-pair energetic descriptors is within reach, for the IQA energetic domain expectation values have been shown to capture the essentials of the well-known, nevertheless subtle inductive and mesomeric effects in the simplest electrophilic aromatic substitutions. 3. Theoretical and Computational Details

All electronic structure calculations and wavefunctions for ELF and IQA analyses were obtained through DFT calculations, via the Gaussian 09 suite (Gaussian, Inc.; Wallingford, CT, USA) [59]. The B3LYP exchange-correlation functional [60,61], in conjunction with the 6-31G(d) basis set [62,63] were used. The method selection was validated in Ref. [14], with respect to the basis set and calculation method.

In order to obtain the domains associated to bonding and non-bonding pairs, we used the Electron Localization Function (ELF) [17] in its simplest formulation. The ELF kernel,_χσ , which was originally defined in terms of the curvature of the Fermi hole [64], can be interpreted as a measure of the excess of local kinetic energy due to the Pauli Principle, relative to a homogeneous electron gas kinetic energy density [65]. The ELF, η(r), is then obtained by mapping its kernel,_χσ, through a Lorentzian (see Equation (2)) in order to obtain a function that ranges from 0 to 1.

η(r→)=11+_χσ2

This scalar field is used to induce an exhaustive topological partition of space into non-overlapping domains (basins) whose properties can be determined by integrating appropriate operator densities over their associated volumes. When the electron density, ρ(r), is integrated, for instance, the average number of electrons in the domain is recovered. ELF basins have been calculated through the TopMod program [66], with a tridimensional grid of 150 points in each direction. The energetic descriptors of these topological basins have been calculated within the IQA energy decomposition scheme [23]. IQA is usually applied to the QTAIM atomic decompositions, leading to Equation (3).

E=∑A(_TA+_VeeAA+_VenAA)︸_{∑A EintraA}+∑A>B(_VenAB+_VneAB+_VnnAB+_VeeAB)︸_{∑A>B EinterAB}

In this expression we sum over atomic regions, and the total energy is decomposed into intra-atomic contributions,_EintraA, and inter-atomic interactions between each pair of atoms,_EinterAB. The former is the sum of the kinetic energy of electrons in atom A,_TA, and the electron-electron and electron-nucleus interactions,_VeeAAand_VenAA, respectively._EinterABis obtained by adding the interaction among the electrons in A and the nucleus of B,_VenAB; its BA symmetric counterpart,_VneAB; and the nucleus-nucleus and electron-electron interactions between A and B,_VnnABand_VeeAB, respectively. The electron repulsion terms are usually divided in IQA into a classical contribution that depends only on the electron density,_VelecAB, and a non-classical or exchange-correlation one,_VXCAB, as shown in Equation (4).

_VeeAB=_VelecAB+_VXCAB

Let us now isolate the_VXCABterm and include_VenAB,_VneAB,_VnnAB, and the classical contribution of_VeeAB(_VelecAB) in the classical Coulombic energy, which is thus defines as indicated in Equation (5):

_VCoulAB=_VenAB+_VneAB+_VnnAB+_VelecAB

The classical Coulomb contribution for the intra-atomic terms is analogous, see Equation (6).

_VCoulAA=_VenAA+_VneAA+_VelecAA

As a result, intra- and inter-atomic energy contributions can be expressed according to Equations (7) and (8).

_EintraAA=_TA+_VCoulAA+_VXCAA

_EinterAB=_VCoulAB+_VXCAB

Since IQA is based on a given exhaustive partition of space, it can also be applied to ELF partitions. We have already shown how this can be done [13]. Notice that this allows us to include some ELF basins related to electron pairs as lone pairs (monosyaptic basins) and bonds (disynaptic and higher order) in addition to atomic-like contributions (the core basins). In this regard, it is necessary to highlight that the former lack nuclei in their interior, so that there are no intra-domain electron-nuclear potential terms for them.

Since we are interested in interactions between C–C bonds, let us take a general V(C_i,C_j) disynaptic basin, and label it as “Bond”. Given that this basin does not involve any nucleus:_VenAA=_VenBond=0, and_EintraBondwill read as shown in Equation (9).

_EintraBond=^TBond+_VCoulBond+_VXCBond

For the “Bond” basin interaction with the core basin of atom A, for instance, it turns out that_VenA-Bondand_VnnA-Bondvanish, and the A-Bond interaction can be calculated according to Equation (10).

_EinterA-Bond=_VneA-Bond+_VCoulA-Bond+_VXCA-Bond

Finally, the interaction between two basins representing two different bonds, say “Bond1” and “Bond2” is given by Equation (11). Notice that, in this case, since any of the interaction entities bear any nucleus:_{VenBond1-Bond2}=_{VnnBond1-Bond2}=_{VenBond1-Bond2}=_VnnA-Bond=0, that is to say, the only non-vanishing term is_{VeeBond1-Bond2}.

_{EinterBond1-Bond2}=_{VCoulBond1-Bond2}+_{VXCBond1-Bond2}

From the previous discussion, it is clear that the use of IQA-ELF leads to an energy decomposition in which the basic interacting units are electron pairs, thus leading to interactions between cores and bonds, or between the possible combinations of bonds and lone pairs. IQA-ELF calculations have been performed with a modified version of the PROMOLDEN code that integrates over ELF instead of over QTAIM basins [13]. ELF representations were performed by means of UCSF Chimera package [67,68].

A set of 28 mono-substituted benzene derivatives (Ph–X) have been optimized at the above-mentioned level of theory. These include inductively (−I) and resonance mediated (−M) electron withdrawing groups (EWG), as well as inductively (+I) and resonance mediated (+M) electron donating groups (EDG). The list includes, in addition to benzene, the Me, Et and Pr alkyl groups, the F, Cl, and Br halides, and the CHF₂, CH₂Cl, and CF₃ species, the NO₂, O^−, NH₂, NH₃⁺, SH, SiH₃, SiH₂Me, CN, CCH, CC–F, CC–Cl, CC–Me, and Ph groups, and also a number of non-standard substituents like Li, BH₂, BMe₂, BF₂, and AlH₂. Common wisdom establishes that mesomeric (M) effects are generally stronger than inductive ones, and that EWGs deactivate the phenyl ring toward electrophilic attack by decreasing its electron density while EDGs activate it. Usually, activation and deactivation is stronger at the ortho/para positions (especially when it is due to mesomeric effects), meaning that, in general, deactivating groups direct electrophiles to the meta position [31]. Note that we did not include OH and OMe substituents in our survey due to integration issues generated by the irregular shape of the V(C,O) basin, which led to rather large errors with the PROMOLDEN code for these two substituents.

4. Conclusions In this contribution, we have shown that inductive and mesomeric effects can be grasped from the behavior of the population of the aromatic valence bond basins. Dominating inductive effects decay exponentially with the distance to the perturbing center, as in insulators, while resonance or mesomeric effects show an oscillatory pattern that recalls known behavior in metallic systems. Delocalization indices between adjacent C–C bonds provided by the IQA-ELF framework have been found as particularly suited expectation values to that end. In agreement with chemical intuition, when a substituent increases the electron flow ability between the two aromatic bonds that share a given carbon atom of the ring, the position tends to be activated toward electrophilic attack and vice versa. This good agreement strongly supports the good performance of IQA-ELF methodology in capturing the chemical behavior and making accurate predictions. Moreover, we have convincingly shown that the basin kinetic energy evolves in consonance with basin populations, and that exchange-correlation energies provide a similar energetic analog to the population fluctuations. The relevance of this work is double, as, on the one hand, we have provided a new energetic descriptor for the prediction of EAS reactions based on the ELF topology. On the other hand, we have proven that the IQA-ELF scheme works well in predicting chemical reactivity. This will be very useful for paving the way to the construction of well-behaved electron-pair energetic descriptors to be used in future force field developments, as they accurately capture the chemical behavior of the systems under study.

Supplementary Materials

The following are available online, Table S1. ELF basin populations, Table S2. Delocalization indices relative to benzene, Table S3. Kinetic energy, Table S4. Exchange-correlation energies relative to benzene, Table S5. Cartesian coordinates, Figure S1. ELF basin populations relative to benzene, Figure S2. Kinetic energy with relative to benzene.

Author Contributions

Conceptualization, Á.M.P., J.M., and J.C.-G.; methodology, J.M., J.C.-G., and Á.M.P.; software, Á.M.P.; validation, formal analysis, data curation, J.M. and M.G.; writing, editing, visualization, J.M. and Á.M.P.; funding acquisition, Á.M.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Spanish MICINN, grant number PGC2018-095953-B-I00, and by the FICyt, grant IDI/2018/000177. MG thankfully acknowledges the Spanish "Ministerio de Ciencia, Innovación y Universidades (MICIU)" for a predoctoral FPU fellowship (FPU19/02903).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Sample Availability

Samples of compounds are not available from the authors.

References

1. Houk, K.N.; Liu, F. Holy Grails for Computational Organic Chemistry and Biochemistry. Acc. Chem. Res. 2017, 50, 539-543.

2. Willems, H.M.G.; De Cesco, S.; Svensson, F. Computational Chemistry on a Budget: Supporting Drug Discovery with Limited Resources. J. Med. Chem. 2020, 63, 10158-10169.

3. Huang, S.-D.; Shang, C.; Zhang, X.-J.; Liu, Z.-P. Material discovery by combining stochastic surface walking global optimi-zation with a neural network. Chem. Sci. 2017, 8, 6327-6337.

4. Butler, K.T.; Frost, J.M.; Skelton, J.M.; Svane, K.L.; Walsh, A. Computational materials design of crystalline solids. Chem. Soc. Rev. 2016, 45, 6138-6146.

5. Curtarolo, S.; Hart, G.L.; Nardelli, M.B.; Mingo, N.; Sanvito, S.; Levy, O. The high-throughput highway to computational materials design. Nat. Mater. 2013, 12, 191-201.

6. Gooneie, A.; Schuschnigg, S.; Holzer, C. A Review of Multiscale Computational Methods in Polymeric Materials. Polymers 2017, 9, 16.

7. Mater, A.C.; Coote, M.L. Deep Learning in Chemistry. J. Chem. Inf. Model. 2019, 59, 2545-2559.

8. Dral, P.O. Quantum Chemistry in the Age of Machine Learning. J. Phys. Chem. Lett. 2020, 11, 2336-2347.

9. Senior, A.W.; Evans, R.; Jumper, J.; Kirkpatrick, J.; Sifre, L.; Green, T.; Qin, C.; Žídek, A.; Nelson, A.W.R.; Bridgland, A.; et al. Improved protein structure prediction using potentials from deep learning. Nature 2020, 577, 706-710.

10. Leach, A. Molecular Modelling: Principles and Applications, 2nd ed.; Prentice Hall: Upper New Jersey River, NJ, USA, 2001.

11. Brázdová, V.; Bowler, D.R. Atomistic Computer Simulations: A Practical Guide; Wiley-VCH: Weinheim, Germany, 2013.

12. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157-1174.

13. Martín Pendás, Á.; Francisco, E.; Blanco, M. Electron-electron interactions between ELF basins. Chem. Phys. Lett. 2008, 454, 396-403.

14. Munarriz, J.; Laplaza, R.; Martín Pendás, Á.; Contreras-García, J. A first step towards quantum energy potentials of electron pairs. Phys. Chem. Chem. Phys. 2019, 21, 4215-4223.

15. Parr, R.G.; Borkman, R.F. Simple Bond-Charge Model for Potential-Energy Curves of Homonuclear Diatomic Molecules. J. Chem. Phys. 1968, 49, 1055-1058.

16. Borkman, R.F.; Parr, R.G. Toward an Understanding of Potential-Energy Functions for Diatomic Molecules. J. Chem. Phys. 1968, 48, 1116-1126.

17. Becke, A.D.; E Edgecombe, K. A simple measure of electron localization in atomic and molecular systems. J. Chem. Phys. 1990, 92, 5397-5403.

18. Silvi, B.; Savin, A. Classification of chemical bonds based on topological analysis of electron localization functions. Nat. Cell Biol. 1994, 371, 683-686.

19. Munarriz, J.; Calatayud, M.; Contreras-García, J. Valence-Shell Electron-Pair Repulsion Theory Revisited: An Explanation for Core Polarization. Chem. A Eur. J. 2019, 25, 10938-10945.

20. Contreras-García, J.; Mori-Sánchez, P.; Silvi, B.; Recio, J.M. A Quantum Chemical Interpretation of Compressibility in Solids. J. Chem. Theory Comput. 2009, 5, 2108-2114.

21. Andrés, J.; González-Navarrete, P.; Safont, V.S.; Silvi, B. Curly arrows, electron flow, and reaction mechanisms from the perspective of the bonding evolution theory. Phys. Chem. Chem. Phys. 2017, 19, 29031-29046.

22. Munárriz, J.; Laplaza, R.; Polo, V. A bonding evolution theory study on the catalytic Noyori hydrogenation reaction. Mol. Phys. 2019, 117, 1315-1324.

23. Blanco, M.A.; Martín Pendás, A.; Francisco, E. Interacting Quantum Atoms: A Correlated Energy Decomposition Scheme Based on the Quantum Theory of Atoms in Molecules. J. Chem. Theory Comput. 2005, 1, 1096-1109.

24. Maxwell, P.; Martín Pendás, Á.; Popelier, P.L.A. Extension of the interacting quantum atoms (IQA) approach to B3LYP level density functional theory (DFT). Phys. Chem. Chem. Phys. 2016, 18, 20986-21000.

25. Guevara-Vela, J.M.; Francisco, E.; Rocha-Rinza, E.; Martín Pendás, A. Interacting Quantum Atoms-A Review. Molecules 2020, 25, 4028.

26. Munarriz, J.; Velez, E.; Casado, M.A.; Polo, V. Understanding the reaction mechanism of the oxidative addition of ammonia by (PXP)Ir(i) complexes: The role of the X group. Phys. Chem. Chem. Phys. 2018, 20, 1105-1113.

27. Fernández-Alarcón, A.; Guevara-Vela, J.M.; Casals-Sainz, J.L.; Costales, A.; Francisco, E.; Martín Pendás, Á.; Rinza, T.R. Photochemistry in Real Space: Batho- and Hypsochromism in the Water Dimer. Chem. A Eur. J. 2020, 26, 17035-17045.

28. Contreras-García, J.; Marqués, M.; Menéndez, J.M.; Recio, J.M. From ELF to Compressibility in Solids. Int. J. Mol. Sci. 2015, 16, 8151-8167.

29. Laplaza, R.; Polo, V.; Contreras-García, J. A Bond Charge Model Ansatz for Intrinsic Bond Energies: Application to C-C Bonds. J. Phys. Chem. A 2020, 124, 176-184.

30. Galabov, B.; Nalbantova, D.; Schleyer, P.V.R.; Schaefer, I.H.F. Electrophilic Aromatic Substitution: New Insights into an Old Class of Reactions. Accounts Chem. Res. 2016, 49, 1191-1199.

31. Vollhardt, K.P.C.; Schore, N.E. Organic Chemistry: Structure and Function; W.H. Freeman, Macmillan Learning: New York, NY, USA, 2018.

32. Taylor, R. Electrophilic Aromatic Substitution; Wiley: New York, NY, USA, 1990.

33. Brown, J.J.; Cockroft, S.L. Aromatic reactivity revealed: Beyond resonance theory and frontier orbitals. Chem. Sci. 2013, 4, 1772-1780.

34. Zhou, Z.; Parr, R.G. Activation hardness: New index for describing the orientation of electrophilic aromatic substitution. J. Am. Chem. Soc. 1990, 112, 5720-5724.

35. Koleva, G.; Galabov, B.; Wu, J.I.; Schaefer, H.F., III; Schleyer, P.v.R. Electrophile Affinity: A Reactivity Measure for Aromatic Substitution. J. Am. Chem. Soc. 2009, 131, 14722-14727.

36. Schnatter, W.F.K.; Rogers, D.W.; Zavitsas, A.A. Electrophilic Aromatic Substitution: Enthalpies of Hydrogenation of the Ring Determine Reactivities of C6H5X. The Direction of the C6H5-X Bond Dipole Determines Orientation of the Substitution. J. Phys. Chem. A 2013, 117, 13079-13088.

37. Franco-Pérez, M.; Polanco-Ramírez, C.A.; Gázquez, J.L.; Ayers, P.W.; Vela, A. Study of organic reactions using chemical re-activity descriptors derived through a temperature-dependent approach. Theor. Chem. Acc. 2020, 139, 44.

38. Hadzic, M.; Braïda, B.; Volatron, F. Wheland Intermediates: An ab Initio Valence Bond Study. Org. Lett. 2011, 13, 1960-1963.

39. Kromann, J.C.; Jensen, J.H.; Kruszyk, M.; Jessing, M.; Jørgensen, M. Fast and accurate prediction of the regioselectivity of electrophilic aromatic substitution reactions. Chem. Sci. 2017, 9, 660-665.

40. Fievez, T.; Pinter, B.; Geerlings, P.; Bickelhaupt, F.M.; De Proft, F. Regioselectivity in Electrophilic Aromatic Substitution: Insights from Interaction Energy Decomposition Potentials. Eur. J. Org. Chem. 2011, 2011, 2958-2968.

41. Cheshmedzhieva, D.; Ilieva, S.; Galabov, B. Computational evaluation of σI and σR substituent constants. J. Mol. Struct. 2010, 976, 427-430.

42. Stuyver, T.; Danovich, D.; De Proft, F.; Shaik, S. Electrophilic Aromatic Substitution Reactions: Mechanistic Landscape, Electrostatic and Electric-Field Control of Reaction Rates, and Mechanistic Crossovers. J. Am. Chem. Soc. 2019, 141, 9719-9730.

43. Hubig, S.M.; Kochi, J.K. Structure and Dynamics of Reactive Intermediates in Reaction Mechanisms. σ- and π-Complexes in Electrophilic Aromatic Substitutions. J. Org. Chem. 2000, 65, 6807-6818.

44. Koleva, G.; Galabov, B.; Kong, J.; Schaefer, H.F.; Schleyer, P.V.R. Electrophilic Aromatic Sulfonation with SO3: Concerted or Classic SEAr Mechanism? J. Am. Chem. Soc. 2011, 133, 19094-19101.

45. Galabov, B.; Koleva, G.; Simova, S.; Hadjieva, B.; Schaefer, H.F., III; Schleyer, P.v.R. Arenium Ions are not Obligatory. In-termediates in Electrophilic Aromatic Substitution. Proc. Natl. Acad. Sci. USA 2014, 111, 10067-10072.

46. Gwaltney, S.R.; Rosokha, S.V.; Head-Gordon, M.; Kochi, J.K. Charge-Transfer Mechanism for Electrophilic Aromatic Nitration and Nitrosation via the Convergence of (ab Initio) Molecular-Orbital and Marcus-Hush Theories with Experiments. J. Am. Chem. Soc. 2003, 125, 3273-3283.

47. Shernyukov, A.; Genaev, A.M.; Salnikov, G.E.; Rzepa, H.S.; Shubin, V.G. Noncatalytic bromination of benzene: A combined computational and experimental study. J. Comput. Chem. 2016, 37, 210-225.

48. Fuster, F.; Sevin, A.A.; Silvi, B. Topological Analysis of the Electron Localization Function (ELF) Applied to the Electrophilic Aromatic Substitution. J. Phys. Chem. A 2000, 104, 852-858.

49. Silvi, B.; Kryachko, E.S.; Tishchenko, O.; Fuster, F.; Nguyen, M.T. Key properties of monohalogen substituted phenols: In-terpretation in terms of the electron localization function. Mol. Phys. 2002, 100, 1659-1675.

50. Fuster, F.; Sevin, A.; Silvi, B. Determination of substitutional sites in heterocycles from the topological analysis of the electron localization function (ELF). J. Comput. Chem. 2000, 21, 509-514.

51. Poater, J.; Duran, M.; Solà, M.; Silvi, B. Theoretical Evaluation of Electron Delocalization in Aromatic Molecules by Means of Atoms in Molecules (AIM) and Electron Localization Function (ELF) Topological Approaches. Chem. Rev. 2005, 105, 3911-3947.

52. Hammett, L.P. Some Relations between Reaction Rates and Equilibrium Constants. Chem. Rev. 1935, 17, 125-136.

53. Gallo-Bueno, A.; Francisco, E.; Martín Pendás, Á. Decay rate of real space delocalization measures: A comparison between analytical and test systems. Phys. Chem. Chem. Phys. 2016, 18, 11772-11780.

54. Gallo-Bueno, A.; Kohout, M.; Martín Pendás, A. Decay Rate of Correlated Real-Space Delocalization Measures: Insights into Chemical Bonding and Mott Transitions from Hydrogen Chains. J. Chem. Theory Comput. 2016, 12, 3053-3062.

55. Kohn, W. Theory of the Insulating State. Phys. Rev. 1964, 133, A171-A181.

56. Hansch, C.; Leo, A.; Taft, R.W. A survey of Hammett substituent constants and resonance and field parameters. Chem. Rev. 1991, 91, 165-195.

57. Francisco, E.; Crespo, D.M.; Costales, A.; Martín Pendás, Á. A multipolar approach to the interatomic covalent interaction energy. J. Comput. Chem. 2017, 38, 816-829.

58. Martín Pendás, A.; Blanco, M.A.; Francisco, E. Two-electron integrations in the quantum theory of atoms in molecules. J. Chem. Phys. 2004, 120, 4581-4592.

59. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09; revision D.01; Gaussian, Inc.: Wallingford, CT, USA, 2009.

60. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785-789.

61. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652.

62. Dunning, T.H. Gaussian Basis Functions for Use in Molecular Calculations. III. Contraction of (10s6p) Atomic Basis Sets for the First-Row Atoms. J. Chem. Phys. 1971, 55, 716-723.

63. Dunning, T.H.; Peterson, K.A.; Wilson, A.K. Gaussian basis sets for use in correlated molecular calculations. X. The atoms aluminum through argon revisited. J. Chem. Phys. 2001, 114, 9244-9253.

64. Dobson, J.F. Interpretation of the Fermi hole curvature. J. Chem. Phys. 1991, 94, 4328-4333.

65. Savin, A.; Jepsen, O.; Flad, J.; Andersen, O.K.; Preuss, H.; Von Schnering, H.G. Electron Localization in Solid-State Structures of the Elements: The Diamond Structure. Angew. Chem. Int. Ed. 1992, 31, 187-188.

66. Noury, S.; Krokidis, X.; Fuster, F.; Silvi, B. Computational tools for the electron localization function topological analysis. Comput. Chem. 1999, 23, 597-604.

67. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera-A visual-ization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605-1612.

68. Goddard, T.D.; Huang, C.C.; Ferrin, T.E. Visualizing density maps with UCSF Chimera. J. Struct. Biol. 2007, 157, 281-287.