Electron localization function at the correlated level

The electron localization function ELF has been proven so far a valuable tool to determine the location of electron pairs. Because of that, the ELF has been widely used to understand the nature of the chemical bonding and to discuss the mechanism of chemical reactions. Up to now, most applications of the ELF have been performed with monodeterminantal methods and only few attempts to calculate this function for correlated wave functions have been carried out. Here, a formulation of ELF valid for monoand multiconfigurational wave functions is given and compared with previous recently reported approaches. The method described does not require the use of the homogeneous electron gas to define the ELF, at variance with the ELF definition given by Becke. The effect of the electron correlation in the ELF, introduced by means of configuration interaction with singles and doubles calculations, is discussed in the light of the results derived from a set of atomic and molecular systems. © 2006 American Institute of Physics. DOI: 10.1063/1.2210473


I. INTRODUCTION
The first to recognize the importance of the electron pair in chemical bonding was Lewis 90 years ago. 1 His idea that chemical bonds are formed by shared electron pairs has been one of the most fruitful concepts in modern chemistry. 2,3The electron pair plays a dominant role in the present day chemistry, being a powerful concept to rationalize and understand the molecular structure ͓through the valence shell electron pair repulsion ͑VSEPR͒ theory 4 ͔, bonding, and chemical reactivity.Because of their importance, quantum chemists have pursued the development of procedures to determine the spatial domains where local electron pairing takes place in molecules.The first attempts in this direction were based on the idea of the so-called localized orbitals. 5The Laplacian of the one-electron density is another function used to denote electronic localization. 6,7Afterwards, with the intention to calculate the electron pair regions with greater precision in the molecular space, theoretical chemists turned their interest to two-electron quantities. 8The simplest of such quantities that describes the pair behavior is the so-called pair density.Other related functions employed to locate the electron pairs have their origin in the seminal works of Artman 9 and Lennard-Jones. 10 Among them we can mention the exchange-correlation density, 11,12 the Fermi hole, 13 the Lennard-Jones function, 14 the Luken-Culberson function, 15 and the conditional pair density. 15,168][19][20] The definition of the conditional probability ͑CP͒, i.e., the probability density of finding electron 2 nearby, r 2 , when electron 1 is at r 1 , reads with ␥ ͑2͒ ͑r 1 , r 2 ͒ and ͑r 1 ͒ being the pair density and the density, respectively.Both the pair density and the CP contain all necessary information about the correlation of electrons.However, a good advantage of the CP function is that it is actually discounting irrelevant information concerning the position of the reference electron.Hence, it is not surprising that Becke and Edgecombe 17 chose the CP of same spin electrons as a measure of electron localization.In particular, they used the spherical average of the CP expanded by Taylor series around the position of the reference electron: 21 ͗e s•١ P ͑r,r + s͒͘ e sٌ dd ͪ ͉P ͑r,r + s͉͒ s=0 = sinh͑s ٌ ͒ sٌ ͉P ͑r,r + s͉͒ s=0 , ͑2͒ which after using the Taylor expansion of sinh around the reference electron ͑s =0͒ generates an expression whose leading term reads ͗e s•١ P ͑r,r + s͒͘ = ͑1 + 1 6 s 2 ٌ s 2 + ... ͉͒P ͑r,r + s͉͒ s=0 Ϸ 1 6 ͉s 2 ٌ s 2 P ͑r,r + s͉͒ s=0 .͑3͒ The right-hand side ͑rhs͒ of Eq. ͑3͒ holds because the Pauli principle forces the same spin pair density at the coalescence point to be zero, ␥ ͑r , r͒ = 0. Thus, the leading term of the Taylor expansion of the spherically averaged CP is proportional to Hence, Becke and Edgecombe 17 used the relative ratio of D with respect to the same quantity for the homogeneous electron gas ͑HEG͒, assuming a monodeterminantal situation where the pair density can be calculated from the density: 22 with c F = ͑3/10͒͑3 2 ͒ 2/3 being the Fermi constant. 23Hence the relative ratio reads The ͑conditional͒ probability of finding an electron with spin when there is already another electron with the same spin nearby is lower when the former is localized.Therefore, the ratio in Eq. ͑6͒ accounts for electron localization; the higher it is, the lower the localization.In order to range the electron localization measure in the interval ͓0,1͔, Becke and Edgecombe chose the following scaling to define the ELF: 17 thus, ELF= 1 corresponds to a completely localized situation, and 0 to its opposite.ELF= 0.5 within a HEG.The ELF helps in clarifying the structure of molecules 24 and the shell structure of atoms, 25,26 and its topological interpretation gives rise to the molecular space splitting in terms of chemical meaningful regions such as bonds, lone pairs, and core regions. 19The ELF has also been calculated separately for alpha and beta orbitals with the aim to study different spin contributions on some radical systems. 27The ELF has been used to gain insight on aromaticity, 20,28 to analyze the changes in electronic structure along a given reaction coordinate, 29 and to distinguish between similar mechanisms of reaction. 30ecke and Edgecombe 17 developed the ELF formula in Eq. ͑7͒ for a monodeterminantal case, where the two-body quantities can be expressed in terms of one-body ones.The fact that Becke and Edgecombe based their ELF on the CP ͑a two-body quantity͒ and later on worked out their expression for the monodeterminantal case had the disadvantage of the loss of generalization, since only monodeterminantal cases were contemplated.However, they overcame an epistemological problem: the natural extension of the ELF at correlated level should also take into account the correlated expression of its reference, the spin-unpolarized HEG.To the best of our knowledge, the exact analytical expression of the pair density for the HEG is unknown.
The investigations of the ELF at correlated level are scarce, mainly due to the considerable computational effort needed for the calculation of the pair density.The purpose of this paper is to investigate the ELF at correlated level of theory.Thus, we firstly explore the definitions of the ELF already given for correlated calculations to show they are indeed equivalent.Secondly, we apply the ELF to the calculation of some atomic systems up to krypton, to analyze its ability to recover shell structure and shell populations at the configuration interaction singles and doubles ͑CISD͒ method.And finally, we take a series of isoelectronic molecules ͑CO, CN − , N 2 , and NO + ͒ and another series of molecules known to be problematic for monodeterminantal approaches ͑F 2 , FCl − , and H 2 O 2 ͒ to show the effect of the inclusion of electron correlation in the population and variance of ELF basins.

