Electron sharing indexes at the correlated level. Application to aromaticity calculations

Electron sharing indexes (ESI) have been applied to numerous bonding situations to provide an insight into the nature of the molecular electronic structures. Some of the most popular ESI given in the literature, namely, the delocalization index (DI), defined in the context of the quantum theory of atoms in molecules (QTAIM), and the Fuzzy-Atom bond order (FBO), are here calculated at a correlated level for a wide set of molecules. Both approaches are based on the same quantity, the exchange–correlation density, to recover the electron sharing extent, and their differences lie in the definition of an atom in a molecule. In addition, while FBO atomic regions enable accurate and fast integrations, QTAIM definition of an atom leads to atomic domains that occasionally make the integration over these ones rather cumbersome. Besides, when working with a many-body wavefunction one can decide whether to calculate the ESI from first-order density matrices, or from second-order ones. The former way is usually preferred, since it avoids the calculation of the second-order density matrix, which is difficult to handle. Results from both definitions are discussed. Although these indexes are quite similar in their definition and give similar descriptions, when analyzed in greater detail, they reproduce different features of the bonding. In this manuscript DI is shown to explain certain bonding situations that FBO fails to cope with. Finally, these indexes are applied to the description of the aromaticity, through the aromatic fluctuation (FLU) and the para-DI (PDI) indexes. FLU and PDI indexes have been successfully applied using the DI measures, but other ESI based on other partitions such as Fuzzy-Atom can be used. The results provided in this manuscript for carbon skeleton molecules encourage the use of FBO for FLU and PDI indexes even at the correlated level.


Introduction
The analysis of electronic distribution in molecules is best understood with the concept of bond order (BO). Its definition is strongly related with the seminal work of Coulson, 1 who provided one of the first attempts to quantify the extent of electronic sharing between two atoms by means of quantum mechanics calculations. It is worth noticing that in this work Coulson had the intention of using these bonds of fractional order-as he stated-to study some polyenes and aromatic molecules. It is thus of remarkable importance that, already from the very beginning, the concepts of aromaticity and bond order were tightly connected.
From then on, the study of bond orders followed its own way, gaining a prominent situation in the descriptors used for the analysis of the electronic structure of molecules. In order to avoid controversy, since not all indexes that account for the electron sharing of a couple of atoms can be regarded as bond orders, we prefer the term electron sharing index (ESI) to comprise them all. Several ESI have been put forward in the literature since Coulson's work. It is not our intention to make here an extensive quotation of all papers treating chemical bonding by the use of ESI, and therefore we will just comment on the relevant works for our present purposes. In this sense, it is worth noticing that for a long time Mulliken's analysis 2 (and its extension to bond orders and valences) [3][4][5][6][7] has been one of the preferred for its simplicity and inexpensive computational cost. On the other hand, since the work of Wiberg 8 there has been a growing interest on ESI based on the exchange-correlation density, 9 the major contributions being those reported for Atoms-in-Molecules 10 and Fuzzy-Atom 11 partitions.
Aromaticity was discovered as a multifold property whose characterization was performed by means of other criteria such as typical magnetic manifestations, 12 bond length equalization, 13 aromatic stabilization energies 14 or more occasionally by characteristic electron delocalization. 15 It is widely known that all these manifestations individually are insufficient for the characterization of the aromaticity of a given species. Since aromaticity is said to be multifold, several aromaticity manifestations need to be taken into account simultaneously to finally assess the aromaticity of a certain compound. 16 Furthermore, the need for continuous examination of existing aromaticity indexes to recover situations whose aromatic character is perfectly known can provide further evidence of deficiencies and weak points of current aromaticity measures to ensure a correct description of aromaticity in molecular systems.
The goal of this work is thus twofold. First of all, we will review the concept of an ESI from Fuzzy-Atom and QTAIM partitions, starting from the calculation of correlated first-and second-order density matrices; while its adequacy to describe electronic structure in molecules is discussed. And secondly, the use of ESI for the definition of aromaticity indexes is revised, while it is tested for correlated calculations. The present manuscript is devoted to open a discussion on ESI and aromaticity in a conference named after Michael Faraday, who in 1825 discovered the first aromatic molecule, benzene. 17

