domain: chim info

text:
Johanna Klein

Paul Fleurat-Lessard

Julien Pilmé
email: julien.pilme@sorbonne-universite.fr

New Insights In Chemical Reactivity From Quantum Chemical Topology

Keywords: Quantum Chemical Topology, electron localization function, Chemical Reactivity, electrophile, nucleophile, Frontier Molecular Orbital Theory

I. Introduction.

 it seems more computationally efficient to mimic the chemist and to use reactivity descriptors to guide the search of possible products.

In the framework of molecular orbitals, this was considered by Fukui who introduced the Frontier Molecular Orbital (FMO) theory, [12] and others. [13][14][15] Overall, Fukui functions are powerful tools to get an insight into the reactivity since they allow selecting the electrophilic and nucleophilic regions of the molecular space. Later on, Parr pioneered the conceptual DFT approach (cDFT) and proposed scales to estimate the electrophilic and nucleophilic sites of a system, as well as to quantify the reactivity of these sites. [16] More recently, the interaction energy between a nucleophile and an electrophile has been expressed as a function of reactivity descriptors involving the chemical potential and Fukui functions. [17,18] However, as recalled recently by Geerling et al., most cDFT studies have focused on interpreting the experimental or computational results of a reaction, and less on the prediction of unknown products. [19] Indeed, to our knowledge, only a handful of studies have used these indices semiquantitatively to predict the most stable products of a reaction. [20,21] Meanwhile, alternative strategies of bonding analysis belonging to the Quantum Chemical Topology (QCT) methodology [22][23][START_REF] Silvi | Comprehensive Inorganic Chemistry II[END_REF][START_REF] Popelier | Intermolecular Forces and Clusters I[END_REF][START_REF] Silvi | [END_REF] have been developed for a long time in the field of non-relativistic quantum calculations. These approaches grew up from the pioneering works of Richard F. W. Bader and coworkers in the 70s which was initially motived by the generalization of the quantum mechanical principle of stationary action to a molecular subsystem. [27,[START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF] The purpose of this methodology is to address questions about the chemical bonding in molecules or solids, and describe or explain the chemical reactivity trends. [START_REF] Popelier | Applications of Topological Methods in Molecular Chemistry[END_REF][START_REF] Popelier | [END_REF][31][32][START_REF] Contreras-García | Modern Charge-Density Analysis[END_REF][START_REF] Silvi | [END_REF] These methods have progressively introduced numerous wellestablished one-density descriptors for both covalent [START_REF] Silvi | [END_REF]35] and non-covalent intermolecular interactions in the non-relativistic context [36][37][38][39] as well as for the relativistic wavefunctions. [40][41][42] However, the development of its chemical reactivity component still remains a fascinating challenge. For this latter purpose, some attempts to rationalize the chemical reactivity have been proposed such as the analysis of topological descriptors behaviors along IRC, [43][44][45][46][47][48][49][50][51][52][53] the search of reactivity sites, [54,55] the analysis of the topology of Fukui functions, [START_REF] Fuentealba | Applications of Topological Methods in Molecular Chemistry[END_REF][START_REF] Tiznado | [END_REF] some recent development based on the distributed electrostatic moments [58,59] or other developments inspired by the set of VSEPR rules. [60,61] In this same spirit, the hardy development of energy decomposition analysis (EDA) methods [62,63] such as the Interacting Quantum Atoms (IQA) [23,64,65] energy scheme applied to the QTAIM partition or more recently to the electron localization function (ELF) partition [66,67] provided an alternative to extend the FMO and conceptual DFT approaches, since the definitions of atomic properties remain still valid under this context.

The present work addresses an original strategy designed to combine on equal footing treatment between QCT and cDFT schemes using the termed ELF x function recently proposed to identify both electrophilic and nucleophilic domains of the molecular space. We aim at showing that this combined approach can be used quantitatively to predict not only the most stable outcome of a reaction, but also to give the shape of the energy profile from the reactants to the products.

II. Quantum Chemical Topology.

