The hardness kernel as the basis for global and local reactivity indices

In a very recent article (Torrent‐Sucarrat et al., J Comput Chem 2007, 28, 574), we have shown how to evaluate the global hardness for polyatomic molecules using a hardness kernel approximation. We present here an extension of this previous work by improving the model used to evaluate the hardness kernel and the Fukui function. In addition, the concept of the local hardness is analyzed in detail, and for the first time, profiles of local hardness with kinetic and exchange‐correlation contributions for polyatomic molecules are reported. Finally, the concept of condensed atomic hardness is introduced and its usefulness as chemical reactivity descriptor is examined. © 2007 Wiley Periodicals, Inc. J Comput Chem, 2008


Introduction
In the last decade, conceptual density functional theory (CDFT) [1][2][3][4][5] has become one of the most popular and relevant branches of density functional theory (DFT). 1,6 CDFT is focused on the development and application of new reactivity descriptors and principles based on DFT concepts, bringing to the scientific community useful tools to analyze, understand, and predict the chemical reactivity. [1][2][3] Reactivity indices in the CDFT framework can be split into three general groups: global, local, and nonlocal. The overall chemical characterization of a molecule is accomplished by the global reactivity indices (e.g., chemical potential and global hardness and softness). Local descriptors like the Fukui function, local hardness, and softness are different in several positionsr of the space, yielding information on the preferential site for a chemical reaction. Finally, the nonlocal descriptors (e.g., hardness and softness kernels and the linear response function) are functions of two positionsr andr 0 in the space. The nonlocal character of these descriptors thus far prevented a simple direct link with chemical reactivity so that the major part of the research devoted to CDFT has been focused on the global and local reactivity indices. [1][2][3] The nonlocal reactivity descriptors have the interesting feature that the local and global indices can be generated from their expressions. [7][8][9][10][11][12][13][14] For instance, the local hardness, gðr Þ, can be obtained by integrating the nonlocal hardness kernel, gðr;r 0 Þ, multiplied by the Fukui function, f ðr 0 Þ. In addition, the global hardness, g, arises from the integration of gðr Þ multiplied again by the Fukui function. The drawback of these methodologies is that analytical expressions for the nonlocal reactivity descriptors must be derived. Fortunately, contrary to other nonlocal descriptors, like the softness kernel and the linear response function, analytical expression for gðr;r 0 Þ can easily be obtained through the second functional derivative of the Hohenberg-Kohn universal density functional, F½qðr Þ, with respect to qðr Þ and qðr 0 Þ.
There have been few attempts in the past to carry out the integration of the gðr;r 0 Þ; all of them limited to atoms [15][16][17][18][19][20][21] or with highly simplified models for the hardness kernel. 22,23 In a very recent article, 24 the present authors have evaluated for the first time the global hardness for polyatomic molecules by inte-gration of gðr;r 0 Þ using the density of the HOMO orbital, q HOMO ðr Þ, as an approximation for the Fukui function and the Thomas-Fermi and Dirac models as approximations of the kinetic and exchange contributions to the hardness kernel, respectively (as well as the leading Coulombic term).
The present work focuses on three main goals. Firstly, the hardness kernel of our previous work is improved by introducing correlation contributions via a Wigner-type local functional and correcting the Thomas-Fermi kinetic energy functional with 1/9th of the Weizsäcker functional. Secondly, profiles of local hardness with kinetic and exchange-correlation contributions for polyatomic molecules are reported, to the best of our knowledge, for the first time. Finally, the concept of condensed atomic hardness is analyzed in detail and its evaluation from the numerical integration of the hardness kernel is presented paying special attention to its potential use as chemical reactivity descriptor.
This article is organized as follows. Theoretical Background and Computational Details section contains the definitions of several CDFT descriptors computed in this work together with the computational methodology used. Then, in Results and Discussion section for selected species, we present profiles of local hardness and we analyze the results obtained for global and condensed atomic hardnesses, and finally, our conclusions are summarized in Conclusions section.