Electron Sharing Indexes (ESI)
The concept of bond order, or ESI, 18,19 according to Fulton's terminology, has gained a prominent role in the quest for the best tools to characterize the electronic structure and chemical bonding of molecular systems. Since the seminal work of Wiberg, 8 different indexes have been put forward in the literature and, despite there is no unique definition available, many of them have a common pattern: they measure the extent the electrons are shared by two (or more) entities, usually atoms. In this way, well-known quantities such as the delocalization index (DI), 20,21 Fuzzy-Atom bond order (FBO), 11 Mayer Bond Order 3-6 or Wiberg Index 8 (MWI), among others, may be included in this definition and can reproduce to some extent the Lewis structures. 22 These ESI can be classified according to the way they define an atom within the molecule. In the pioneer works, 2,7,8,23 the atoms were formally defined by the Hilbert-space basis functions assigned to them. Obviously, this definition is restricted to the LCAO-MO approximation and may lead to severe difficulties if basis functions lacking atomic character are included, such as diffuse or bond functions. Another way is to partition the 3D physical space into regions assigned in one way or another to each atom. These atomic regions or domains have usually sharp boundaries, like the Voronoi cells, 24 the muffin-tin spheres, 25 the Wigner-Seitz cells, 26 the Daudel loges, 27 the ELF basins, 28 or Bader's atomic basins; 10 however, they can also ''exhibit a continuous transition from one to another'', 11 like the Fuzzy-Atoms. The 3D space atomic definition is more general, in the sense that it can be potentially applied beyond the LCAO, contrary to the Hilbert-space case. Since Hilbert-space approaches have been widely studied, in the present work, we will only focus on the family of ESI defined on the 3D-space.
Since the advent of the QTAIM 10 (quantum theory of atoms in molecules) the preferred way of most researchers to define an atom in a molecule has been the Bader atomic basin: a region of the space containing an attractor surrounded by a zero-flux surface or by infinity. Such atomic regions are non-interpenetrating 10 and can often exhibit complicated shapes, which makes the integration over the atomic domain rather cumbersome.
In contrast, a Fuzzy-Atom partitioning of the space provides atomic domains that can overlap and hence part of the physical space is assumed to be shared to some extent by two or more atoms. There are many ways to define the Fuzzy-Atomic domains. All of them are described by a non-negative weight function w A (r) defined for each atom A at every point of the space r in such a way that w A (r) is close to one in the vicinity of the nucleus of atom A and gradually tends to zero as the distance to the given nucleus increases.
However, the numerical integrations to be carried out within the Fuzzy-Atom framework are noticeably less expensive and straightforward than those involving the QTAIM atomic basins. On the other hand, there is indeed some arbitrariness concerning the definition of the weight factors that characterize the Fuzzy-Atoms. In previous studies 11,29,30 we have used Becke's 31 method for multicentric integration, where the weight factors have been obtained through an algebraic function for each atom that is equal to one in the vicinity of the respective nucleus and progressively decreases to zero with distance. This function depends upon some set of atomic radii that control the size of the Voronoi cells 24 and upon a stiffness parameter that controls the ''shape'' of the atom.
We have also tested other weight functions, particularly by adjusting the ratio of atomic radii using the position of the extremum of the density along the line connecting the atoms or even those derived from Hirshfeld's original idea of promolecular densities. Numerical experience 11,30 seems to indicate that reasonable results can be obtained using the simplest Becke's scheme and a balanced set of covalent atomic radii, like those from Suresh and Koga 32 or Slater. 33 The results for covalent unpolarized bonds are very similar to those obtained using a QTAIM partitioning of the space. For ionic systems or bonds involving B, Be or Al atoms the deviation is more important, basically because the relative atomic sizes predicted by the QTAIM partitioning are very different from those used to define the Fuzzy-Atoms. As a result the chemical bonds are apparently strongly polarized in QTAIM, 34 presenting exaggerated partial atomic charges but also the side effect of a much lower bond order, according to the predicted ionicity of the bond.
In the present work we have also explored the possibility of using the topology of the electron density to determine the set of atomic radii for each molecule to be used in Becke's integration scheme. As will be further discussed, typical features of the QTAIM-based ESI can be easily recovered using this mixed scheme with lower computational cost.
Even though there is a large amount of literature behind the concept of ESI, the most accepted approaches to the calculation of such quantities are somehow based on the Ruedenberg's generalized exchange density term. 9 Originally introduced as a decomposition of the pair density, this term is more generally known as the exchange-correlation density (XCD). The physical meaning of the XCD becomes evident if such function is defined as the difference between the pair density of the system and a fictitious pair distribution constructed from two independent oneparticle distributions, where g (1) (x 1 ) and g (2) (x 1 ,x 2 ) are the diagonal parts of the first-(1-RDM) and the second-order (2-RDM) reduced density matrices, respectively. They are commonly known simply as the density and the pair density.
In the MO approximation, both the pair density and the density can be written in terms of the molecular spin orbitals as and where G kl ij and G j i expand the second-and first-order density matrices. For wavefunctions expanded in terms of Slater determinants, the 1-RDM and 2-RDM are computed from the coefficients of the wavefunction expansion. For monodeterminantal wavefunctions (HF or DFT within the Kohn-Sham formalism) these density matrices are diagonal at any order, and therefore the 2-RDM can be expressed in terms of the 1-RDM; that is, the XCD can be written entirely in terms of the 1-RDM: This formula can be used for correlated wave functions too, even though, as an approximation to the true XCD. It will be referred hereafter as HF-like exchangecorrelation density (HF-XCD). The electron density contains information about the position of a single electron regardless the position of the other N À 1 electrons. In the same way, the pair density is a two-electron quantity that accounts for the distribution of an electron pair, irrespective of the position of the remaining N À 2 ones. Hence, the XCD carries the information about the motion of a pair of electrons contained in the pair density, while it avoids that of the individual distribution of the electrons inherent in the one-particle density. Such a function is expected to be close to zero for two points distant in space and thus it is not surprising that it integrates to N, the number of auto-interactions artificially contained in the independent electron pair distribution, It is easy to understand why the XCD has become so popular in the definition of ESI. The integration of the XCD is decomposed over all pairs of atomic domains providing a measure of the extent of correlation (interaction) between a given pair of regions. Indeed, back in 1975, Bader and Stephens 21 pointed out that the integration of the XCD upon two different regions of space A and B leads to a ''measure of the correlative interactions between the electrons in these two moieties'': which at the same time is a measure of the extent to which electrons are delocalized over these regions. 20