II. METHODS
This section is organized as follows.First, we analyze the previous definitions of correlated ELF given by different authors and we prove that all of them are equivalent; second, we provide the definitions of different populations and their variances; and finally, we discuss the computational tools used in this study.

An alternative ELF with no need for arbitrary references to HEG
The ELF thus defined, Eqs.͑6͒ and ͑7͒, has advocated several criticisms because of the arbitrariness of choosing the HEG as a reference.However, recently one of us 31 derived a quantity to account for the electron localization that, in turn, was identical to the ELF, and thus served as an alternative derivation to ELF without the arbitrariness of measuring with respect to the HEG.A similar process was followed by Kohout 32,33 to achieve identical results.Here we demonstrate the equivalence of both formulations, which lead to the same ELF definition ͑vide supra͒.
The interpretation of the ELF in terms of the second derivative of the pair density is given in the papers of Dobson. 34The ELF is indeed the angularly averaged ͑spherical averaged͒ measure of the Fermi hole curvature, and in turn, the Fermi hole curvature can be also understood as the density of the kinetic energy of relative motion pairs of the same spin electrons centered at r. Dobson suggests the integration of the CP ͓Eq.͑3͔͒ over a small sphere of radius R: ͉s 2 ٌ s 2 P ͑r,r + s͉͒ s=0 4s 2 ds Ϸ 2R 5 15 ٌ͉ s 2 P ͑r,r + s͉͒ s=0 .͑8͒ When R is close to zero the charge inside the sphere can be approximated as the density of its central point multiplied by its volume, and therefore the density is proportional to the inverse of the third power of R. In this case we are left with where we have supposed a closed-shell situation for the sake of simplicity.The reader should note that, except for constant factors, we have reached the expression already given by Eq. ͑6͒ in terms of two-body quantities without using the arbitrary reference of the HEG, or in other words, one can give another interpretation of Eq. ͑6͒ used to define the ELF.
Hence the ELF can be understood as the number of electrons with spin near r when there is another electron with the same spin by r.It is worth noticing that the ELF is actually an approximation to this quantity.
The expression in the numerator of Eq. ͑9͒ is a measure of the spin composition of the system 31 and is usually referred to as the Fermi hole curvature.It is not the real Fermi hole curvature, since the ␥ ␣␣ is not a Fermi hole, but a pair density of the same spin electrons.However, one can easily demonstrate that the curvature of the Fermi hole, is indeed very similar to the pair density curvature.Hence The rhs of the latter equation holds for volumes where the density is small enough to be considered constant.Here, the advantage of using the curvature instead of the function or its gradient clearly shows; the former is the only one that does not vanish when r 1 = r 2 . 35Therefore, one could get an alternative ELF definition in terms of the Fermi hole curvature, which would be equal to the current definition except for the exponent of the density.