Theoretical Background and Computational Details
The hardness kernel is defined as 7,8 F½qðr Þ is a functional that contains the classical electron-electron Coulomb repulsion functional, J½qðr Þ, the kinetic energy density functional, T½qðr Þ, the exchange density functional, E X ½qðr Þ, and the correlation density functional, E C ½qðr Þ (note that the present partitioning of F½qðr Þ does not follow the Kohn-Sham scheme and that the exchange and correlation functionals do not contain a kinetic contribution). Then, one can write where the only part that can be written exactly is Explicit exact expressions are unknown for the other terms in F½qðr Þ, although approximate expressions exist, functionals of qðr Þ and its derivatives, which can be inserted in eq. (1). It is important to remark that a hardness kernel in the Kohn-Sham context has been recently derived by Fuentealba et al. 13 In this work, the kinetic energy functional will be approximated using the Thomas-Fermi functional, T 0 ½qðr Þ, 25,26 with C T 5 (3/10)(3p 2 ) 2/3 , corrected by adding 1/9th of the Weizsäcker functional, T w ½qðr Þ, 27 In addition, the Dirac, E X ½qðr Þ, 28 and Wigner-type local correlation, E C ½qðr Þ, 29 functionals will be used as approximations of the exchange and correlation functionals, respectively. and where C X 5 (3/4)(3/p) 1/3 , B 5 1/21.437, and A 5 9.81/21.437. Using this approximate model of F½qðr Þ in eq. (1), one gets the following expression for gðr;r 0 Þ 9,14 g Àr where the first square bracket contains the hardness kernel contributions from the Thomas-Fermi, Dirac, and Wigner functionals, while the second square bracket is the contribution of the 1/9th of the Weizsäcker functional (a detailed functional derivation of the 1/9th Weizsäcker functional term can be found in ref. 14). Finally, the remaining term of eq. (8) is the Coulomb contribution. As mentioned before, the local hardness function can be evaluated from the hardness kernel as Combining this model of gðr;r 0 Þ with eq. (9), one obtains an analytical expression for gðr Þ.
*The original definitions of the global, local, and nonlocal hardnesses contain an arbitrary factor of 1/2.
In this context, it is worth noting that, according to Harbola et al., 12 the local hardness becomes a constant (and therefore equal to the global hardness) when the exact ground state density is considered, just as the chemical potential equalizes in the ground state. 30 This conclusion is also recovered when a variational principle for determining the Fukui function and the global hardness is applied to an electronic system, as put forward by Ayers and Parr, and Chattaraj et al. 4,9 In the Results and Discussion section, we will show that for practical calculations the local hardness function is far from being constant in all points of space. To connect the global and local hardnesses, it is necessary to use the following expression Then, the analytical expression for the global hardness using the model of gðr;r 0 Þ [see eq. (8)] becomes As discussed in our previous work, 24 the density of the HOMO (highest occupied molecular orbital) can be used as approximation of the Fukui function. Other approximations of the Fukui function like q LUMO ðr Þ and 1=2ðq HOMO ðr Þþq LUMO ðr ÞÞ can be also used, although spurious values have been obtained for some systems. 24 In particular, when using large basis sets, especially including diffuse functions, it is possible to obtain a spurious Rydberg-type LUMO, which leads to unphysical, even negative, values of the global hardness when calculated from eq. (12). Equations (10) and (12) can be formally written in a more compact way as and where the expression between brackets represents the approximation used to evaluate the Fukui function. The Eqs. (10) and (12) illustrate the use of the hardness kernel as the basis for global and local reactivity indices. Given a model for the Hohenberg-Kohn universal functional and for the Fukui function, the local and global hardnesses can be evaluated using the information of qðr Þ and its derivatives ðrqðr Þ; r 2 qðr Þ; . . .Þ. The only computational requirement is the calculation of the numerical integrals in eqs. (10) and (12). As we have done in our previous work, 24 the corresponding 3D and 6D numerical integrations are evaluated using the Becke's multicenter numerical integration scheme, 31 as implemented in a program developed in our laboratory.
The key point of this method is the definition of a weight function, w i ðr Þ, assigned to each nucleus i in every point of the space, that allows to express any function, Fðr Þ, as the sum of different atomic contributions Then the numerical integration of Fðr Þ can be reduced to a sum of single-center integrations, which can be carried out using standard single center numerical techniques in spherical polar coordinates The atomic weights w i ðr Þ used in this work are derived from the fuzzy Voronoi polyhedra proposed by Becke, 31 tuned by the Bragg-Slater set of atomic radii 32 and following Becke's suggestion to increase the radius of hydrogen to 0.35 Å . In the present work, each atom has been integrated using Chebyshev's integration for the radial part (40 points) and Levedev's quadrature 33 (146 points) for the angular part. The routine for the Levedev quadrature has been downloaded from ref. 34. This level of methodology allows achieving accuracies of the order of 10 25 and 10 23 a.u. for the 3D and 6D integrations, respectively. The local and global hardnesses of eqs. (10) and (12) have been evaluated for a large set of molecular and atomic systems at the B3LYP level 35,36 and 6-31G, 6-311G(d), and 6-31111G(2d,2p) basis sets, 37,38 using the wavefunctions and densities obtained from the Gaussian 03 package. 39