ESI based on the HF-XCD and XCD
The first ESI based on the QTAIM space partitioning was introduced by Cioslowski and Mixon in 1991. 35 They defined the so-called covalent bond order, which was based on localized orbitals obtained by an isopycnic transformation, 36,37 and whose formula reads: where the S ii (A) is the diagonal element of the atomic overlap matrix (AOM) of the localized spin orbitals integrated over the domain of atom A so that P i l 2 i S ii ðAÞ 2 is maximized, and l i is the corresponding occupancy. These localized spin orbitals must be related to the original natural spin orbitals through a unitary transformation, which guarantees that the same density is reproduced from both (isos pyknos).
Two years later, Fulton 19 proposed another ESI in the framework of the QTAIM by using the 1-RDM in following manner: whereupon the pseudo-square root of the 1-RDM reads l k and Z k (x) being the natural occupancies and natural spin orbitals, respectively. It is more convenient to express eqn (8) in the following manner: where the S ij (A) are the elements of the overlap matrix of the natural spin orbitals integrated over the domain of atom A.
One year later, 38 Á ngya´n et al. proposed a mapping between the Mayer's Hilbertspace based bond order and the QTAIM formalism that yet leads to a different ESI, which was again using the HF-XCD, and turned out to coincide with Bader's definition eqn (6) for single-determinant wavefunctions: As a curiosity, it is worth pointing out that Fulton originally proposed eqn (11), but he rejected it afterwards because with this formulation the sum of all electron sharing indexes for a given atom A is not proportional to its population, N(A). It is easy to prove this relationship is fulfilled in eqn (10). Finally in 1999 Fradera et al. 20 reexamined Bader's 1975 original work based on the XCD 21,39 to use it explicitly as an ESI, which nowadays is known as the delocalization index (DI): where the 1-RDM and 2-RDM are integrated over the atomic domains to provide the atomic and atomic pair populations, respectively. More recently, Bochicchio et al. 40 have studied in detail Fradera's DI and have showed that in the case of a correlated wavefunction it can be expressed as the sum of two terms, one of each being Á ngya´n's ESI 38,41 (which they call exchange term) and the other coming entirely from the second order density matrix (cumulant term): The cumulant term measures the difference between the true XCD and the HF-XCD one; hence for monodeterminantal wavefunctions d cum (A,B) = 0. In most cases the magnitude of the cumulant contribution is rather small. However, in a recent work by Bochicchio et al. 40 high cumulant contributions have been reported, which unfortunately could not be reproduced in the present work (vide infra).
It can be shown easily that Fulton's, Á ngya´n's and Fradera's definitions coincide for single-determinant wavefunctions. It is only interesting to see how they perform at the correlated level and what is the effect of using the true XCD instead of the HF-XCD approximated one.
On the other hand, all those ESI have always been computed within the framework of the QTAIM theory. However, any other 3D-space partitioning could in principle be applied. In the present work, we explore the possibility of computing such ESI at the correlated level in the framework of the Fuzzy-Atoms, for which the necessary numerical integrations are much more efficient from a computational point of view.