One of us recently 31 proved that
D͑r;q i ͒ Ϸ N ¯2/3 ͑r;q i ͒c ͑r͒, ͑12͒ where D is the size-dependent pair composition of the system ͓our Eq. ͑6͒ considered a small sphere centered at r with charge q i ͔, N ¯the population of the volume enclosed in the sphere centered at r, and c the size-independent spin composition given by ٌ 2 2 P ͑r 1 , r 2 ͒ /3 5/3 ͑r 1 ͒.This way, an empirical proof for the approximation from Eq. ͑8͒ to the current definition of the ELF was given.In Ref. 31 it is shown how the approximation works well for small volumes, differing slightly for large volumes where no electron homogeneity is granted; in Eq. 30 of Ref. 31 a generalization of ELF for open-shell cases is provided as follows: where we have taken proper normalization in order to reproduce our Eq.͑6͒ for closed-shell species.On his side, Kohout 33 decided to split the molecular space into regions of fixed charge as follows: 11 N ¯͑r;⍀ q,i ͒ = q = ͵ ⍀ q,i ͑r͒dr.͑14͒ Taking the same spin pair density as given by Eq. ͑10͒, we have In the situation of fixed charge ͓Eq.͑14͔͒ we get So, in this context, Eq. ͑16͒ is measuring not only the number of expected pairs of electrons in ⍀, but provided the population in these regions is constant, we are indeed measuring F, which is to say-as previously demonstrated by Fradera et al. 3 and Bader 7 -a quantity that accounts for the electron localization.This q -restricted pair population is what Kohout suggests ͑after proper scaling͒ as a measure of electron localizability.In fact, it is not the quantity itself, but the first nonvanishing term in the expansion of Eq. ͑16͒ around position r 1 : Expression ͑17͒ is indeed equal to the latter expression ͑8͒ except for a factor proportional to the density, which, in turn, is constant in this context.Then, Kohout uses an arbitrary factor to normalize the latter expression and fit the one already given by Becke and Edgecombe. 17It is worth noting that Kohout et al. did some preliminary complete active space self-consistent field ͑CASSCF͒ calculations using Eq.͑17͒ as an expression 32 for the correlated ELF.Moreover, Kohout et al. have recently proposed an antiparallel version of the localization function. 36Finally, the paper of Scemama et al. 37

B. Variance and covariance basin populations
The partition of the molecular space enables basinrelated properties to be calculated by integrating the density of a certain property over the volume of the basins.Thus, for a basin labeled ⍀ i , one can define the average population as and its variance One can also define likewise the pair population:

͑20͒
The variance is a measure of the quantum mechanical uncertainty of the basin population, which can be interpreted as a consequence of the electron delocalization.The variance 2 ͑⍀ i ͒ can also be split in terms of contributions from other basins, where V͑⍀ i , ⍀ j ͒ is the covariance and is calculated as follows:

͑22͒
Finally, the integrated spin density may also be important for open-shell systems:

C. Computational details
Since coupled cluster and MP2 methods do not fulfill the Hellmann-Feynman theorem, neither the density nor the pair density can be obtained as an expectation value.Therefore, the pair density and the density thus obtained are relaxed and, in general, not N representable 38 ͑i.e., derivable from an N-electron wave function͒.For this reason, here we have preferred to deal with the CISD unrelaxed density and pair density that have been employed in Eq. ͑13͒ to obtained correlated ELFs.The complete active space ͑CAS͒ unrelaxed density and pair density could also be used, but in this work we will focus our study on CISD results.The density functional theory ͑DFT͒ calculations in this paper are performed with Kohn-Sham formalism, and building up an approximate HF-like pair density, 39 thus reducing the ELF expression to that originally given by Becke and Edgecombe. 17n order to compute the ELF from Eqs. ͑13͒ and ͑7͒, one must calculate the curvature of the pair density, that in terms of molecular orbitals reads In case we are dealing with a linear combination of Gaussian primitives, the latter expression turns into derivatives of Gaussian functions, whose analytical expressions are straightforward.In this work, we have worked with such functions, obtained with the GAUSSIAN03 program. 40From the same program the coefficients of the expansion were obtained and used by our own set of programs to construct the ⌫ matrix in Eq. ͑24͒.The wave function together with the ⌫ matrix were read by our own modification of the TOPMOD package of programs 41 to calculate the curvature of the pair density and finally, the ELF.
In this study we have focused on configuration interaction ͑CI͒ calculations to introduce the electron correlation.The electron correlation of parallel spins was introduced by mixing the Hartree-Fock ͑HF͒ determinant with other monoexcited and biexcited determinants, which is known as a configuration interaction with singles and doubles ͑CISD͒. 42nlike atoms, CISD calculations for molecules were performed within the restriction of frozen core.All molecules were optimized with the CISD, HF, B3LYP, and BP86 levels of theory with Pople's 6-311+ G͑2df ,2p͒ basis set.Additional single point calculations were performed for all methods at the HF geometry, with the aim to study the purely electronic effects in the inclusion of electron correlation.
The coefficients of the CI expansion have been obtained from energy calculations performed with the GAUSSIAN03 program; 40 only coefficients with values above 10 −6 were taken.The density and the pair density thus obtained are unrelaxed.It is well known that relaxed CISD densities are better for the calculation of one-electron response properties, 43 but this is not the goal of the present study, and since they are not N representable, 38 the unrelaxed ones have been preferred for this study.Only values of density and pair density above 10 −12 a.u.have been considered.

III. RESULTS AND DISCUSSION
With the goal to briefly explore the influence of the electron correlation in the ELF expression, we have limited ourselves to calculate some atomic systems ͑atoms from Li to Kr͒, four diatomic isoelectronic molecules ͑CO, N 2 , NO + , and CN − ͒, and a series of molecules which are expected to be more affected by the inclusion of the electron correlation ͑F 2 , FCl − , and H 2 O 2 ͒.

A. Atomic calculations
The ELF has been previously demonstrated as a good tool to discern the shell structure of atoms. 25,26,44Unlike the

024301-4
Laplacian of the density, which is unable to recover quantitatively the number of electrons per atomic shell 45 and even the shell structure for heavy atoms, 46 the ELF reproduces the shell structure and gives a quite good approach ͑0.1-0.2 electrons of deviation 25 ͒ to the number of electrons per shell.Recently Kohout demonstrated that the one-electron potential also reveals with a reasonable accuracy the electron numbers for the occupation of the shells. 47n order to check whether ELF provides better shell numbers with the improvement of the wave function, 25 ELF calculations for the ground state of atoms from Li to Kr with CISD/ 6-311+ G͑2df ,2p͒ level of theory have been performed.The basin boundaries have been assumed to be spherical and the integration of the atomic electron density has been carried out analytically ͑see Appendix͒.Table I contains the data from ELF calculated with CISD and the 6 -311+ G͑2df ,2p͒ basis set.The radius of the K shell is the result of the competition between the attractive ͑electron-nucleus͒ and repulsive forces ͑electron-electron͒.Since the attractive potential increases with Z, the net effect is the reduction of the shell radii, and the increase of the shell number with Z.For the same reason, for atoms from Li to Sc, L shell populations lie below the expected on chemical grounds, decreasing shell population with Z.Because of the screening effect of d electrons on the L shell, the opposite situation is found for atoms from Ti to Kr.
Shell radii and population of atoms with electrons in 3d orbitals are quite unpredictable, especially concerning the M and N shell numbers; their configuration is no longer ruled exclusively by the Aufbau procedure and Hund's rule plays a prominent role.One should also take into account that real orbitals p and d do not guarantee the spherical symmetry except for empty, half-occupied and full occupied configurations.