Quality of the Hardness Kernel
In Table 1, we report the global hardness obtained with eq. (12) and the two most popular approximations for the global hardness* and where I and A are the first vertical ionization potential and electron affinity of the neutral molecule, respectively, and e LUMO and e HOMO are the energies of the lowest unoccupied and the highest occupied molecular orbitals, respectively. The g exp 1 values have been evaluated with the experimental values of I and A, 1,40,41 while g 2 has been calculated using the orbital energies obtained at B3LYP/6-31111G(2d,2p) level.
It is worth to mention that the reported values of g J ðr Þ½q HOMO ðr Þ, g T0 ðr Þ½q HOMO ðr Þ, and g Ex ðr Þ½q HOMO ðr Þ are slightly different from those reported in Table 1 of ref. 24 whenever the HOMO orbital is degenerate. The difference resides in the way to evaluate the Fukui function. In our previous article, the Fukui function was taken simply as the density of the last occupied orbital printed in the checkpoint file of Gaussian, while in the present work it has been evaluated as the average of the degenerated HOMO orbitals, that is where d is the degeneracy of the HOMO level.
As one can see in Table 1, the two new terms (Wigner correlation and 1/9th Weizsäcker functionals) not considered in our previous works 24,42 are very small and represent less than 5% of the total global hardness. Then, the conclusions derived from Table 1 are similar to those of our previous studies, i.e., the Coulombic and Thomas-Fermi contributions are the dominant terms in eqs. (12) and (14). The Coulombic term represents between 35 and 88% of the total global hardness (from 65 to 88% if excluding Li 1 , Be 21 , and B 31 ), whereas, the Thomas-Fermi contribution gives from 15 to 68% (excluding Li 1 , Be 21 , and B 31 : from 37 to 68%). The exchange term (negative contribution) approximated by the Dirac functional is also relevant and it lies in the range from 1.5 to 7.2% of the total global hardness and is particularly relevant in systems like H 2 (21.7 eV), H 2 O (21.23 eV), or CS (20.96 eV). It is important to remark that for the Li 1 , Be 21 , B 31 , and H 2 systems the Weizsäcker correction term vanishes, because in these systems q HOMO ðr Þ is equal to qðr Þ. 14