Statistical interpretation of the ESI
Another interesting property of the ESI from eqn (10) to (13) is that half of its sameatom contribution can be regarded as the number of electrons localized in the atomic basin; according to Fradera et al. this quantity is called localization index (LI). 20 Taking into account that the expected number of electrons in the basin of atom A, the population of A, is calculated as: its variance reads, which some authors recognize as the atomic valence for closed-shell systems. 7 Thus, the LI for an atom is the expected number of electrons minus its uncertainty (variance) in its atomic basin. On the other hand, recalling the formula for the covariance of different atoms population, one can easily see that ESI between two atoms A and B is proportional to the covariance of their electron populations. Finally, the sum of all ESI values involving a given atom, X may be regarded as a valence of the atom in a closed-shell system, or the number of electrons of atom A shared with the rest of the atoms in the molecule. This statistical interpretation holds for any ESI given in eqn (10) to (13) calculated through a monodeterminantal wavefunction. In addition, for the ESI in eqn (13) this interpretation is also valid for many-body wavefunctions, since the latter is based on the true XCD, and therefore on the true pair density. From the previous analysis it can be deduced that all ESI discussed in this work are non-negative quantities. From eqn (4) and (12), it is clear that the ESI thus defined cannot be negative for a single determinant wavefunction. For correlated wavefunctions the same is true for Cioslowski's, Fulton's and Á ngya´n's because they also come from the square of a real value. Finally, the ESI as defined in eqn (12) for a given pair of atoms A and B must also be non-negative because, when the population of one atom N(A) increases it is necessarily at the expense of the decrease of the population on the other one N(B). Hence the covariance 42 of both quantities is negative. To put it in a nutshell, all ESI reported here are non-negative defined.

ESI spin cases
If the ESI is constructed entirely from the 1-RDM it can be decomposed according to the two spin cases the 1-RDM is built of: alpha and beta. These two contributions will only differ for an open-shell calculation. On the contrary, when the ESI is calculated from the 2-RDM, it can be decomposed into three contributions, g ss xc ðr 1 ;r 2 Þdr 1 dr 2 and accounting for the three spin cases (a, a), (b, b) and (a, b) that are included in the true XCD. The first two are equivalent for a restricted closed-shell calculation. Therefore, it can be interesting to decompose the ESI into the same spin, d ss (A,B) + d s 0 s 0 (A,B), and different spin, d ss 0 (A,B), contributions. The natural names for these quantities are Fermi-ESI and Coulomb-ESI, respectively, named after Fermi and Coulomb pair densities. Although the ESI is a strictly positive quantity, it is worth noticing that ESI ss 0 can be negative, since the Coulomb pair density is not necessarily defined positive.
Having introduced the spin decomposition, let us focus on the sum rules fulfilled by each aforementioned ESI. The most general sum rule is as follows: where an exhaustive partition of the molecular space {A i } i=1...n has been supposed. It is easily shown that eqn (21) must be fulfilled by all indexes from (10) to (13) (13)). This is one of the shortcomings of the d A (A,B) index: for many-body calculations it does not provide a decomposition of the N electrons in the system. As will be shown below, the number of electrons missing to count up to N corresponds roughly to the unlike spin pairs. In order to avoid this shortcoming, Bochicchio et al. 40 introduced the unshared population, l A (A). The unshared populations correspond to those LI which would accomplish eqn (21) when using d A (A,B) instead of d(A,B), say: It is worth noticing that these l A (A) or unshared populations no longer fulfil eqn (14). Additionally, some sum rules are obeyed by the spin cases in which d(A,B) is decomposed according to eqn (19): Therefore, one can expect d ss 0 (A,B) to vary from negative to positive values, while d ss (A,B) is always non-negative.