B. The isoelectronic series: CO, CN − , N 2 , and NO +
The four molecular systems CO, CN − , N 2 , and NO + are isoelectronic, but present different bonding situations.Although the four systems are covalent, one would advance nitrogen molecule to have nonpolar covalent character, while NO + , CN − , and CO are somewhat polarized.In the formal Lewis sense one would expect a triple bond for all species and some lone pairs around the nuclei.This kind of mol- ecules with large electron charge in the bonding regions are especially susceptible to the inclusion of electron correlation in the calculations. 48able II contains the basin population and its variance for the isoelectronic series, with the four methods studied at the optimized geometry of each method and also at the HF geometry.Focusing either in the results derived from the calculations at different geometries ͑which gather simultaneously the electronic and geometric effects of electron cor-relation͒ or in the results with HF geometry, we can generally assess that the population of monosynaptic valence basins increases at the expense of populations in disynaptic bonding basins, which is consistent with the general expectation of smaller force constants and larger internuclear distances at the CI or DFT levels of theory.At the optimal geometry of each method, the calculated decrease of the V͑A , B͒ basin population with respect to the HF value is larger for the BP86 functional than for the B3LYP and CISD approaches, the latter yielding the smallest change.The net charge transfers towards the lone pair due to electron correlation are often correlated with the bond lengthening: this is always the case for the BP86 and B3LYP results, and most of CISD results, whereas CISD value for V͑N͒ of CN − shows an unexpected opposite behavior.
The loge theory, as introduced by Daudel et al., 49 suggests the minimization of electron fluctuation in individual domains as the way to reach the most likely decomposition of molecular space: the smaller the lack of localization, the better the decomposition is.It suggests the idea of taking the sum of variances ͑⌺ 2 in Table I͒ of basin populations as a measure of the best partitioning achieved. 50In this context, Gallegos et al. 51 explicitly calculated the maximal probability domains for N 2 and CO.Since the best partition found by Gallegos et al. for N 2 and CO yields smaller populations for the bonding regions than the ELF partition at the HF level, one would expect the populations of bonding regions for N 2 and CO to decrease with respect to the values given by HF with the inclusion of electron correlation.This is indeed what we have found in our DFT and CISD calculations.Furthermore, in all cases the CISD partitioning yields the smallest sum of variances.It is worth noting that for DFT the variance calculation is not exact since an approximate expression of TABLE II.Basin populations and its variance for CO, CN − , N 2 , and NO + calculated at the HF, B3LYP, BP886, and CISD levels.The ͑HF͒ columns correspond to single point calculations at a given method with the HF geometry.Bond lengths in Å and populations and variances in electrons.Since core basins do not significantly change from one method to the other, only valence basins ͑V͒ are included.the two-particle density is used instead of the ͑unknown͒ exact one.These results seem to indicate that further correlated ELF partitioning can lead to better agreement with loge theory.However, one should take into account that only a real space exploration of maximal probability spaces against ELF partitioning will lead to the definite answer about the similarity of these two approaches.