Profiles of Local Hardness
Local hardness, gðr Þ, has been a source of great controversy in the literature 11,12,43,44 (a very recent review of this issue can be found in ref. 14), although it has become a reactivity index very useful in predicting the regioselectivity of chemical reactions, [45][46][47][48][49][50][51][52] particularly when gðr Þ is approximated by the electronic part of the molecular electrostatic potential 2 * This equation can be derived from eq. (9) if one considers only the Coulombic contribution to the hardness kernel and qðr Þ=N as approximation of the Fukui function. As far as we know, Figures 1-3 are the first profiles of gðr Þ with kinetic and exchangecorrelation contributions for polyatomic molecules, based on eq. (10). These figures recover some of the features displayed in the previous section. For instance, they illustrate that the Coulombic term is the dominant contribution in gðr Þ, with exception of the Li 1 case, and they also pinpoint the crucial role of the Thomas-Fermi term in the shape of the total local hardness. For instance, note the differences between the Coulombic and Coulombic-Thomas-Fermi gðr Þ profiles of Na 1 and K 1 in Figures 1b and  1c, respectively, where the shell structure of the atoms is only It is well-known that the largest value of f ðr Þ indicates the preferred site for interaction of one system with another, because in ''two different sites with generally similar dispositions for reacting with a given reagent, the reagent prefers the one which on the reagent's approach is associated with the maximum response of the system's chemical potential.'' 1,53 A similar expression to eq. (21) exists for gðr Þ 7 *  Then, in analogy to f ðr Þ, see eq. (21), one can conclude that gðr Þ indicates the most reactive site of a molecule induced by a change in its density, see eq. (22). This fact explains the maximum of gðr Þ displayed in Figures 2a, 2b, 3a, and 3b, which shows that the most reactive positions of the water and ethylene molecules are the regions of the lone pairs of the oxygen atom and the p-bonding orbital between the two carbons, respectively. In the gðr Þ profile parallel to and at a distance of 1 Bohr from the internuclear axis of HCN molecule, one can see that gðr Þ profile is larger at the carbon (located at 20.9392 a.u.) and nitrogen (located at 1.2268 a.u.) atoms than that at the hydrogen (located at 22.9523 a.u.) atom, indicating that two former atoms are the reactive places of the system. In addition, the nitrogen atom is harder than the carbon and gðr Þ presents a minimum approximately at the center of the CBN bond (at the origin of coordinates); although this minimum only appears when the Thomas-Fermi contribution in gðr Þ is included. The fact that the maximum of gðr Þ is associated with the most reactive positions of the molecule might seem to be in contradiction with the maximum hardness principle, 54,55 which affirms that the systems tend to a state of maximum hardness at constant temperature, external potential, and chemical potential. This principle implies that the hardest species are also the less reactive. From our point of view, gðr Þ indicates the region of the molecule that contributes more to the global hardness, which is at the same time the most reactive place upon a change in its density. Then, it is important to distinguish between a large value of g indicating a large stability of the total system (and a low global reactivity) and the maximum of its local counterpart, gðr Þ, representing in this stable molecule the most reactive region. .9392, and 1.226 Bohr, respectively. Note that in the hardness profile of water, the oxygen atom is at the origin and the axis is perpendicular to the molecular plane. All values are given in a.u.
As noted in the Theoretical Background and Computational Details section, Harbola et al. 12 argued that if the density of the ground state and the exact functional are used in eq. (9), the equivalent of gðr Þ of eq. (10) becomes constant in all the positions of the space and equal to the global hardness. In the results displayed in the Figures 1-3, one can see that the considered gðr Þ model of eq. (10) is quite far from being a constant in all positions of space. It should be noticed that the functional used to calculate eq. (10) is also far from being the exact functional and that there is an inconsistency in our calculations due to the differences between the functional (Thomas-Fermi-Dirac-1/9th Weizsäcker-Wigner) used to evaluate eq. (10) and the functional used to calculate qðr Þ and q HOMO ðr Þ and their derivatives (B3LYP level). Nevertheless, looking at Figures 1-3, it is difficult to envisage that the kinetic and exchange-correlation contributions will cancel the dominant and exact Coulombic term yielding a constant local hardness in all positions of the space. This would happen only if the normalized function that minimizes the variational principle of Chattaraj et al. 9 leading to both constant local gðr Þ and the ''best'' global hardness is very different from q HOMO ðr Þ. We think that such a possibility is unlikely. The present approach provides at different point of view of the local hardness with an approximate model, which is far from being exact, but allows us to obtain results with a chemical meaning, which can be successfully interpreted.