The density relaxation
Any property can be defined as either an expectation value of the corresponding operator or as a response of the molecule with respect to a perturbation (a derivative). The Hellmann-Feynman theorem ensures that both methods lead to the same result. This is true for exact wavefunctions, and some approximate theoretical methods, such as the self-consistent field approach. Unhopefully, Hellmann-Feynman theorem does not hold for most many-body methods, with the exception of CASSCF and Full-CI calculations, being the latter prohibitive for most molecular systems. This fact leads to the ambiguity of whether calculating a given property from the corresponding energy derivative or from its expectation value. In this sense, the density matrices can be determined either as the expectation value of the proper Dirac's delta based operator or as a wavefunction response. The density matrix obtained in the latter case is known as response density or relaxed density, since it includes the effects due to orbital relaxation within the coupled-perturbed Hartree-Fock (CPHF) calculations.
One of the epistemological problems around the relaxed density (in practice, the wavefunction which reproduces such density) is that the density obtained in this way is not N-representable. As a consequence, the population of the natural orbitals can be larger than one or negative. Such a situation is especially uncomfortable for Fulton's ESI, where the square root of the natural occupancies appears. With this aim, to avoid this physically unsound situation, we believe one should always use the unrelaxed density for the calculations of ESI.
The relaxed 2-RDM has also been defined, 43 but unfortunately they are generally not available in standard computational packages. On the other hand, the second order reduced density matrices for methods such as MP2 and CCSD are not straightforwardly computed, since N-electron wavefunctions cannot be produced from these methods. Although energy-derivative second-order reduced density matrices can be obtained, its use is discouraged because such 2-RDM is not Nrepresentable inasmuch as it does not add up to the number of electron pairs. 44 Therefore, exact ESI as defined in eqn (12) are only obtainable from HF, CI and CASSCF calculations. This fact, together with the computational expenses needed for CASSCF calculations, has led us to restrict our study to the CISD results presented here, whereas CASSCF calculations have only been calculated for the systems where such a level of theory is needed to properly describe it (vide infra). It is worth mentioning in this context that 2-RDM have only been computed for CISD calculations.

Aromaticity indexes from ESI
Despite the electron delocalization is one of the key features of aromatic compounds, there have been few attempts to give quantitative aromaticity measures from ESI. Since these indexes measure the electron sharing between pairs of atoms, one would expect to obtain from them a reliable measure of the degree of delocalization in a given ring. Recently, two aromaticity measures based upon ESI have been put forward, the para-DI (PDI) 47

and the aromatic fluctuation index (FLU). 45
The para-delocalization index (PDI) Bader and coworkers 46 reported that for QTAIM-ESI the extent of electron sharing in benzene is greater for para-related carbons (para-ESI) than for meta-related ones (meta-ESI). Poater et al. 47 undertook an analysis of the QTAIM para-ESI in a series of polycyclic aromatic hydrocarbons (PAH) with six-membered rings (6-MRs). They concluded a close relationship between PDI and other measures of aromaticity such as ASE, HOMA or NICS, thus validating PDI as a reliable measure of aromaticity in 6-MRs.

The aromatic fluctuation index (FLU)
The FLU index 45 is constructed following HOMA 48 philosophy, i.e., measuring divergences from certain aromatic molecules chosen as references. The formula given for FLU reads: where the summation runs over all adjacent pairs of atoms around the ring, n is equal to the number of atoms of the ring and d ref (A,B) is the ESI of the atomic pair A and B for the aromatic molecule chosen as a reference, and a is a constant to ensure that the ratio of atomic valences is greater than one, so that this ratio penalizes species with high electron localization, which leads to lesser delocalization in the ring. The reference ESI (QTAIM, 45 Fuzzy-Atom 29 ) values for C-C and C-N bonds were obtained from benzene (1.4e, 1.4e) and pyridine (1.2e, 1.6e), respectively, at the HF/6-31G* level of theory. The larger the difference, the larger FLU; ergo, low FLU values correspond to aromatic molecules. It is fair to say that Bird 49 compared Gordy bond orders of all bonded pairs in a given ring to define a measure of aromaticity that has some resemblance to our FLU index. In this context it is also worth mentioning the work of Matta et al. 50,51 who have also attempted to construct an HOMA-like 48 index from QTAIM by substituting the bond length by the total electron delocalization. Although Bird's and Matta's indexes resemble FLU, and actually should be recovering the same aromaticity characteristics, one should notice the three indexes are indeed different. The differential feature of the FLU index is that it weights the divergence of electron sharing for each pair of bonded atoms with a quantity that accounts for electron localization cf. eqn (25). Interestingly, another aromaticity criteria has recently been put forward, the multicenter index (MCI) of Bultinck et al., 52 which is also appealing as an electronic measure of aromaticity.