In this article, we assume that the reader is familiar with the quantum chemical topology (QCT) of scalar fields, numerous presentations and applications being previously published in the literature. [22,[START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF][START_REF] Popelier | Applications of Topological Methods in Molecular Chemistry[END_REF] Briefly, a QCT gives a partition of the molecular space into subsystems (the so-called basins) achieved by applying the theory of dynamical gradient systems to the properties of the one-electron density function. The topology of the gradient field is characterized by its critical points (CP), where the gradient of the electron density is zero, and by their connectivity. CP can be either maxima (3, -3), minima (3, 3) or saddle points. The basins are localized around the (3, -3) CP of the function and are separated by the zero flux surfaces. In the QTAIM framework, the function considered is the electron density and these basins are associated with each atom in the molecule. [27,[START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF] Note that non-nuclear attractors (NNAs) can be found to very flat electron density maxima occurring exceptionally in the internuclear regions of metal clusters. [68] Other scalar fields can be used such as ELF [35] which is typically interpreted as a signature of the electronic-pair distribution. [START_REF] Silvi | [END_REF] The functional form of ELF [69,70] is based on the total electron density (closed-shell or open-shell systems), its gradient and the kinetic energy density τ. The kernel of ELF is then defined as: where τ is the positive definite kinetic energy density and (r) is the total electron density.

Electron localization from density components. One of us has recently introduced the termed ELF x function in order to identify the electrophilic and nucleophilic domains which are spread out over the molecular space [54] . The function was defined from ELF as follows:

Where x(r) is a normalized dimensionless quantity that can be expressed from the field of the frontier molecular orbitals, [12,71] τ c and

(1)



is the total electron density of a molecular system with N electrons.  is the total electron density of the molecular system with N+1 electrons with the same geometry and the same orbitals that are obtained for the system with N electrons. When the HOMO or LUMO orbitals are (quasi)degenerated, x(r) can be augmented as:

  
The kernel x is a meaningful measure of electron localization for only a subpart of the total density system. The localization is inversely proportional to the number of electrons delocalized in the subsystem, i.e. if x(r) 0, the function x (r) goes to infinity and then 1/(1+ x (r) 2 ) goes to zero. The ELF x basins are thus expected to stand for the contribution of frontier orbitals to the electron pair regions in the studied systems. Within a frozen molecular orbital approximation (FMO), i.e. orbitals remaining unchanged upon addition and subtraction of one electron,  (resp.  ) reduces to the well known Fukui functions for nucleophilicity (resp. electrophilicity). [START_REF] Fuentealba | Applications of Topological Methods in Molecular Chemistry[END_REF][71][START_REF] Fukui | Orientation and Stereoselection[END_REF][START_REF] Yang | Proc. Natl. Acad. Sci[END_REF] The integration of the HOMO/LUMO density computed over the ELF x basin volumes provides populations for each nucleophile/electrophile basin. For instance, Figure 1 depicts the ELF x localization domains for the phenol molecule. The topological analysis yields valence basins accounting for eight electrophilic basins (red domains) and several nucleophilic basins (blue domains). The ortho-and meta-carbon atoms seem to be electrophilic while ortho-and para-carbon atoms as well as the oxygen lone-pairs have clearly a nucleophilic character. These findings are perfectly consistent with the well-known reactivity of the phenol: the electrophilic aromatic substitutions preferentially involve the ortho-and the parareactivity. [START_REF] Fievez | [END_REF]    

III. Theory.

In a reaction path where two molecular reactants MA and MB approach each other, the intermolecular interaction energy between MA and MB can be obtained from a supermolecular point of view [START_REF] Waller | Handbook of Computational Chemistry[END_REF] :

- - (1)
where is the total energy of the MA--MB complex, and, and correspond to the energies of reactants MA and MB calculated in their relevant isolated states, respectively.

To explore the potential energy surface (PES) of the complex conformations, the sum remains unchanged but the total energy needs to be computed again which entails an untenable computational cost for large or complex molecular system where local minima abound. Thus, the search of different local minima generally requires an adequate guess or previous knowledge of the reaction mechanism. It is interesting to explore another route in order to obtain the variation for different conformations with respect to the lowest energy structure.

Interestingly, the total energy can be directly calculated within the framework of the Interacting Quantum Atoms (IQA). [64,[START_REF] Maxwell | [END_REF] IQA provides a partition of the total energy of the complex into intra-atomic and inter-atomic energy terms o all di erent topological basins Ω:

where , , and are, respectively, the kinetic energy of the electrons, the electron-nuclear and electron-electron potential energies in one basin. and are, respectively, the electronnuclear, nuclear-nuclear and electron-electron potential energies when the first particle is in basin and the second in basin . This approach has been initially developed for the QTAIM partition. Some attempts to extend the IQA approach to the ELF partition have been also proposed. [66,67,77] Within the IQA methodology, -can be further partitioned from a supermolecular point of view as

- -
. Here,

is an intermolecular interaction energy that will be specified more accurately soon, and we define the total energies and of the reactants by aggregating all the intra-molecular terms and only the potential terms:

The sum is often considered as a self-energy quantity which corresponds to the internal basin energy. The other terms of equation ( 3) are the potential energies of the basins belonging to MA.

Even if the geometry of MA is frozen in the geometry of complex, generally differs from where the reactant is computed in its relevant isolated state. This difference can be considered as a deformation energy defined itself as . When the reactants are located very far from each other, goes to zero. Thereafter, we can rewrite as:

Thus, we obtain:

- - - (4)
Where -.

Note that equation (4) proposed in the QCT framework can be brought closer to the activation strain model of chemical reactivity proposed by Bickelhaupt et al. [78] Finally, we note that the variation of the supermolecular interaction energy int for different conformations of the MA--MB complex with respect to the lowest energy structure, is simply given by the sum:

int --

.

We can now move on the expression of the interaction energy -. From equations ( 2) and ( 3), this latter can be partitioned as:

where Ω Ω is the inter-nuclear repulsion between the nuclei of MA and MB , Ω Ω and Ω Ω are the electron-nuclear interactions between the basins of MA and MB, and finally, Ω Ω recovers the total interatomic electron-electron potential energy between the basins of MA and MB.

Ω Ω can be further split into coulomb and non-classical exchange-correlation term Ω Ω . Finally, -int can be written as the sum of coulomb and exchange-correlation energies:

- - - (6)
Where, being the distance between an electron in the basin Ω A and an electron in the basin Ω B , respectively. R A and R B are the nuclear locations of atoms Z A and Z B belonging to Ω and Ω domains.

It has long been well-known that the electrostatic contribution, here termed oul , accounts for a large fraction of the interaction energy in most intermolecular interactions even though other contributions are not necessary small or very small. [79][80][START_REF] Bertrand | Advances in Molecular and Cell Biology[END_REF][START_REF] Orangi | [END_REF][83] In the context of this work, it is interesting to quantify the role and the magnitude of the terms involved in equation ( 8). We review different examples to get a better understanding of the relative evolution of the three terms for various MA-MB intermolecular distances and different conformations when MA and MB interact. Two illustrative reactions have been selected. First, the formation of a typical dative bond in the borazane NH 3 BH 3 from NH 3 and BH 3 reactants and the formation of a hydrogen-bond between two water molecules leading to the formation of the canonical water dimer. For NH 3 BH 3 , various B-N distances have been selected (3, 4 and 5 Å) and two different conformations, displayed in Figure 2, have been considered. Three different conformations have been considered for the water dimer (H 2 O) 2 . The distance between the two oxygen atoms has been frozen to 4 Å. The results are gathered in Table 1 and the selected conformations are depicted in Figure 2. Coulomb and Exchange-correlation energies have been calculated with AIMALL software. [84] - A detailed analysis of results presented in Table 1 leads to the non-ambiguous pithy conclusions:

(a) - (b) - (c) - (c) MA = NH 3 ; MB = BH 3 d N-B = 5 Å GS -0.
1. The results confirm that the coulomb intermolecular interaction energy oul remains strongly correlated to int whatever the geometry of the selected molecular system. This will be further illustrated in the following section:

oul remains the driving term of the chemical reactivity when the two molecules interact. In which we used the fact that the total variation Δ because the total system is isolated.

We can reasonably assume that the electron density in the domains of MA will be mainly perturbed by variation of the number of electrons NA, i.e. the crossed derivatives such as or are expected to be negligible. Thus, from the expression of oul , we get:

Which reflects the interaction between the unchanged electron density and the response to changes in the number of electrons.

And then,

Where and are the Fukui functions [71] typically associated with reactive nucleophilic or electrophilic sites of the reactants. In addition, the link between Fukui functions and the frontier molecular orbital theory (FMO theory) [85] description of the reactivity arises when interacting reactants are considered as electron donors (here MA for example), or otherwise electron acceptors (here MB for example). For the sake of clarity, we will first consider the case in which MA is a pure nucleophile, while MB is a pure electrophile. Thus, by neglecting the orbital relaxation terms and .

Moreover, the total density can be partitioned as:

So that the term in equation (10) can be written i i

The second term of this latter sum corresponds to a stabilization of the LUMO energy of the electrophile due to the inner electronic density of the nucleophile. However, within the context of the FMO theory we assume that the frontier interactions, here , are the dominating terms. This typically matches with the topological partition of Fukui functions [86] or with the electron localization function ELF x which will used in the next sections of this work. Therefore, Equation ( 11) can be also expanded as:

Interestingly, the interaction between the charge domains can be rewritten from equation (12) as:

In the framework of the FMO theory, is the termed Dual-Descriptor of the molecule MB [87] . To conclude, it is clear from equation ( 12) that the variation of the intermolecular interaction energy between interacting reactants MA and MB is fully driven by E dual . [87] Practical calculations. Expressions (11) or (12) are not really appropriated for practical applications due to their computational cost. However, they can be numerically evaluated by means of a multipole expansion (ME). The latter requires precise computations of the charge distributions at any point of the molecular space. The ME is based on a Taylor expansion of the |R AB -(r A -r B )| -1 term, which leads to typical form: [88][89][90][91] where l l is an interaction tensor obtained from stable recurrence relations. The indexes are

computed from m A = -l A … l A (m B = -l B … l B ) and l A ... ∞ (l B ... ∞ .
are the distributed multipole moments on the topological partition. Using only the first terms of the ME (that is only the monopoles) and considering equation (11) in terms of electrophilic (Ω El ) and nucleophilic (Ω Nu ) domains, reads:

Where N Ω u l are the populations of nucleophile/electrophile domains. The populations are obtained from the usual condensation of the HOMO/LUMO density computed over the ELF x basin volumes:

where ω is a weight which still needs to be de ined.

And, and are the charges of the nucleophilic and domains, respectively.

It is worth noting that the way to condense a function is arbitrary as far as that an electron domain or the definition of an atom in a molecule remains arbitrary. [START_REF] Tiznado | [END_REF]92] In chemistry, one usually prefers to assign a given scalar to the weight in the condensation scheme . Thus, based on a equal footing treatment of nucleophilic and electrophilic molecular regions ω was assigned to or depending on the single or double occupation of the HOMO orbital.

A simplified expression can be also used when it is relevant to only consider the HOMO-LUMO electron stabilizing interaction between the MA electron donor and the MB electron acceptor (typically a Lewis acid-base reactivity):

This latter equation is typically well suited for the Lewis acid-base reactivity such as the H 3 -X 3 (X = H, F, Cl, Br, I) series where the N-bonding sche e is driven by the charge trans er ro the interaction between the lone pair located on the nitrogen and the e pty p orbital at the boron. [93] Interestingly, the latter equation ( 14) can be related to the maximum matching criterion previously defined [18] in terms of Fukui interactions in the whole molecular space. Except for the protonated basins which typically contain the hydrogen atoms, it is worth noting that nuclei do not belong to ELF x Ω and Ω domains and consequently, the terms Finally, we can extend equation ( 13) to the general case in which both molecules exhibit some nucleophilic and electrophilic sites distributed in their respective molecular space. This leads to the more comprehensive expression,

IV. Implementation and Algorithms.

In this section, we propose an original algorithm in order to perform an automated search of the favored chemical reaction paths corresponding to the minima of E dual . This latter can be either calculated from the equations ( 13), ( 14) or ( 15), depending on the studied reactivity. Figure 3 provides the flow diagram of the algorithm used. The automated search algorithm can be roughly split in four distinct steps:

Step 1. The ELF x localization domains. At beginning of the process, the ELF x nucleophilic (resp. electrophilic) basins are determined for each molecule taken in their isolated states by using the frontier orbitals of MA (resp. MB). After this initial process, the basins are obtained for each molecule separately taken in the geometry of the complex according to the step of splitting of the supermolecular MA-MB density described in the Step 4. The ELF x function is also computed by using the frontier orbitals of the supermolecule MA-MB. Last, the basin populations and are thereafter computed using the TopChem2 package. [94] Step 2. The rotation scheme. Controlling the orientation and the rotation of single molecules MA and MB with each other is key to ensure a full scan of the whole molecular space. As illustrated in Figure 4, we used a time-tested process based on the frozen center of mass separation of MA and MB and rotation angles. [95] For a frozen MA geometry, the process selects by carrying out all the small geometrical rotations of MB, the best relative orientation of MB associated to the lowest value of E dual .

All geometries of MA are tested, each can be defined with a set of rotation angles (θ, φ, ψ) around cartesian axis (Ox, Oy, Oz), respectively. The initial structure of MA is defined as (θ, φ, ψ) = (0, 0, 0). At the end of the process, the lowest value of E dual , termed , is obtained.

The optimization process of rotation steps ( , , dψ) is detailed in section I of the supplementary information. Step 3. Update the supermolecular MA-MB density. Once the optimized structure associated with is determined, a new single point DFT calculation of the supermolecule is computed and then, an updated supermolecular MA-MB density is obtained. Thereafter, the whole space is partitioned into MA and MB molecular domains as explained in the step 4.

Step 4. Splitting the supermolecular MA-MB density. In this step, we explain how the total electron density of the supermolecule MA-MB is split into MA and MB molecular densities, needed in Step 1 to compute the ELF x basins.

For a long time, several approaches were used to split the electron density into atomic contributions that has typically led to the population analysis. [96][97][98] In the context of our work, the electron density of the supermolecular system MA-MB can be properly split as:

n φ occ φ
Where  i (r) are the orbital expanded using atom-centered basis functions, µ (r). P µν are elements of the total electron density matrix defined as follows:

n i c i occ i c i
n i is the occupation of  i (r) and c iµ and c iν are the real expansion coefficients.

When MA and MB are very distant from one another, quickly goes to zero and the sum between ρ(r) MA and ρ(r) MB quickly goes to total density of the supermolecular system ρ(r). Overall, for intermolecular interactions, is expected to re ain s all with respect to the r) MA and r) MB density components, as illustrated in Figure 5(1) for the case of the formation of the covalent bond B-N bond in BH 3 NH 3 . The splitting of the supermolecular MA-MB density into MA and MB molecular densities then allows to compute the ELF x basins of MA and MB while preserving the ELF x topology of the supermolecule MA-MB. This is illustrated with the case of NO 2 + in interaction with the benzene molecule in the first step of the benzene nitration process. [99] We have separately calculated ELF x [ (r) NO2+ ] and ELF x [ (r) C6H6 ] from respectively the and densities obtained from the total density of the supermolecule.

As shown on Figure 5(2), the topologies of ELF x of each molecule NO 2 + or C 6 H 6 , remains almost unchanged with respect to the ELF x topology obtained from the supermolecular electron density. Thus, replacing by into the MA domains or in the same way, replacing by into the MB domains is a relevant approximation. In fact, this latter looks like to an improved well-known promolecular approximation [100,101] where the polarization of each molecule is nonetheless taken into account. Thus, the topologies of ELF x [ (r) MA ] and ELF x [ (r) MB ] are separately computed and then, we go back to Step 1.

V. Computational Details. The B3LYP hybrid functional level with the Gaussian09 software [102] was used for all calculations of total and intermolecular interaction energies. For NH 3 +BH 3 , the standard allelectron 6-31+G(d,p) was used for all atoms, while the 6-311+G(3df,2p) basis sets was used for the water dimer. The TopChem2 package [94] was used for all QCT analyses. The ELF x isosurfaces were displayed by means of the Molekel v4.3 software. [103] Coulomb and Exchange-correlation energies given in Table 1 have been calculated within the QTAIM partition using the AIMALL [84] software. Note that the total molecular energy obtained with the IQA partition scheme remains very close to the energy computed from Gaussian09 owing to the well-known compatibility between IQA and B3LYP. [START_REF] Maxwell | [END_REF] VI. Application to selected examples. In principle, our methodology can be applied to any chemical reaction. However, in this work we review only two examples where its applicability is illustrated and evaluated. First, we considered the typical dative covalent bond formation in the borazane NH 3 + BH 3 → NH 3 BH 3 . Thereafter, we have studied the water dimer, where the water molecules are held together by a typical hydrogen bond which remains a subject of intensive studies. [65] A. The borazane formation. Let us consider the case of the formation of the BH 3 NH 3 molecule. Its chemical path has long been known: Fujimoto et al. determined in 1974 that the covalent bond formation comes from a well-known electron transfer between the HOMO of NH 3 and the LUMO of BH 3 as an archetypical Lewis acid/base mechanism. [93,104] The robustness of our methodology has been evaluated for this reaction path, E dual being computed from equation (14). The nucleophilic sites of NH 3 and electrophilic sites of BH 3 were first determined according to the ELF x topology. We then explored the conformational space and we looked for the minima on the potential energy surface of E dual using the algorithm described in the previous section. For example, Figure 6 displays a two di ensional colored map of E dual , the distance between the centers of mass of NH 3 and BH 3 being frozen to 5 Å. For a given BH 3 orientation, the process selects the orientation of NH 3 associated to the lowest value of E dual .

All orientations of BH 3 are tested, each geometry being only defined by two rotation angles (θ, φ) On the one hand, the overall picture perfectly agrees with the well-known reactivity of this system since the expected minima between the Lewis-acid NH 3 and the Lewis-base BH 3 (Fig. 6c) is clearly identified on the surface. On the other hand, the displayed maps highlight an excellent mapping of E dual (Fig. 6a) and the DFT intermolecular interaction energy (Fig. 6b). Hence, the locations of critical points (minima and maxima) of E dual , notably the location of the minima, appear in excellent agreement with the DFT intermolecular interaction energy surface. We noted that the structure associated to the minimum is consistent with the topology of ELF x domains of BH 3 and NH 3 displayed on Figure 7. For BH 3 , this topological analysis of ELF x yields valence basins accounting mainly for two electrophilic basins symmetrically distributed around the boron atom (red domains with a population of 1 e). For NH 3 , it yields two nucleophilic domains (blue domains with a total population close to 2 e) associated to the nitrogen lone pair while additional small nucleophilic domains have been also found close to hydrogen atoms of NH 3 . From this spatial arrangement of basins in the molecular space, we predict the main stabilizing contribution between the boron electrophilic domains (red domains on Figure 7) and the nitrogen lone pairs of the Lewis-acid (blue domains on Figure 7).

Thereafter, the automated search algorithm aimed to follow the favored chemical reaction path was applied to the formation of the Borazane. We have looked for the evolution of along the path for different distances between the center of mass of NH 3 and BH 3 selected from 2.5 to 5.0 Å. Figure 8 gathers the obtained plots of (Figure 8b) and the evolution of DFT intermolecular interaction energy (Figure 8a) computed as the difference between the total energy of NH 3 BH 3 and the corresponding monomers energies calculated in their relevant isolated states. We tested two strategies: in the first one, the BH 3 and NH 3 moieties are frozen to their isolated geometries for all distances: they are flat irrespective of the fact that both are pyramidal in the borazane molecule (in black on Figure 8b). In the second one, we let both fragments relax: BH 3 and NH 3 can thus become pyramidal as the distance decreases (in red on Figure 8b). The analysis of the plots leads to a unambiguous conclusion: the evolution clearly matches the evolution of the DFT interaction energy. It decreases when the reactants become closer. When the pyramidalization of BH 3 is taken into account (red plot in Fig. 8b), the energy profile of is modified and becomes very similar to the DFT profile. These last results noticeably depict a remarkable mapping of the E dual onto the intermolecular interaction energy surface.

B. Formation of the hydrogen bond in the water dimer. Another example is the canonical water dimer. Numerous theoretical as well as experimental studies on the nature and the strength of hydrogen bonds (HB) can be found in the scientific literature, [65,[105][106][107][108][109] which testifies to the importance of this interaction in any field in chemistry and in biology. Thus, efforts to better understand HB interactions remain of interest. The canonical water dimer remains an archetypal example. An orbital analysis of the water dimer has shown that the HOMO-1 orbital must be considered in addition to the HOMO in the donor water molecule to correctly represent the situation of interaction between the two monomers. Thus, the nucleophilic basins of ELF x have been computed for x From this topology, it is easy to predict a stabilizing contribution between oxygen nucleophilic basins and protonated electrophilic ones and on the other hand, a repulsive contribution between oxygen lone-pairs themselves.

In search of the global and local minima of (H 2 O) 2 . We explored the conformational space of the interaction scheme between the water molecule and we looked for the minima on the potential energy surface of E dual computed from equation (15). As before, the distance between the two centers of mass is first fixed to 5 Å. Comparing the two maps highlights a faithful mapping of E dual and the DFT intermolecular interaction energy of the dimer water. The expected linear non-planar (so called canonical) minima is clearly identified on both surfaces. Because of the high symmetry of the dimer, each monomer having the same C 2v geometry, the linear non-planar minimum is found degenerated on several locations on the map. Moreover, several other critical points corresponding to the well-known typical structures are also identified. Figure 10c notably highlights the two conformers known as the bifurcated and cyclic stationary points of the water dimer. [110,111] The presence of these structures are confirmed on the DFT interaction energy surface where the same critical points are found. The results presented here demonstrate that energetic properties as well as structural parameters of the different geometries of the water dimer can be reproduced quickly and efficiently. In spite of numerous approximations used in this work, equation (15) used to compute E dual seems to be reasonably sufficient to seize not only the global minima but also structures higher in energy.

Thereafter, the automated search algorithm aimed to follow the favored chemical reaction path was applied for the formation of the hydrogen bond in the water dimer. We have looked for the evolution of along the path for different distances between the center of mass of monomers selected from 1.5 to 5.0 Å. Figure 11 VII. Concluding Remarks. By combining conceptual DFT and quantum chemical topology, we proposed here an original point of view leading to an efficient methodology able to describe and predict the chemical paths from the topological analysis of the modified electron localization function ELF x . We have shown, in a rigorous way, how the first-order variation in the intermolecular coulomb interaction energy expressed in terms of the response to changes in the number of electrons can be a relevant quantity to drive the chemical reactivity between the reactants. An automated reaction-path algorithm aimed to determine the most favorable relative orientations has been presented. This algorithm was applied to the formation of a dative covalent bond in NH 3 -BH 3 and to the analysis of the water dimer minima bound by a weak hydrogen bond. For both systems, our approach unveils a noticeable mimicking of E dual onto the DFT intermolecular interaction energy . We show that the global minima, but also structures higher in energy, can be clearly identified in order to predict the chemical reaction paths. As we have only presented a few examples, future works will focus on the description and the prediction of a large panel of chemical reactions. We are also currently exploring an improved expression of ELF x constructed from orbital-weighted Fukui functions previously proposed elsewhere. [112][113][114] Additional Supporting Information may be found in the online version of this article.

Figure 1 .

 1 Figure 1. ELF x localization domains and their populations obtained for the phenol optimized at the B3LYP/6-31+G(d,p) level of theory. Color Code: Red: electrophilic regions; Blue: nucleophilic regions.

  of the nuclei of one molecule with the charge domains of the other molecule. Finally, we go back to the variation of the supermolecular interaction energy -

  (a) Intermolecular interaction energy (kcal/mol) calculated with respect to the isolated reactants frozen in the geometry of complex. (b) Deformation energy calculated from Equation (8). (c) Coulomb and Exchange-correlation energies have been evaluated within the QTAIM partition.

Figure 2 .

 2 Figure 2. NH 3 BH 3 and (H 2 O) 2 selected conformations.

2 .

 2 When the reactants are located far from each other, the repulsive deformation energy de tends to be fully compensated by the attractive exchange-correlation part of the interaction energy. For shorter MA-MB distances, this fundamental behavior seems to be partially sustained and thus, explains why oul remains the dominant contribution when reactants are approaching each others. As these examples clearly confirm that the coulomb interaction largely contribute to the driving force of the intermolecular interaction between MA and MB, we will now see how it can be estimated from cDFT to predict reaction pathways.Chemical reactivity theory applied to the IQA coulomb interaction term.A relationship between the chemical DFT description and the quantum chemical topology can be explored by the calculation of the first-order variation in the coulomb intermolecular interaction energy expressed in terms of the response to changes in the number of electrons Δ A or Δ B where the external potential remains unchanged. Thus,

  the time zero.

Figure 3 .

 3 Figure 3. Flow diagram of the automated search algorithm of favored chemical reaction paths.

Figure 4 -

 4 Figure 4 -Illustrative example of the rotation scheme described in the step 2. For a frozen geometry MA (here NH 3 ), the best relative orientation of MB (here BH 3 ) is selected. All geometries of MA are tested, each can be defined with a set o rotation angles θ φ ψ .

Figure 5 . ( 1 )

 51 Figure 5. (1) B3LYP/6-31+G(d, p) plots of the electron density of NH 3 BH 3 (e.bohr -3 ) oriented along the x axis between N and B. The N-B distance was frozen to 2 . Some electron density contributions given in equation (16). Color code: black: total electron density r); blue: r) NH3 and r) BH3 ; red: r) NH3-BH3 (2) ELF x,HOMO localization domains of NO 2 + in interaction with the benzene molecule in the first step of the benzene nitration process calculated at the B3LYP/6-31+G(d, p) level of theory. The distance between the mass centers of NO 2 + and C 6 H 6 was moved away from 3 . (a) ELF x,HOMO is only calculated ro r) NO2 . (b) ELF x,HOMO is only calculated ro r) C6H6 . (c) ELF x,HOMO is calculated ro the total electron density r).

  around each cartesian axis (Ox, Oy) since the rotation around the O z axis, oriented along the N-B bond, has been frozen. The obtained map E dual is compared with the corresponding map of the DFT intermolecular interaction energy computed from the relevant isolated states of NH 3 and BH 3 reactants.