Condensed Atomic Hardness
The most popular CDFT reactivity indices, which contain relative information about different regions in a given molecule are the Fukui function, f ðr Þ, 53 and the local softness, sðr Þ. 56 The majority of applications of f ðr Þ and sðr Þ are reported 57-64 on a condensed form, 65 that is, using the finite difference approximation and integrating over atomic regions by means of some population analysis techniques. In the literature, one can find a large amount of techniques in population analysis (Mulliken,66 CHELPG, 67 MK, 68 natural population analysis, 69 atom in molecules, 70 Hirshfeld 71 ) and the choice of the correct partitioning scheme to describe the condensed Fukui function and local softness has been analyzed in detail and numerous works have been reported. [72][73][74][75][76][77][78][79][80][81][82][83][84][85] The two main problems of these reactivity indices are the basis set dependence often associated with atomic population analysis and their physical meaning when negative values are obtained.
Here, we introduce a new atomic condensed index that is easily derived from the methodology described above. In eq. (11), the global hardness is obtained by integration over the whole space of the product of two local functions, namely the local hardness and the Fukui function. One can perform the integration over a particular region of the space and obtain its contribution to the global value. Then, by taking appropriate atomic domains defined in one or another way one can define the condensed atomic hardness contributions, g i , to the global hardness value as in such a way that where the g i quantities are examples of the so-called fragment of a molecular response type. 85 In this work, we have used the fuzzy Voronoi polyhedra 30 to define the atomic domains, which are readily obtained from our numerical integration of eq. (12). Such atomic definition has already been successfully applied for the calculation of overlap populations, bond orders, valences, 86 or in several molecular energy decomposition schemes. 87,88 Any other disjoint (Atoms in Molecules 70 ) or fuzzy (Hirshfeld 71 ) decomposition of the 3D space could also be used. In any case, the condensed atomic hardness values are expected to be fairly basis set independent, as usual for 3D space partitioning methods. Table 2 contains the condensed atomic hardness of the oxygen and nitrogen atoms for seven molecules [H 2 O, CH 3 OH, CH 3 OCH 3 , NH 3 , NH 2 CH 3 , NH(CH 3 ) 2 , and N(CH 3 ) 3 ] evaluated at B3LYP level with three different basis sets: 6-31G, 6-311G(d), and 6-31111G(2d,2p). Considering the reference values as the results obtained with the largest basis set, one can realize that the differences in the condensed atomic hardness are smaller than 20% in 6-31G basis set and smaller than 3% in 6-311G(d), showing little basis set effects of the condensed atomic hardness. Furthermore, g X values are always positive and they show the same tendencies for the three basis sets, i.e., g X decreases with the number of methyl groups on oxygen or nitrogen in analogy to the global hardness. Note that the oxygen values are always larger than the nitrogen for the same number of methyl groups carried by them in line with the global hardness of water and ammonia.
Then, one can conclude that the condensed atomic hardness can be efficiently used as a CDFT local reactivity index and additional applications of the condensed atomic hardness are presently under study.

Conclusions
In this article we have tested the combination of approximate hardness kernel and Fukui function for the evaluation of global and local hardnesses. The hardness kernel and Fukui function have been approximated using the second-order derivative of the Coulombic-Thomas-Fermi-1/9th Weizsäcker-Dirac-Wigner functional with respect to the density and the density of the HOMO, respectively. We have shown that the dominant terms are the Coulombic and Thomas-Fermi contributions, although the remaining contributions cannot be neglected and must be included to obtain accurate values of global and local hardnesses.
In addition, 2D and 3D local hardness profiles have been reported and an analogy between the Fukui function and the local hardness has been established, showing that the local hardness indicates the regions of the molecule that contribute more to the global hardness and which are the most reactive sites of a molecule due to a change in its density. Finally, the concept of condensed atomic hardness has been introduced, paying special attention to its basis independent character and its potential as CDFT local reactivity index. Condensed atomic hardness obtained from the numerical integration of eq. (12) in the fuzzy Voronoi polyhedra of oxygen and nitrogen regions.