Computational details
A wide set of molecules (H 2 , N 2 , F 2 , LiF, CO, CN À , NO + , He 2 , Ar 2 , LiH, HF, NaH, HCl, NH 3 , CH 4 , H 2 O, CO 2 , B 2 H 6 , C 6 H 6 and Mg 3 2À in its lowest-lying singlet electronic state) have been studied to analyze the effect of the inclusion of electron correlation and the differences between the different indexes in eqn (8)- (13). For this set of molecules geometry optimizations and correlated wavefunction calculations were obtained with Gaussian03 53 at the CISD/6-311++G(2d,2p) level of theory using Cartesian d and f orbitals (6d, 10f). Only the valence electrons were correlated (frozen core). In the case of benzene, a smaller basis set, 6-311G(d,p), has been used. For the sake of clarity, the nomenclature of the atoms for B 2 H 6 has been given in Scheme 1.
Unrelaxed density matrices may be obtained from the Gaussian03 package with the keywords density = rhoci and use = l916. A program designed in our laboratory 54 was used to generate the 2-RDM and 1-RDM matrices from the coefficients of the CISD expansion (the matrices thus computed are always unrelaxed). In the present work, only coefficients with values larger than 10 À15 were considered, and the elements of 2-RDM greater than 10 À12 were used for subsequent calculations. 55 The atomic overlap matrices for AIM and Fuzzy-Atom domains were generated with AIMPAC 56 and FUZZY 57 programs, respectively. For the latter, the Becke's multicenter integration algorithm with Chebyshev and Lebedev radial and angular quadratures has been used. A grid of 46 radial by 140 angular points per atom has been used in all cases. We have employed the recommended stiffness parameter k = 3 and adjusted the size of the Voronoi cells with the set of atomic radii determined by Koga, 32 except otherwise indicated. Finally, other software of our group 54 was used to generate the different ESI from the atomic overlap matrices and the corresponding density matrices.
The aromatic molecules studied consist of a subset of that already studied in ref. 29 and 45. Due to the computational cost of CISD calculations we have limited ourselves to molecules containing a maximum of three rings, that is, those of Scheme 2. For this set of molecules, geometry optimization and PDI and FLU calculations have been performed at the CISD/6-31G(d) level of theory. For such systems the calculation of the second order density matrix is prohibitive, and therefore only

