Relations among Several Nuclear and Electronic Density Functional Reactivity Indexes

An expansion of the energy functional in terms of the total number of electrons and the normal coordinates within the canonical ensemble is presented. A comparison of this expansion with the expansion of the energy in terms of the total number of electrons and the external potential leads to new relations among common density functional reactivity descriptors. The formulas obtained provide explicit links between important quantities related to the chemical reactivity of a system. In particular, the relation between the nuclear and the electronic Fukui functions is recovered. The connection between the derivatives of the electronic energy and the nuclear repulsion energy with respect to the external potential offers a proof for the ''Quantum Chemical le Chatelier Principle.'' Finally, the nuclear linear response function is defined and the relation of this function with the electronic linear response function is given.


I. INTRODUCTION
According to Density Functional Theory ͑DFT͒, 1 the electronic energy of a many-electron system is given by E͓͔ϭF͓͔ϩ ͵ ͑r ជ ͒͑ r ជ ͒dr ជ,

͑1͒
where F͓͔ is the universal Hohenberg-Kohn functional defined as the sum of the electronic kinetic energy functional, T͓͔, and the electron-electron interaction energy functional, V ee ͓͔, and (r ជ) is the external potential that for an isolated molecule is just the potential due to atomic nuclei, but in general it may also include contributions from external fields or neighboring molecules.Minimization of the electronic energy with respect to the electron density with the constraint of constant total number of electrons, N, leads to an Euler equation of the form ϭ͑r ជ ͒ϩ ␦F͓͔ ␦͑r ជ ͒ .

͑2͒
The solution to this equation yields the ground state electron density, from which one can determine the ground state energy.
Since the number of electrons and the external potential completely determine the Hamiltonian of the system, one can write the electronic energy of the system as a functional of N and (r ជ) ͑i.e., EϭE͓N,͔).On the other hand, Eq. ͑1͒ establishes that E is a functional of (r ជ) and (r ជ) ͑i.e., E ϭE͓,͔).Differential expansion of the electronic energy, E, for these two energy functionals allows us to determine the physical significance of the Lagrange multiplier, , which is found to be the electronic chemical potential of the system, that is, and to find that the density is the functional first derivative of the energy with respect to the external potential at a constant total number of electrons, i.e., ͑r Second-order derivatives of the energy with respect to N and (r ជ) also have clear physical meanings. 2Thus, the second derivative with respect to the number of electrons at a fixed external potential (r ជ) is the hardness: 3 ͑r ជ,r ជЈ͒ϭ which measures the variation of the electron density at point r when the external potential of the system is locally perturbed at rЈ.Finally, the Fukui function, first defined by Parr and Yang, 5 is the second derivative with respect to external potential (r ជ) and the number of electrons, The Fukui function describes the local changes in the electron density of the system, due to the perturbation in the global number of electrons, so it reflects the character of a molecule to accept ͑donate͒ electrons from ͑to͒ an another system.For a molecular or atomic system, the density is a discrete function of the number of electrons.Because of that, Parr and Yang 5 associated different physical meanings to the left, right, and central derivative of the density with respect to N, corresponding to a reactivity index for a nucleophilic f ϩ (r ជ), electrophilic f Ϫ (r ជ), and radical f 0 (r ជ) attacks, respectively.By applying a finite difference approximation to Eq. ͑7͒, these three approximate Fukui functions can be written as where Nϩ1 (r ជ), N (r ជ), and NϪ1 (r ជ) are the electronic densities of the system with Nϩ1, N, and NϪ1 electrons, respectively.
Upon the course of a chemical reaction, molecules adjust their number of electrons and their external potential.This implies changes in both the electronic structure, characterized by the electron density (r ជ), and the geometric structure defined by the position vectors ͕R ␣ ͖ of the nuclei.To describe the electronic changes, one can use the Fukui function, the hardness, and the linear response function.These are electronic indexes that measure the electron density and electronic chemical potential responses to a change in N or (r ជ).A Kohn-Sham formulation of the chemical electronic responses has been provided by Senet. 6To describe the geometrical changes, a set of nuclear reactivity indexes that characterize the response of the nuclei due to changes in N or (r ជ) of the system have been recently defined. 7,8In particular, the nuclear Fukui function ͑NFF͒ has been defined by Cohen and coworkers 7 as the electronic force F acting on the nuclei due to the perturbation in N at a constant (r ជ): In analogy with the electronic Fukui function, it is possible to define three NFFs corresponding to nucleophilic, electrophilic, and radical attacks as 9 The last equality in Eqs.͑12͒ and ͑13͒ is only valid when the N electron system of reference is at its equilibrium geometry ͓F ជ ␣ (N)ϭ0͔.Baekelandt 8 via a Maxwell relation has shown that this NFF can be interpreted as the conformational contribution to the change in the electronic chemical potential: From Eq. ͑15͒, three definitions of the NFF corresponding to the left, right, and central derivatives are also possible: 10 where the Koopmans' theorem 11 has been used in Eqs.͑16͒ and ͑17͒, and 2 and 1 in Eq. ͑18͒ refer to the following operational approximations for the electronic chemical potential: It is remarkable the similarity between Eqs. ͑7͒ and ͑15͒.
Geerlings and coworkers have computed NFFs for diatomic molecules using Eqs.͑12͒ and ͑13͒, as well as from numerical calculations of the negative derivative of the chemical potential with respect to the atomic coordinates ͓Eqs.͑16͒ and ͑17͔͒ 9 and with an approximate analytical method 10 developed in analogy with Komorowski and Balawender's 12 coupled Hartree-Fock approach to the electronic Fukui function.They have found that analytical and numerical results show a high correlation.6][17] These studies 8 -12,14 -17 have reinforced the idea that a complete representation of the total chemical response to a given perturbation must involve the analysis of both electronic and nuclear reactivity descriptors.Because of the coupling between the electronic and nuclear responses to external perturbations, one can expect that relationships between nuclear and reactivity indexes exist.In the present paper we apply the general functional expansion scheme to expand the potential energy, U, and the total electronic energy, E, of the systems in terms of the total number of electrons, N, and normal coordinates displacements, ͕Q k ͖, to derive relationships between electronic and nuclear indexes that connect the electronic and nuclear responses in front of a given perturbation.The energy expansion is carried out within the canonical ensemble, for which the natural variables are the global variable N and the local variable (r ជ).