Figure 6 .

 6 Figure 6. Two-dimensional map E dual ( (a) vs. the DFT intermolecular interaction energy (b) surfaces for NH 3 BH 3 calculated at the B3LYP/6-31+G(d, p) level of theory. The distance between the centers of mass of NH 3 and BH 3 is fixed at 5 Å. θ and φ correspond to the orientation o the 3 molecule around the x and y axis (see text for details). (c) Coordinates (θ, φ) of the minimum the both surfaces.

Figure 7 .

 7 Figure 7. ELF x localization domains and their populations for NH 3 and BH 3 calculated at the B3LYP/6-31+G(d, p) level of theory. Color Code: blue: nucleophilic regions; red: electrophilic regions. (a) Localization domains of the optimized reactants taken in their isolated states (b) Localization domains of reactants obtained from the splitting of the supermolecular electron density, the distance between the centers of mass of NH 3 and BH 3 being frozen to 5 Å.

Figure 8 .

 8 Comparative study of the variation of (a.u.) vs. the DFT intermolecular interaction energy E 0 int (a.u.) calculated at the B3LYP/6-31+G(d, p) level of theory according to the distance between the center of mass of NH 3 and BH 3 . On the E dual picture (b), the black plot corresponds to a constrained planar geometry of BH 3 whereas the red plot corresponds to the partially optimized BH 3 .

22 9

 22 depicts the ELF x localization domains for a single optimized water molecule and for the water molecules embedded in the water dimer, the centers of mass of molecules being moved away from 5 Å.For both cases, the topological analysis of each water molecule remains unchanged and yields valence basins accounting for two electrophilic basins close to the hydrogen atoms (red domains) and two nucleophilic (blue domains) for the oxygen lone pair basins.

Figure 9 .

 9 Figure 9. Main localization domains of ELF x and their populations for the water molecules computed at the B3LYP/6-311+G(3df, 2p) level of theory. Color Code: red: Electrophilic regions; blue: Nucleophilic regions. (a) Localization domains of a single optimized water molecule (b) Localization domains of water molecules obtained from the splitting of the electron density of the water dimer, the centers of mass of two water molecules being moved away from 5 Å. The dimer is taken in the initial geometry used to build the map of E dual (see next part).

Figure 10

 10 displays the map of E dual together with the corresponding surface of the DFT intermolecular interaction energy computed as the difference between the total energy of the dimer system and the monomers calculated in their relevant isolated states. For a given orientation of one H 2 O molecule (on the left in Figure 10c), the process selects the best orientation of the other H 2 O molecule associated to the lowest value of E dual . All orientations of the frozen H 2 O are tested, each geometry being only defined by two rotation angles (θ, φ) around each cartesian axis (Ox, Oy) since the rotation around the O z axis, oriented along the O-O bond, has been frozen. For each rotation step , the structural parameters associated with E dual have been extracted.

Figure 10 .

 10 Figure 10. Two-dimensional map E dual ( (a) vs. the DFT intermolecular interaction energy (b) surfaces for the water dimer calculated at the B3LYP/6-311+G(3df, 2p) level of theory. The distance between the centers of mass is fixed at 5 Å. θ and φ correspond to the orientation o the le t water olecule denoted by around the x and y axis (see text for details). (c) Coordinates (θ, φ) of structures found on the both maps.

  gathers the obtained plot of compared to the evolution of DFT intermolecular interaction energy . Comparing both plots reveals that the evolution clearly matches the evolution of the DFT interaction energy. It decreases when the reactants become closer and it reaches a minimum for the equilibrium distance (close to 3 Å). When the distance is further reduced, increases. In the same way as for NH 3 BH 3 , these latter results noticeably reveal a remarkable isomorphic mapping of the E dual onto the intermolecular interaction energy .

Figure 11 .

 11 Figure 11. Intermolecular interaction energy (a) calculated at the B3LYP/6-311+G(3df, 2p) level of theory v.s (b) with respect to the distance between the centers of mass of each water molecule.

Table 1 .

 1 B3LYP

		461	0.324	-0.321	-0.464
	Conformation A	0.029	0.419	-0.396	0.007
	Δ=conf A -GS	0.490	0.095	-0.075	0.471
	d N-B = 4 Å
	GS	-1.318	2.023	-2.010	-1.331
	Conformation A	0.371	6.076	-5.894	0.189
	Δ=conf A -GS	1.689	4.053	-3.883	1.520
	d N-B = 3Å
	GS	-5.681	11.801	-13.621	-3.861
	Conformation A	10.044	41.242	-33.960	2.762
	Δ=conf A -GS	15.725	29.441	-20.339	6.623
	MA = H 2 O ; MB = H 2 O ; d O-O = 4 Å
	GS	-2.461	1.779	-1.654	-2.586
	Conformation A	-1.947	1.531	-1.361	-2.117
	Conformation B	1.683	0.907	-0.730	1.506
	Δ=conf A -GS	0.514	-0.248	0.293	0.469
	Δ=conf B -GS	4.144	-0.872	0.924	4.092

/6-31+G(d, p) intermolecular interaction energies.