Results and discussion ESI at correlated level
Since the effects of inclusion of the electron correlation in the topology of the density have already been addressed in the literature, 58,59 we will not discuss this topic further. Rather, we will focus on the influence that the Fermi and Coulomb correlation have on the ESI under study. In order to purely evaluate the effects of inclusion of the electron correlation on the wavefunction, reference HF and B3LYP calculations have been performed at the CISD/6-311++G(2d,2p) geometry with the same basis set. Tables 1a and 1b contain the calculations for the ESI with QTAIM molecular space decomposition, whereas Tables 2a and 2b contain the Fuzzy-Atom counterpart. Comparison of the QTAIM-and Fuzzy-Atom-ESI shows that for equally- a Inaccurate PROAIM integration (DN = 1.1e) the QTAIM approach, HCl being an exception to this general trend. Interestingly, for the isoelectronic series N 2 , NO + , CN À , and CO, the differences in d(A,B) computed from Fuzzy-Atom and QTAIM partitions widen with the increase in the electronegativity differences between the two bonded atoms (from 0.081 e in N 2 to   species. Some authors have provided support to the idea that the QTAIM partitioning exaggerates electron densities for electronegative atoms, 34 although others have opposed this view. [60][61][62] Our results indicate that, as compared to the Fuzzy-Atom approach, QTAIM overestimates the charge separation in polar and ionic bonds. However, both charges and bond orders being non observable quantities, one cannot tell whether QTAIM results are better from a chemical point of view than our particular Fuzzy-Atom approach ones to partition the space or vice versa. Let us focus now on the results for single-determinant wavefunction methods. As already reported for QTAIM, 63 in spite of the correlation included within B3LYP, d B3LYP (A,B) results are much closer to d HF (A,B) than to CISD ones, which leads us to notice that in general Coulomb correlation is not fully reflected in these d B3LYP (A,B) calculations, for both QTAIM and Fuzzy-Atom partitions. If we were about to decide whether d HF (A,B) or d B3LYP (A,B) is closer to the CISD one, based on a least squares fitting procedure, one concludes that HF is slightly closer to CISD than B3LYP for QTAIM, and just the opposite for Fuzzy-Atom. The reader may see how d B3LYP (A,B) values for the QTAIM partition are surprisingly larger for LiH, NaH, H 2 O or CO 2 , whereas CISD and HF values almost coincide. Such a finding is not reproduced by Fuzzy-Atom partition, whose B3LYP results are actually much closer to the CISD than the HF ones.
Providing Let us make it clear by reviewing the formulae already given. In eqn (4) as a result of HF-XCD being approximated by the 1-RDM, the Coulomb contribution is not contemplated. Indeed, if we separate spin cases already in eqn (4), it is pretty clear it must be g ab xc (r 1 ,r 2 ) = 0, since the 1-RDM does not contain cross-spin terms. The same holds for ESI based on HF-XCD, namely, d HF (A,B), d A (A,B) and d B3LYP (A,B), since for the latter exact 2-RDM is not available within the Kohn-Sham formalism, and XCD is also approximated as a HF-XCD. On the other hand, when XCD is calculated exactly from a correlated 2-RDM the Coulomb contribution arises, g ab xc (r 1 ,r 2 ) = g (1)a (r 1 )g (1)b (r 2 ) À g (2)ab (r 1 ,r 2 ), in particular d(A,B) for the correlated calculation contains such a contribution.
A different case is that of the Fulton index. As noticed by Wang and Werstiuk, 64 d F (A,B) can be regarded as an approximation to d (A,B), and its origin comes from the work of Buijse and Baerends 65,66 who proposed the square of eqn (9) as an approximation to the actual XCD for correlated calculations. Hence, unlike HF-XCD we may expect d F (A,B) to partially include the Coulomb information of the 2-RDM, supposing a different and closer 65 approach to the actual XCD.
As a result, when using a HF-XCD (even if we perform a correlated calculation) we are (a) not explicitly considering the Coulomb interaction and (b) approximating parallel spin XCD using eqn (4). In order to measure the first effect we simply must check d ss 0 (A,B), which will be equal to zero for any ESI calculated from HF-XCD. In order to judge the approximation of HF-XCD to recover the actual XCD for correlated calculations, as already mentioned, we must focus on the cumulant term.
But since we want to isolate the effect of Coulomb correlation from the actual approximation given by HF-XCD we have collected the difference d ss (A,B) À l A (A,B) in the last column of Tables 1 and 2; it measures the difference between XCD and HF-XCD for like spins. The reader can check that for systems with merely electrostatic or weak interactions, d ss À d A can be more important that the unlike spins pairs corrections, d ss 0 (A,B), (as is the case of LiF, He 2 , and Ar 2 or LiH).
In summary, as expected, the inclusion of electron correlation in the calculation of ESI plays a marginal role in non-bonded or weakly bonded species, has a slight influence on ionic species, and a prominent role in equally-shared or polar covalent bonding. Such is the case for the isoelectronic series N 2 , NO + , CN À , and CO, where the inclusion of correlation by means of CISD drastically changes the value of the ESI.

Desirable properties of an ESI
Providing a universal unique definition of an ESI does not (and maybe cannot) exist, one can only use the knowledge gathered on chemical grounds to impose or try to define an ESI so that a certain set of desirable properties are fulfilled. For instance, reproducing the growing tendency with lower saturation, as in the series: ethane, ethene, ethyne; an intermediate ESI value for C-C in benzene between formal single and double bonds, asymptotically tending to zero as the distance between the two atoms increase to infinity, are among the desirable properties one may impose on an ESI in order to give truthful results. For instance, some of us 67 recently discussed the dissociation of the H 2 molecule with QTAIM-ESI.
Another desirable property of an ESI is always being non-negative, and to tend towards zero for a no electron sharing situation. If the index can achieve negative values, probably the index can hardly be called a sharing index, since a negative sharing is rather unphysical. One should take into account that this is the case of Wiberg-Mayer Bond Order, because it is based on Mulliken's partition, whereas, as already noticed, ESI from Fuzzy-Atom and QTAIM is always positive.
A recent study by Ponec 68 has brought our attention to another interesting property for an ESI to be fulfilled. This is the case of the dissociation of ionic molecules into neutral species. 69 Such molecules should have a small number of electrons shared in their equilibrium geometry, however, when going to larger atomic distances one of the electrons should go from the negative charged atom to the positive one. Therefore, at some given distance, there must be a maximum of electron sharing, since we start from a rather electrostatic interaction and move towards a more covalent situation. Ponec showed that indeed, QTAIM-ESI-or simply the delocalization index (DI) 20,21 exhibits such a maximum, whereas a Mulliken-like scheme clearly fails to reproduce it.
Since Ponec did not study this dissociation process with the Fuzzy-Atom approach, we decided to recalculate the d F (A,B) and d A (A,B) indexes with both QTAIM and our Fuzzy-Atom partition for the LiH dissociation by using the same method, CAS(4,4) with a rather large basis set, 6-311G**. Results are reported in Fig. 1. One can easily see how, as already reported by Ponec, QTAIM ESI detects a maximum for both ESI, whereas the simplest Fuzzy-Atom approach based on Becke's partitioning fails to reproduce this maximum. The fact that the ratio between the atomic radius of the Li and H atom is kept fixed during the dissociation is clearly responsible for this situation. Hence, we have decided to dynamically tune the radii used in the Fuzzy-Atom approach 11 to reproduce the maximum of electron sharing which in principle one would predict for the dissociation of this ionic species. Because of the success of the QTAIM approach, based on the topology of the electron density, a natural way to define the radii of an atom in a molecule may be its distance to the bond critical point at which it is connected. Whereas for polyatomic molecules this definition can be controversial (a given nucleus can be connected through different gradient paths to several critical points), such a definition for diatomic molecules is quite simple. The results obtained for such a Fuzzy-Atom approach are striking, since a noticeable improvement in the results is gathered (cf. Fig. 1). With these radii the Fuzzy-Atom predicts a slight maximum near the region where QTAIM-ESI displays it.
This finding encourages the possibility of using the information from the topology of the density as a way to dynamically adjust the ratio of atomic radii for a Fuzzy-Atom calculation. Needless to say, these atomic basins, unlike QTAIM's, would not possess the important features related with the fulfilment of the atomic virial theorem and other interesting properties. However, the numerical values of the ESI will be quantitatively very similar to the true QTAIM ones. This gives the possibility of using the Fuzzy-Atom approach when the calculation of a QTAIM ESI is prohibitive from a computational point of view, or even to avoid the presence of some spurious non-nuclear attractors. 70,71 Research in this line is underway in our laboratory.

Aromaticity indexes at correlated level. PDI and FLU
So far, the FLU and PDI aromaticity indexes have been reported for monodeterminantal wavefunctions (either HF or Kohn-Sham), but no attempt to explore these indexes for many-body methods has yet been carried out. Since the calculation of 2-RDM is prohibitive for these systems, here we propose to use the d F (A,B) and d A (A,B) that need only the 1-RDM to construct these aromaticity indexes. The computational cost of both ESI is the same, and thus one is tempted to consider Fulton's ESI which is the closest one to the fully correlated ESI, d(A,B) (vide supra).
Let us notice an important feature of the results for benzene reported in Tables 1(b) and 2(b). Contrary to what is found with the QTAIM partition, regardless of the ESI used, the Fuzzy-Atom approach shows that the electron sharing is smaller for para-related atoms than for meta-related atoms in benzene. Therefore, one would anticipate Fuzzy-Atom PDI to fail to reproduce this aromaticity trend found by QTAIM. However, as already demonstrated, 29 Fuzzy-Atom does not suffer this drawback. In addition it is worth commenting that using a stiffness value of k = 4 (by default we use k = 3) we reproduce a larger para-related sharing than the metarelated one, and for k = 6 the results mimic the QTAIM ones. The results reported in Table 3 show that in general, for both partitions QTAIM and Fuzzy-Atom, the inclusion of the electron correlation generally leads to smaller values of FLU and PDI, as compared to the HF values. It is specially well known for Fulton's ESI, whose PDI values are sensibly smaller than HF ones, whereas Á ngyan's are somewhat closer to the HF values. As compared to previous work where the effect of the basis set was discussed, 29 the effects due to the inclusion of electron correlation are more well known. Nonetheless, the general trends described hold, so that these electronic aromaticity measures may be regularly applied to correlated calculations.

Conclusions
The effect of inclusion of the electron correlation has been analyzed for the ESI from different partition schemes (QTAIM and Fuzzy-Atom), based on both the first-and the second-order density matrices. It has been shown that for weakly bonded molecules most indexes coincide and the effect of the electron correlation is minor; whereas for covalent and polar covalent molecules the electron correlation effects are well known, and the ESI based on the 1-RDM do not wholly reproduce the trends shown by the ESI based on the 2-RDM.
The Fuzzy-Atom approach for a fixed set of radii (we have tested Koga and Slater ones) has been shown not to reproduce a maximum of electron sharing in the dissociation of LiH species in neutral atoms and to yield larger meta-than pararelated sharing in benzene. Nevertheless, an improvement in the description is gathered if the relative radius of both atoms is dynamically adjusted according to the distance of the atoms to the bond critical point connecting them.  Finally, aromaticity indexes have been studied for the first time for correlated calculations, leading to a general trend similar to that already reported for monodeterminantal methods. Nevertheless, it is worth mentioning that these indexes depend more on the methodology employed than in the basis set used. PDI and FLU for correlated calculations can be calculated yielding reasonable results.