II. ENERGY REPRESENTATIONS
In the Born-Oppenheimer approximation, the Potential Energy Surface ͑PES͒ of a system can be divided into an electronic (E) and nuclear repulsion (V nn ) energy terms as U͓,͔ϭE͓,͔ϩV nn ͓͔.͑21͒ As stated for the electronic energy, one can also write the PES of the system as a functional of N and (r ជ), i.e., U ϭU͓N,͔:

U͓N,͔ϭE͓N,͔ϩV nn ͓͔. ͑22͒
A Taylor series functional expansion of the PES around a reference external potential o (r ជ) and total number of electrons N, retaining terms up to second order, acquires the form ͑note that the nuclear repulsion term does not depend on N)

͑23͒
Equation ͑23͒ is the same used by Ayers and Parr to explore the effect that changing the external potential has on chemical reactivity. 18Now, substituting Eqs.͑3͒ to ͑7͒ into Eq.͑23͒ one finds that

͑24͒
Let us now consider a molecule in an isotropic and homogeneous environment.In this case, the external potential depends only upon the relative nuclear positions.These relative positions can be specified using a set of 3N-6 ͑or 3N-5 for linear molecules͒ coordinates, the most natural choice being the use of the normal coordinates, ͕Q k ͖.Given that the external potential for an isolated molecule can be defined from a set of ͕Q k ͖ coordinates, one can write the energy of this system as a functional of N and ͕Q k ͖, i.e., U ϭU͓N,͕Q k ͖͔.In this functional form, a Taylor series expan- sion of the total energy up to second order yields