C. The series F 2 , FCl − , and H 2 O 2
Finally, the results for F 2 , FCl − , and H 2 O 2 are given in Table III.These molecules are supposed to be more affected by the inclusion of electron correlation. 52Indeed, F 2 exhibits two symmetric V͑F,F͒ basins for correlated methods, instead of a single one as in the HF case, holding less electrons than in the single V͑F,F͒ situation; this number is specially small for DFT methods.H 2 O 2 is partially well described by the HF method, yielding lower V͑O,O͒ populations and variances for CISD and DFT calculations; for the latter this basin is even split into two symmetric ones.
The FCl − is a special case, where the inclusion of electron correlation leads to notorious change in both the structure and the electron distribution.First, one must notice that HF equilibrium distance is about 1 Å larger than that obtained with more accurate methods; DFT methods approach better this distance, while our CISD result agrees with that of 2.111 Å given by the EOMEA-CCSD/aug-cc-PV5Z method. 53Specifically, it is worthy to comment that the spin density integrated in each basin, cf.Eq. ͑23͒, is zero for Cl basin within HF.This is an indication that the HF picture of the molecule is a closed-shell fragment given by the anion Cl − interacting with a neutral fluorine atom.Therefore, the ELF results reported for the HF wave function differ considerably from those calculated with DFT or CISD methods.In addition, are the differences of the populations and variances between CISD and DFT noticeable.The inclusion of the electron correlation leads to not only shorter bond distances, but also to partial electron transfer from the chlorine atom towards the fluorine one; an electron transfer slightly exaggerated by the DFT methods, as one can deduce from both the populations and the spin integrated densities.

IV. CONCLUSIONS
In this paper we review the ELF definition for correlated wave functions, showing the equivalence of the most recently proposed definitions of ELF.This ELF definition reproduces that of Becke and Edgecombe for monodeterminantal wave functions.We give the first results for CISD calculations of ELF, together with the first populations and variance analysis for a correlated partitioning according to the ELF function.Our analysis shows a clear reduction of the bonding region populations as expected from previous results derived from the loge theory.It is worth noting that there is a qualitative agreement between the CISD and DFT results in all the systems investigated here.Although more molecular systems need to be explored to assess the effect of the electron correlation in the ELF and its consequences for the description of the molecular structure, results for F 2 , ClF, and HOOH indicate that correlated ELF can differ dramatically in some cases from the uncorrelated one.Since the description of some chemical reactions may also be very sensitive to the inclusion of electron correlation, we are currently investigating the ability of the correlated ELF to deal with such chemical processes.We hope to be able to report the results on chemical reactions analyzed using correlated ELF in due time.

ACKNOWLEDGMENTS
Financial help has been furnished by the Spanish MEC Project Nos.CTQ2005-08797-C02-01/BQU and CTQ2005-02698/BQU.The authors are also indebted to the Department d'Universitats, Recerca i Societat de la Informació ͑DURSI͒ of the Generalitat de Catalunya for financial support through the Project No. 2005SGR-00238.One of the authors ͑E.M.͒ thanks the Secretaría de Estado de Educación y Universidades of the MECD for the doctoral fellowship No. AP2002-0581.The authors also thank the Centre de Supercomputació de Catalunya ͑CESCA͒ for partial funding of computer time.

APPENDIX: INTEGRATION OF ATOMIC POPULATIONS
In order to compute the populations given in Table I one must find the ELF minima along the radial direction, since they indicate the separation between shells.Afterwards, one must integrate the electron density over a spherical basin of radius r a , say, to compute integrals of the form

TABLE I .
Shell radii ͑r X ͒ and shell numbers populations ͑N X ͒ for the ground state of atoms from Li to Kr.All quantities in a.u.

TABLE III .
Basin populations and its variance for H 2 O 2 , F 2 , and FCl − calculated at the HF, B3LYP, BP86, and CISD levels within the corresponding geometry.Distances in Å, angles in degrees, and populations and variances in electrons.Since core basins do not significantly change from one method to the other, only valence basins ͑V͒ are included.OH = 0.950 Å, r OO = 1.475Å, and ЄHOO = 94.8 ͑Ref.54͒.
a Experimental values are r b Experimental bond length is 1.417 Å ͑Ref.54͒.