͑25͒
⌬Q k being the displacement of the normal mode k from the equilibrium position.It is worth noting that since the external potential can be written as a function of the normal coordinates, it is possible to transform the derivatives with respect to the external potential to derivatives with respect to normal coordinates by using the chain rule.In this way, it is also possible to derive Eq. ͑25͒ from Eq. ͑23͒.
If we take as the reference state a molecule in its equilibrium geometry, then the term in Eq. ͑25͒ and the term in Eq. ͑24͒ are zero.We note in passing that Eq. ͑27͒ leads to the following interesting relationship that holds for an isolated molecule in its equilibrium geometry: This equation is the result of the fact that, in the equilibrium geometry, the electronic energy change due to a perturbation in (r ជ) equals the negative of nuclear repulsion energy change, providing a relation between electronic and nuclear reactivity descriptors.Further, it is a nice proof to the firstorder and within the framework of conceptual DFT 19 of the ''Quantum Chemical le Chatelier Principle'' formulated by Mezey some years ago. 20This principle states that the electronic energy change induced by an external perturbation is compensated by a modification of the nuclear repulsion energy by a similar amount ͑in absolute value͒, in such a way that the value of the total energy (electronicϩnuclear) is approximately conserved.Now, from Eq. ͑25͒, using Eqs.͑3͒ and ͑5͒ and taking into account that and Equation ͑31͒ holds for a molecule in its equilibrium geometry.In Eqs.͑29͒ and ͑31͒, k are the NFFs in terms of the normal coordinates and, in Eqs.͑30͒ and ͑31͒, k are the force constants of each normal mode.
A comparison of Eqs.͑24͒ and ͑31͒ leads to Eqs. ͑32͒ and ͑34͒.Equation ͑32͒ is obtained, considering the independent terms with respect to ⌬N: and defining the nuclear linear response function as while Eq.͑34͒ appears matching the terms that depend on ⌬N: Equation ͑32͒ affords a relation for the nuclear and electronic linear response functions.For small ͕Q k ͖ displace- ments, the term (1/2)͚ kϭ1 3NϪ6 k ⌬Q k 2 in Eq. ͑32͒ gives the total change in the potential energy due to nuclear rearrangement.Equation ͑32͒ shows that this change can be decomposed in two parts: the first term, (1/2)͐͐X(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ, represents the change of the nuclear potential energy due to the distortion while the remaining term is the change in electronic energy.
Equations ͑32͒ and ͑34͒ are the result of the Taylor series expansion of the PES functionals UϭU͓N,͔ and U ϭU͓N,͕Q k ͖͔.Similarly, one can expand the functionals of the electronic energy EϭE͓N,͔ and EϭE͓N,͕Q k ͖͔.Fol- lowing the same procedure one recovers Eq. ͑34͒ and also a modified version of Eq. ͑32͒ that reads as 1 2 ͵ ͵ ͑r ជ,r ជЈ͒⌬͑r ជ ͒⌬͑ r ជЈ͒dr ជdr ជЈ Considering now the relation between the electronic energy and the PES given by Eq. ͑22͒, one obtains This equation provides a way to compute the (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ integral from the evaluation of the other two terms in Eq. ͑36͒.Remarkably, using Eqs.͑36͒ and ͑32͒, it is possible to compute numerically all terms present in Eq. ͑24͒.By means of Eq. ͑24͒, one could carry out an energy decomposition analysis of the interaction between two given systems.More research on the value of such an energy decomposition analysis is underway in our laboratory.
Coming back to Eq. ͑34͒, this equation represents a basic expression that connects the electronic and nuclear Fukui functions.It is a quantitative expression for the coupling of nuclear and electronic responses to perturbations.It was first found by Baekeland 8 in differential form and in a different way, starting from the definition of the NFF and using the chain rule to express the derivative of the electronic chemical potential with respect to the external potential.According to this author the term Ϫ ͚ kϭ1 3NϪ6 k ⌬Q k corresponds to the nuclear-dependent part of the change in the electronic chemical potential.Interestingly, we have recovered this important connection between electronic and nuclear Fukui functions through a Taylor functional expansion of the energy as a functional of N and ͕Q k ͖.

III. CALCULATION OF THE TERMS INVOLVED IN THE RELATIONS BETWEEN NUCLEAR AND ELECTRONIC REACTIVITY INDEXES
Throughout this paper we have assumed that all Taylor series converge.Moreover, we expect that if the external potential of the final system resembles that of the reference state sufficiently closely, low-order truncations of these functional Taylor series will be accurate.Our aim in this section is first to check this hypothesis by carrying out calculations of the two terms of Eq. ͑34͒.Second, the calculation of the (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ and (1/2)͐͐X(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ integrals will be performed for a series of diatomic molecules.
All calculations have been carried out with the GAUSS- IAN 98 package 21 at the Hartree-Fock ͑HF͒ 22 level using the aug-cc-pVQZ basis set, 23 except for the lithium atom, for which the cc-pVQZ has been used.The energy and gradient of the cationic and anionic doublet species has been computed within the unrestricted methodology at the geometry of the neutral systems, while the neutral singlet molecules have been calculated using the restricted formalism.
The numerical differentiation of HOMO , LUMO , 1 , I, A, and 2 needed to evaluate Eqs.͑16͒-͑18͒ has been carried out performing displacements of the equilibrium geometry of Ϯ(1,2,4,8,16,32,64)ϫ10 Ϫn (nϭ2 to 5͒ bohrs.Then, the smallest magnitude displacement that produced a stable derivative has been selected using a Romberg method triangle. 24efore starting the evaluation of the Ϫ ͚ kϭ1 3NϪ6 k ⌬Q k term in Eq. ͑34͒, we have checked for 13   I. As already pointed out by Baekelandt 8 and Geerlings and coworkers, 9,10 NFF data can be related to the Berlin's function 8,13,14 to analyze bonding in molecules.Thus, in diatomic molecules a positive nucleophilic NFF is associated to an increasing bond length upon the addition of an electron, which increases the total electron density in antibonding regions.All species studied show positive nucleophilic NFF as indicated by the values obtained from the (␦A/␦R ជ ␣ ) N and (‫ץ‬F ជ ␣ /‫ץ‬N) ϩ derivatives.Likewise, a negative electrophilic NFF is also related to an increase of the bond length during ionization, which means a reduction of the total electron density in bonding regions upon ionization.We found negative electrophilic NFF for all species studied, except for BH, F 2 , CO, NF, and BF, for which positive electrophilic NFFs are obtained, indicating that the ionization should lead to a decrease of the bond length in these molecules.
Let us now numerically test the validity of Eq. ͑34͒.We have checked Eq. ͑34͒ for different displacements using the NFF obtained from Eqs. ͑12͒ to ͑14͒, although we only show in Table II the results for ⌬Q k ϭ1ϫ10 Ϫ5 ϫͱ a.u., where is the reduced mass, because with this displacement the error percentage becomes stable.We note that, from Eq. ͑8͒, the first term appearing in Eq. ͑34͒ is

͑38͒
where N 0 (r ជ) and Nϩ1 0 (r ជ) are the electronic density at the equilibrium geometry for a system with N and Nϩ1 electrons, respectively, while 0 (r ជ) is the external potential at the equilibrium geometry and (r ជ) is the external potential obtained increasing or decreasing the equilibrium bond length.The ͐ Nϩ1 0 (r ជ) 0 (r ជ)dr ជ and ͐ N 0 (r ជ) 0 (r ជ)dr ជ integrals are explicitly the electron-nucleus potential energy of the system with Nϩ1 and N electrons at the equilibrium geometry, respectively.The ͐ Nϩ1 0 (r ជ)(r ជ)dr ជ and ͐ N 0 (r ជ)(r ជ)dr ជ integrals in Eq. ͑38͒ are the electron-nucleus potential energy introducing the modified external potential (r ជ), but using the Nϩ1 and N electronic density at the equilibrium geometry.All integrals present in Eq. ͑38͒ have been obtained from electron-nucleus potential energies computed with the GAUSSIAN 98 program. 21,25Equivalent expressions to Eq. ͑38͒ can be used to evaluate the integrals containing f Ϫ (r ជ) and f 0 (r ជ).
As can be seen in Table II, the difference between the first and the second term of Eq. ͑34͒ is usually quite small for all systems studied and for the different nucleophilic, electrophilic, and radical approximations to the electronic and nuclear Fukui functions.The similarity between the values of the NFF computed using the two most reliable approximations make us confident about the accuracy of the calculated NFFs, while operational equations that provide the electronic Fukui function ͓Eqs.͑8͒ to ͑10͔͒ are unavoidably approximate.Thus, in our opinion, the numerical errors that produce the differences observed may be larger in the computation of the integrals of the type ͐ f (r ជ)⌬(r ជ)dr ជ than in the evaluation of the Ϫ ͚ kϭ1 3NϪ6 k ⌬Q k term.Figure 1 depicts the good correlation observed between the ͐ f Ϫ (r ជ)⌬(r ជ)dr ជ and the Ϫ ͚ kϭ1 3NϪ6 k Ϫ ⌬Q k values giving an equation with a slope very close to 1, which is the expected value according to Eq. ͑34͒.Similar results are obtained for the analogous relations involving the nucleophilic and radical Fukui functions.We have checked Eq. ͑34͒ using different basis sets and we have found that, in general, the correlation between the ͐ f (r ជ)⌬(r ជ)dr ជ and the Ϫ ͚ kϭ1 3NϪ6 k ⌬Q k values improves when the quality of the basis set used increases.Finally, we have also tested Eq. ͑34͒ using the B3LYP method and we have found similar results to those given in Table II.
As a final point, Table III gathers the computed values for the (1/2) ͚ kϭ1 3NϪ6 k ⌬Q k 2 and (1/2) ͚ kϭ1 3NϪ6 ͚ lϭ1 3NϪ6 ‫ץ(‬ 2 V nn / ‫ץ‬Q k ‫ץ‬Q l ) N ⌬Q k ⌬Q l terms appearing in Eq. ͑36͒.The integral (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ is evaluated as the difference between the two previous calculated terms, according to Eq. ͑36͒.It is also worthwhile noting that in line with Eq. ͑37͒, the second term corresponds to the integral (1/2)͐͐X(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ.From Eq. ͑32͒, one can see that, starting from a system in its equilibrium geometry, the change ͑up to second order͒ in the energy of the molecular system due to a molecular distortion can be split into an electronic and a nuclear term.The electronic term, (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ, is the change in the electronic energy due a change in the external potential, while the nuclear part, (1/2)͐͐X(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ, gives the change in nuclear repulsion energy due to the same perturbation.Interestingly, for the diatomic molecules studied, the change in electronic energy due to molecular distortion is always negative and smaller in absolute value than the nuclear energy variation, which is always positive ͑see Table III͒.It would be interesting to investigate whether the sign of these integrals are preserved in polyatomic molecules.Remarkably, those molecules that suffer a large nuclear energetic change in the distortion ͑for instance, the N 2 , F 2 , CO and NO ϩ species͒, also experience a large electronic energetic change.Let us finish by mentioning that the values of (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ in Table III can be a good starting point for discussing the validity of possible analytical expressions for the electronic linear response function.

IV. CONCLUSIONS
In this paper, we have derived a set of connections among several nuclear and electronic indexes of reactivity in the framework of the conceptual Density Functional Theory by using an expansion of the energy functional in terms of the total number of electrons and the normal coordinates within a canonical ensemble.The relations obtained provide explicit links between important quantities related to the chemical reactivity of a system.
First, we have demonstrated that the derivative of the electronic energy with respect to the external potential of a system in its equilibrium geometry is equal to the negative of the nuclear repulsion derivative with respect to the external potential, giving a proof up to the first order for the so-called ''Quantum Chemical le Chatelier Principle.''Second, following this particular approach, we have given an alternative proof the relation between the nuclear and the electronic Fukui functions, and, for the first time, we have numerically checked this relation.Finally, the nuclear linear response function has been introduced and a relation of this function with the electronic linear response function and the force constant has been provided.
diatomic molecules the consistency of the ␣ ͑12͒ to ͑14͒ and ͑16͒ to ͑18͒.The results obtained are listed in TableI.It is remarkable to see that calculations of the nucleophilic NFF, ␣ ϩ , ␣ Ϫ , and ␣ 0 nuclear Fukui functions calculated using the different approximations repre-sented by Eqs.ϩ , by numerical differentiation of (␦A/␦R ជ ␣ ) N or (‫ץ‬F ជ ␣ /‫ץ‬N) ϩ yield practically the same values providing reliability on the computed values; while the calculation of ␣ ϩ using the orbital approximation Ϫ(␦ LUMO /␦R ជ ␣ ) N results in somewhat different numbers.The same is found for the electrophilic and radical NFFs.It is worth noting that the results obtained by Balawender and Geerlings 10 for ␣ ϩ and ␣ Ϫ of H 2 , LiH, BH, HF, HCl, N 2 , F 2 , CO, BF, and LiF at the HF/6-31ϩϩG** are, in general, in good agreement with those reported in Table

TABLE II .
The calculation of the two terms given in Eq. ͑34͒ for the three approximations of the electronic and nuclear Fukui functions for 13 diatomic molecules.All values are given in atomic units.The values of nuclear Fukui functions are obtained using Eqs.͑12͒, ͑13͒, and ͑14͒.The ⌬Q k employed in all these cases is 1ϫ10 Ϫ5 ϫͱ a.u., where is the reduced mass. a

TABLE III .
The calculation of the (1/2)͐͐(r ជ,r ជЈ)⌬(r ជ)⌬(r ជЈ)dr ជdr ជЈ integral from the evaluation of the remaining two terms in Eq. ͑36͒ for the 13 diatomic molecules studied.All values are given in atomic units.a The ⌬Q k employed in all these cases is 1ϫ10 Ϫ5 ϫͱ a.u., where is the reduced mass. a