From the Electron Density Gradient to the Quantitative Reactivity Indicators: Local Softness and the Fukui Function.

Important reactivity measures such as the local softness, the Fukui function, and the global hardness have been calculated directly from first principles with the use of the electron density function, beyond the finite difference approximation. Our recently derived density gradient theorem and the principle of nearsightedness of the electronic matter have been instrumental in obtaining the original, albeit approximate, result on the local softness of an atom. By integration of the local softness s(r), we obtain the global softness S and the Fukui function f(r) = s(r)/S. Local and global softness values have also been calculated analytically for the basic hydrogenic orbitals; the general relation to the atomic number S = σZ -2 has been demonstrated, with constants σ characteristic for each orbital type. Global hardness η = 1/S calculated for atoms and ions has been favorably tested against its conventional measure given by the finite difference approximation: (I - A). Calculated test results for atoms and ions in rows 1-4 of the periodic table have been presented.

Introduction

The essence of a chemical reaction is the variation of a molecular structure accompanied with the flow of electrons from initial positions to those in the final state. Two approaches have been in use for describing the reactivity. The physical one deals with the transformation of a supermolecule—a system created from two or more molecules (reactants) that approach each other to form the final state, the molecule being the product. Quantum mechanical calculations allow for revealing the structure of the product, the transition state, and the path through the potential energy surface. Recently, we have developed a new theoretical procedure that contributes to this physical approach. We have introduced the reaction fragility method for monitoring the alternation of bonds and to reveal the mechanism of the chemical reaction.[1−4]

The chemical approach deals with a molecule being attacked by another object; reactivity indexes are in use to predict the direction of the transformation. The chemical way is based on the concept of regional (site) reactivity, and density functional theory (DFT) has strengthened this approach, since electron density, according to Hohenberg and Kohn theorems, contains all information on the properties of the electron gas including various types of responses to the external perturbation. Here, in this paper, we focus mainly on this chemical perspective.

The lasting step toward the quantitative measure of site reactivity in terms of the electronic configuration of interacting objects has been done by Fukui in his MO theory of reactivity, by indicating the role of frontier orbitals.[5] DFT provided theoretical support for this concept: the Fukui function index has been proposed by Parr and Yang[6] as the response function of the electron density function under the perturbation of the number of electrons, linking early DFT to theory by Fukui. Another important idea of reactivity encompassed by conceptual DFT was the chemical hardness first proposed by Pearson who introduced the hard and soft acid and base (HSAB) principle.[7] Parr and Pearson[8] have identified molecular hardness in cDFT as the second derivative of energy with respect to the number of electrons, and softness as the inverse of hardness. Nalewajski has analyzed the creation of a chemical bond with the HSAB principle by substituting electrostatic energy to the formula of Parr and Pearson.[9] This was the origin of the charge sensitivity analysis later developed by Nalewajski and coworkers.[10,11] These ideas have attracted attention for the following decades, both from computational and chemical perspectives, as testified by the reviews of the subject.[12−15] However, despite many efforts for attribution of numerical values to the Fukui function and hardness, the results have been discouragingly far from expectations of practical chemistry; the necessity for differentiating versus the integer number of electrons presented an insurmountable obstacle inherent in chemical applications of the cDFT formalism.[16,17]

Further development proved that the average of the directional derivatives f– and f represents exact values of the Fukui function at zero temperature.[18] Describing the site reactivity in these terms requires building an ensemble (grand canonical) containing the states of N, N – 1, and N + 1 electrons. Methodological consequences are significant: one cannot describe the site reactivity relying only on just ground-state density. In this present work, an alternative approximate method has been proposed. A working formula for the Fukui function which uses only the ground-state density (gradient) has been proposed and confronted with the directional derivatives f– and f at the same computational level of accuracy.

The only approximation for the Fukui function appreciated by chemists is that based on the densities of the frontier orbitals ρHOMO and ρLUMO as proposed by the original authors,[8] and the sole practical way to characterize the hardness scale for atoms remained the oversimplified model of the parabolic E(N) function that leads to absolute hardness proportional to (I–A).[8] Chemical applications of the derivatives of the energy/density have also been hampered by the increasing complexity in their calculation.[19] Yang et al. have proposed a scheme based on the response function for the construction of the fractional number of electrons to perform differentiation with respect to N.[20,21] Parr and Yang[6] realized that a chemical reaction involves either an increase or decrease in the number of electrons in a reactant as a result of electrophilic or nucleophilic attack, respectively; this fact has motivated the use of right- and left-hand derivatives for the Fukui function.

The potential hidden in the electron density concept has been reviewed by Bader in his late paper,[22] presenting an intriguing legacy: quantum chemistry has not yet accomplished the basic expectation of practical chemistry to quantify properties of its fundamental objects, atoms, and bonds that hold them together in molecules, except the widely accepted view that their nature must be determined by the electron density around the nuclei. Another paper by the same author brought a report from his much interesting dispute with the personality on the scene of chemical theory.[23] An experienced author had declared his belief: “chemical observations and hence the conceptual framework of chemistry are fully described, predicted, and accounted for by quantum mechanics–the physics governing the behavior of electrons and nuclei” as opposed to the cited declaration by Hoffmann who believed that “the most interesting ideas of chemistry are not reducible to physics”. Both arguments can be translated to the DFT language: by the Hohenberg Kohn theorem,[24] the electron density function is sufficient to quantify local reactivity properties in molecules, however, in some indirect and complicated manner.

The review on bond and bonding published few years later by Politzer and Murray added a significant point to this debate.[25] They noted its clear relevance to the concept of the bonded atom: “if we cannot unambiguously define atoms in molecules, we cannot rigorously define bonds between them”.[26] The authors share the view of Ruedenberg et al.: “the theoretical identification of atoms in molecules and a corresponding analysis of their interactions have however proved to be challenging. This is because the resolution of molecules in terms of atoms is not fundamental to rigorous physical theory”.[27]

Remarkably, these skeptical views were voiced in the period when conceptual DFT was reaching its maturity. The focus of theory has been recapitulated in a synthetic way by Liu, largely in tune with the optimism of Bader’s remarks:[28] “according to the basic theorems of DFT, the electron density of a system should contain adequate information to determine everything in the ground state, including all properties related to stability, bonding, and chemical reactivity”. Perspectives and drawbacks of the theory have been collected in the latest joint review by a group of authors active in the field.[29] The long-desired applications of the cDFT concepts in chemistry are also the subject of their analysis: the authors admit that the current role of theory remains predominantly interpretative, not predictive. Specifically, their discussion of the chemical hardness has unveiled the lasting weakness of the theory: the hard and soft acid and base principle has been recognized as a practical guide in chemistry for more than half a century,[30,31] and yet the theory has been unable to provide a coherent measure of hardness based on the electron density for atoms or molecules even for this apparently well-defined quantity.

An original approach has been published by Rahm and Hoffmann who introduced some special ways of partitioning molecular energy and defined a reactivity index as an average binding energy of electrons.[32] The index was called electronegativity, and the whole scheme has been introduced under the name of experimental quantum chemistry. The authors developed this scheme even further as to produce electronegativity of atoms under pressure.[33] The drawback of this approach is the rather arbitrary use of the term electronegativity and the lack of any relationship to other reactivity measures, for example, the Fukui function within cDFT that gives precise definitions of the reactivity descriptors derived by phenomenological procedures with a single obstacle: the inevitable differentiating with respect to N.

This has motivated the goal of this present work: presenting the unambiguous measure of the recognized potential reactivity index—the local softness, global softness (inverse hardness), and consequently the Fukui function. The necessary theoretical analysis has been based on the results of monitoring atoms in reacting molecular systems by the reaction fragility method,[4] the practical tool for describing both the reactivity and mechanism of the chemical reaction.[1−3] It is successfully applicable to various types of chemical processes;[34] however, for obtaining the mechanism of the chemical reaction with the reaction fragility method, one needs to actually perform the IRC calculation of a reaction path. Results of the method are parallel to those by Kraka who studied normal mode evolution along the reaction path.[35]

In this present work, we have extended the method toward predicting the reactivity instead of describing the evolution of atoms and bonds involved in the process. The basic elements for this prediction have been provided in the form of the precise maps of the local softness and Fukui function, the local reactivity indices. Our aim was to fulfill the lack of rules derived from cDFT to be used by chemists (both theoreticians and experimentalists) as a reliable and appropriate interpretational and predictive tool.[29] Our present work delivers direct and accurate, although not exact, link between the electron density function and local softness, without a need for differentiating over N. Taking the results and remarks from Chattaraj, Cedillo, and Parr[36] on obtaining the Fukui function from the softness kernel, we come up with a solution to the problem left open by these authors. We applied our gradient theorem (eq ) and with the local approximation suggested by Vela and Gazquez,[37] we come up with a novel working formula for the local softness. Once local softness is achieved, we gain access to most (if not all) reactivity indices of cDFT. Presenting the global hardness and the Fukui function obtained from the local softness, the first step is made to bring cDFT of chemical reactivity closer to chemical practice: the basic tools for indexing chemical reactivity are ready for experimental tests.

Fukui Function

The Fukui function index f(r) = [∂ρ(r)/∂N]v is linearly related to the local softness: s(r) ≡ −[δN/δv(r)]μ, and the inverse of global softness is the linear coefficient: f(r) = s(r)/S. The Fukui function has been considered to be the cornerstone for the reactivity descriptors in cDFT in a similar manner as the chemical potential μ = [∂E/∂N]v provided a ground for theoretical interpretation of the electronegativity concept.[38] Both quantities need differentiation versus the integer variable, that is, the number of electrons. This is usually done with finite difference approximation, which results in the analysis of the derivatives having serious limitation since the three states needed for this procedure are far away from each other (ΔN = ±1). Efforts to overcome this by calculating the Fukui function solely from the electron density function of a system were undertaken by many authors: by local density approximation,[39] gradient expansion,[40] analytic methods,[41,42] ab initio calculations,[43,44] and polarizability analysis.[37,45−48] These theoretical studies have not found their way into chemical practice due to weak theoretical and conceptual basis and lack of any link to the well-established orbital-based picture of chemical interactions between atoms.

Unlike that Bader’s QTAIM theory focused on atoms,[49] cDFT as founded by Parr has been dealing predominantly with the local analysis:[50] the electron density function ρ(r), the external potential ν(r), and the derivatives of energy E[N,ν(r)] and electron density (the tentative measures of reactivity) are determined at each point in space around atoms.[51] Contracting them to atomic quantities has been an unresolved problem, despite a number of occasional ingenious inventions by many authors.[19,52−55]

The alternative, atomic perspective in DFT, first proposed by Cohen et al., has been based on the analysis of the Hellmann–Feynman (H–F) force:[56,57]The idea was further developed by present authors.[58−65] This approach has been recently appreciated as the valuable source of information on bonded atoms changing their roles in reacting systems.[4] The atomic positions {R} have been used as basic parameters for the energy function E[N,{R}], and consequently, the electric field of atomic nuclei {(r)} has replaced the external potential ν(r) as the key local quantity. Bonded atoms are most naturally identified by the H–F forces acting on their nuclei. The divergences of the H–F force vector over the nuclear displacements ∇·F ≡ ∂F/∂R have been proved to contain the integrated information on the electron density around nuclei, equivalent to the local derivatives of energy and density, thus opening the access to monitoring individual atoms along a reaction path,[2,4] with no need of artificial separation of atoms in molecules.

The abovementioned approach has been fundamental for the method of this present work. The basic property of the H–F force divergence, vanishing the nuclear term ,[2,66] was essential for proving the link between the density gradient ∇ρ(r) and the corresponding density gradients upon the shift of atoms: ∇ρ(r)≡∂ρ(r)/∂R.[4] Together with the existing cDFT relation between ∇ρ(r) and the linear response function ω(r,r′), this allowed for developing the new density gradient theorem and opened the way to overcome the obstacle of differentiating over the integer function of the number of electrons N.

Method: from the Electron Density Gradient to Local and Global Softness

When a molecule is considered with no other external fields other than those generated by its nuclei (r), the divergences of H–F forces on atoms fulfill two significant conditions proved in our recent work[2]

The relation of ∇ρ(r) and the electronic response function ω(r,r′) has been well-founded:[2]Since the sum of all atomic forces vanishes, the sum of their divergences vanishes as well: ∑∇·F = 0. By combining this sum rule with eqs and 2, the novel condition for the stationary electron density is

The solution for ρ(r) in the electronic stationary state must be unique by the Hohenberg and Kohn theorem. Since eq holds separately for every atom (A) in a system, electron density in the external field from all nuclei must comply to the condition

Equation represents the important relation between the density gradient and the external electric field generated by the nuclei, with the linear response function playing a key role. Using the Berkowitz and Parr relation,[67] it is straightforward to show the identity[4]Here, (r) = ∑(r) stands for the total electric field and s(r,r′) is known as the softness kernel. This exact density gradient theorem holds for any system of atoms in the absence of external fields other than from atoms within. Equation allows for exploration of the specific properties of ω(r,r′) and s(r,r′) kernels that have recently been a target for theoretical analyses, and valuable difference between both has been disclosed. By Kohn’s principle of nearsightedness of electronic matter,[68] discussed by Bader[69] and recently reminded by Ayers et al.,[70] the softness kernel has been pointed out as nearsighted, while ω(r,r′) is not. The difference is demonstrated already in their integrals: ∫ω(r,r′)dr′ = 0, and s(r,r′) is intergrated to local softness ∫s(r,r′)dr′ = s(r). This substantiates using local approximation s(r,r′) = s(r)δ(r – r′) for the softness kernel; it has been well-established and explored in many computational procedures aiming at obtaining the Fukui function index.[37,45,46,71] Using this second equality in eq , the straightforward route was opened, leading to the local softness s(r) = [∂ρ(r)/∂μ]ν and hence also to the much desired Fukui Function f(r): the density gradient is approximately expressed by a simple vector equation

Equation represents an approximation to the exact second equality in eq . The new working formula for local softness reads

Both the electric field from all nuclei (r) and the density gradient ∇ρ(r) are readily computable; thus, we suggest eq to be the operating formula for the local softness index s(r). The solution for global softness is also within reach

Once S is found, the solution for the Fukui function is

The gradient theorem allowed the efficient solution to the problem of finding the Fukui function outlined by Chattaraj et al.[36]

Application of the abovementioned scheme to orbital types for a hydrogen atom and the hydrogenic ions has been presented first analytically in order to assess the value of the results obtained at the level of the local approximation. Numerical computations for atoms and ions follow. The results will be confronted with existing experimental data of global hardness for atoms and the results for the Fukui function calculated by other methods, before the method is adapted to monitoring reactivity in molecules, the real challenge in chemistry.

Equation reads that the electron density gradient is approximately proportional to the total electric field. It reflects direct nearsightedness[68] and results from Vela Gazquez local approximation.[37] The site, within the molecular space, where the density gradient is not parallel to the electric field constitutes the region where the Vela Gazquez approximation breaks down, and the deviation of the direction of the density gradient from the direction of the electric field seems to be a measure for the correlation neglected by this approximation. The practical solution to the application of eq to molecules is to introduce a softness tensor, and the local softness function will be obtained as a trace of this tensor. Additionally, performing the integral in eq for molecules requires regularization of the points where the electric field tends to be zero.

The novel theoretical approach to the Fukui function outlined above must undergo the test of its coherence with the established properties of chemical objects; atoms have typically been explored as the first testing ground for computational methods proposed for these indices in the past.[40−46] The unique advantage of this present method is in the a priori reproduction of the global softness (inverse global hardness). The rich literature focusing on reproducing the hardness of atoms[30] provides an opportunity to anchor the results of this work on the basic properties of atoms, before the method could be reasonably applied to describing the site selectivity in molecules.[47]

Fukui Function and Softness Indices for the Hydrogenic Orbitals

This procedure of obtaining the Fukui function from the local softness function (eqs , 8a, and 8b) can be applied to any molecular system. However, the simplest application is the spherically symmetric system—the atom. In fact, the much appealing simplicity of eq for local softness prompts us to test the result for orbitals in a hydrogen atom and hydrogenic ions. Should the local softness s(r) be a local function, eq implies that the gradient vector is parallel to the electric field vector in every point. This is the case for isolated atoms and ions; ∇ρ(r) and ε(r) vectors are parallel according to the spherical symmetry of the entities. The radial distribution of local softness for the 1s orbital is simply

The global softness by integration is . Radial distribution for the Fukui function for the 1s orbital in any hydrogenic species (eq ) is properly integrated to unity for any Z value

Results for the 1s orbital in eq are quantitatively close to those demonstrated by other methods of electron density analysis in a more sophisticated manner.[37,39,45,46] Analytical results for the local softness of other orbitals have also been calculated and are presented in the Appendix.

Radial distributions of local softness s(r) for canonical hydrogenic orbitals are demonstrated in Figure for a hydrogen atom (Z = 1). The global softness index, obtained by integration of the s(r) functions for orbitals, varies proportionally with the inverse of the square of the atomic number, as demonstrated in Table . These results allowed for calculation of the radial distribution for the Fukui function themselves for each hydrogenic type of orbitals, Figure . The abovementioned analysis proves that our procedure for local softness proposed in eq reflects the properties properly correlated to the electron density for orbitals in hydrogenic ions.

Figure 1

Radial distribution function of the local softness s(r) [in a.u] for orbitals in a hydrogen atom. (A): Orbitals 1s, 2s, and 3s; (B): orbitals 2p, 3p, and 3d.

Table 1

Collection of the Analytical Results from Integration of the Radial Distribution of Local Softness (cf. Appendix) in the Hydrogen-Like Cations: The Global Orbital Softness S [in a.u.] for the Hydrogenic Orbitalsa

orbital	1s	2s	2p	3s	3p	3d
S [a.u.]

Z is the atomic number.

Figure 2

Radial distribution of the Fukui function f(r) in [a.u.] for orbitals in a hydrogen atom (Z = 1, S = 6). (A): Orbitals 1s, 2s, and 3s; (B): orbitals 2p, 3p, and 3d.

Results for Atoms and Ions

The local softness index (eq ) has been calculated for 119 atoms and ions of 36 elements from the first to the fourth periods of the periodic table. The global softness has been obtained by the numerical integration according to eq . The radial distribution of the Fukui functions for atoms has been confronted with the commonly used density-based approximations of the finite-difference type: f(r), f–(r), and f0(r). The entire body of the results has been presented graphically in the Supporting Information.

Computational Methods

Numerical analysis was executed with the Gaussian 16 code.[72] The B3LYP method has been chosen, following its former successful tests in the electron density calculations in atoms.[73−75] The aug-cc-pvqz basis set has been routinely applied, except for the Ca atom (cc-pvqz). Using the pVTZ basis set recommended by Sadlej[76] was necessary for the potassium atom. The ground states of atoms have been identified by the specification of atomic electronic terms.[77]

Two calculation procedures have been applied for the local softness values: pseudoanalytical and numerical methods. The first one has been applied for a group of small atoms and univalent ions with 1–4 electrons and orbitals 1s and 2s only (H ÷ B+). The analytical function of the electron density has been formed with the use of the wave function coefficients together with basis set functions constructed from primitive Gaussians. Numerical results are obtained with analytical operations on the analytical form of the electron density function. The analytical density function has been integrated analytically to reproduce the expected number of electrons. The obtained accuracy was about 10–8 a.u. Analytical derivation of this density led to the electron density gradient and then to local softness by eq . The global softness values resulting from analytical integration of local softness for a hydrogen atom and He+ ions were 6.1209 and 1.5156 a.u, respectively, thus exceeding the values expected from the analysis of orbital softness. This is due to the fact that Gaussian-type orbitals do not reproduce correctly the electron density of the H atom.

For practical reasons, the numerical DFT method has been used in order to form the electron density and electron density gradient to overview the local softness of atoms and ions in rows 1–4 of the periodic table. The effect of degeneracy of frontier orbitals has been avoided using the integral electron density for every atom and ion. As it has been proved by Kohn, the ground-state density is unique even for systems where frontier orbitals are degenerate.[78,79]

A simple practical method has been applied to circumvent another well-known difficulty in reproducing spherical symmetry of the electron density in atoms. Two steps have been involved in this procedure.

Integration of the electron density has been routinely made by the spherical algorithm, 4πr2ρ(r), for a density ρ(r) variable in one direction only. The result of integration to the proper number of electrons served as evidence of a spherical symmetry of the density. The accuracy of this test has been reported below, separately for atoms grouped by periods.

For atoms whose basic electronic terms were other than the S type, the numerical averaging procedure over the principal coordinate axis and all diagonal directions has been applied to the raw computational results. The subsequent integration (i) provided proof for the sufficient quality (symmetry) of the averaged density for the purpose of this study.

The density gradient has only been calculated in one direction from the spatially averaged numerical density data.

The integral electron density has been obtained with the grid produced within the Gaussian package. The grid was fine, and the accuracy was controlled by the integration of the resulting radial distribution of the density to the proper number of electrons. The distance between the calculated density points for atoms/ions in the first and second periods was 0.05a0, and the radius of the grid was 15a0. For the third and fourth period, the grid was fine by 0.02a0 within the 20a0 of the grid radius. Standard deviations between the nominal number of electrons and the result of integration of the density function determined for the main group elements in the first rows of the periodic table (1 + 2), 3, and 4 were 0.009, 0.055, and 0.024, respectively. Standard deviation was considerably higher for the group of 3d elements in period 4 (0.42); however, given the large number of electrons in these atoms, the relative accuracy of the procedure was at the same level for all groups of atoms (ca. 0.1–0.5%). The global hardness for a hydrogen atom served as a test for the efficiency of the method and the basis set, leading to SH = 5.83 a.u., reasonably close to the analytic result of 6.0 a.u. (Table ).

Local Softness of Atoms

Results for the local softness of atoms have been presented as the spherical radial distribution s(r) = 4πr2s(r). Collections of the local softness for atoms in the second and third row of the periodic table are shown in Figure together with the results of integration of s(r) curves (the global softness).

Figure 3

Calculated radial distribution function of local softness s(r) in [a.u.] for atoms in the second row (A) and the third row (B) of the periodic table. Calculated global softness for atoms has been inserted: S [a.u.].

An example of meaningful, alternative presentation of the results is provided in Figure as the Fukui functions for a group of ions isoelectronic with neon atom (N = 10) and argon atom (N = 18) electronic configurations. Since the Fukui function integrates to unity, the separate roles of global hardness and the Fukui function itself in describing properties of ions may be appreciated.

Figure 4

Calculated radial distribution of Fukui function f(r) in [a.u.]: (A) for the neon atom and the isoelectronic ions and (B) for the argon atom and the isoelectronic ions. Calculated global softness has been inserted: S [a.u.].

Global Softness of Atoms and Ions

The global softness S resulting from the integration of the radial distribution of the local softness s(r) for 36 atoms is presented in Figure as a function of the atomic number Z. The diagram is clearly divided into sections, corresponding to the valence electron type of atoms: 1s, 2s, 2p, 3s, 3p, 4s, 3d, and 4p. These results are parallel to the analytic results for the hydrogenic orbitals presented in Section . Analogous diagrams for univalent cations and anions are presented in Figure ; several bivalent anions and cations have also been marked in these diagrams.

Figure 5

Global softness of atoms S [a.u.] calculated by integration of the radial distribution of softness density s(r) in atoms.

Figure 6

Global softness S [a.u.] calculated for univalent cations and for univalent anions of elements in rows 1–4 of the periodic table. Selected (2+) cations and (2−) anions have also been marked. The scale of the ordinate axis has been expanded by 2:1 (for cations) and contracted by 1:2 (for anions) as compared to Figure (for atoms).

Calculation of the global softness index has opened access to numerical data for global hardness of all atoms and ions selected for this study, η=1/S. This result has been confronted with the measure of absolute hardness originally proposed by Parr and Pearson,[7], and subsequently adjusted by dropping out the arbitrary coefficient,[13] η = I – A. This formula has been recognized as the working definition of hardness. From our present point of view, this is only an approximation which originates from finite difference approximation. Therefore, we tested correlation between the inverse of softness computed with our procedure with I – A to establish the appropriate coefficient: . Following the finding of nonidentical global softness for various types of orbitals (Table ), the correlation has been tested separately for atoms in rows of the periodic table. Analogous correlation has been investigated for univalent cations, as the higher ionization potentials are known with sufficient accuracy (I2); the correlation has been assumed accordingly: The low accuracy of electron affinity data would not allow for an analogous test for anions. The resulting correlation parameters are presented in Table .

Table 2

Numerical Correlation Parameters in the Relationship between the Conventional Measure of Hardness and the Inverse Softness Resulting from This Work for Atoms and Cations in Rows 1–4a

row no in the periodic table	valence orbitals	αa	αc
1 + 2	s,p	2.10	2.10
2	s,p	1.70	1.63
3	s,p	1.02	1.08
4	s,p	0.71	0.94
4	s,p,d	0.86	0.85

The entries have been separated into valence electron sections. The noble gas atoms (He, Ne, Ar, and Kr) and the group IA cations have been omitted (see text). The most recent values for first and second ionization energies (I, I2) and the electron affinity data (A) data are from refs (80) and (81), respectively, αa: atoms, αc: cations.

Another test for the global softness values resulting from our new approach was provided by the analysis of atomic/ionic radii. The linear relation between the atomic hardness and inverse atomic radii has first been suggested by Gazquez and Ortiz.[82] The concept has been developed by Komorowski within the chemical approximation;[83,84] the van der Waals radii have been found to correlate with the inverse global hardness (I – A)−1. Consequently, the correlations between the calculated global softness for atoms and cations and the selected measures of their radii have been sought. The correlations have been limited to atoms grouped by periods, and three reliable sources for atomic radii have been explored (Table ).

Table 3

Correlation Parameters for the Relation between the Atomic Radii and Calculated Atomic Softness S Assumed as rat= β S (rat and S in a.u.)a

row no	β1	β2	β3
2	0.340	0.369	0.424
3	0.620	0.621	0.554
4	0.828	0.734	0.658

Covalent radii are from the Cambridge database (β1).[85] Two types of van der Waals radii are from the recent analysis by Batsanov (β2: crystal, β3: equilibrium).[86]

Fukui Function Index

Results for the Fukui function indices have been calculated from the local and global softness f(r) = s(r)/S and have been presented separately for each atom as the corresponding radial distribution function f(r) = 4πr2f(r). For the sake of clarity in comparing the results, f(r) has been presented jointly with that for the corresponding cations and anions. An inspection of the Fukui function in isoelectronic atoms and ions has also been provided in a separate diagram for each atom. Finally, the comparative diagrams for the Fukui function with the normalized radial density ρ(r)/N and with the widely used approximation (eq –11c) have been demonstrated[87]

The combination of the abovementioned data in the joint compact picture for the carbon atom is shown in Figure . Analogous diagrams for all atoms under study have been shown in the collection given as the pdf presentation attached in the Supporting Information. Using the spherical integral density for atoms (eq –11c) allows for avoiding the degeneracy of (frontier) orbitals, which must be taken into account for molecules, when only ρLUMO(r) and ρHOMO(r) are considered in an approximate version of eq –11c.[88,89]

Figure 7

Radial distribution for the Fukui functions f(r) calculated for the carbon atom and its ions (A) and for the corresponding isoelectronic ions (B). The global softness values S in [a.u.] for atoms and ions have been inserted accordingly (A,B). Comparison to the typically used approximations for the Fukui function (eq –11c) and the normalized electron density function ρ(r)/N are shown in figures (C,D). Analogous diagrams for all 36 atoms are available in the presentation enclosed as the Supporting Information.

Discussion

Three reactivity indicators quantified by the present method are closely related: the local softness s(r), the global softness S, and the Fukui function f(r). The Fukui function is designed to characterize how much the electron density changes locally with the global change of the number of electrons. On the other hand, global softness gives the change of the total number of electrons with the change of the external potential. The global softness is the integral of the local softness over the whole space (eq ), and the local softness is a product of the global softness and the Fukui function (eq ). The predictive roles of the Fukui function and local softness are remarkably different. Although f(r) and s(r) provide basically the same information on the relative site reactivity within a molecular system, s(r), due to the incorporated information on the value of the global softness,[12] is more suited for intermolecular reactivity since the electron flow is inversely proportional to the sum of the total hardness of the reacting objects[9] (and global softness is the inverse of hardness). The second remark is crucial: the ultimate goal for indexing atoms is in describing their reactions and hence their interaction with other atoms (in molecules). Thus, the local softness function provides complete information on molecular reactivity that can be decoded from the electron density function.

The arguments indicating the role of local softness as the leading parameter in indexing atomic reactivity follow directly from the definitions of the Fukui functionand local softness

The first equality in eq describes a response of the density to ionization in the rigid molecule (constant external potential—rigid geometry constraint), while the second one refers to the polarization of a closed system (constant number of electrons). Both effects are internal in their character; they describe the response of a system to a change introduced by ΔN or Δν(r). Here, some important distinction follows from the grand canonical ensemble approach.[18] (Valuable considerations on this topic have also been presented by Malek and Balawender.[90]) Chemical potential in eq is the electronic chemical potential often denoted as μe. On the other hand, the variational derivative in the first equality of eq is over the chemical potential of the surrounding reservoir, denoted as μbath. In the zero-temperature limit when the entropic term vanishes, μbath = μe. Our discussion within this paper is constrained to this zero-temperature limit.

Equation provides deep insights into the description of chemical reactivity. The first equality is the response of the density to a change that may not be as evident as ionization: changing the chemical potential by Δμ at constant external potential (hence the position of all nuclei). By Hohenberg—Kohn theorem[24]where F[ρ] denotes the universal Hohenberg–Kohn functional that includes kinetic energy functional plus electron–electron repulsion. Thus, the change of μ with the constant v(r) is equivalent to the change of the functional derivative of the universal functional often denoted by u(r) = δF/δρ. This means that local softness describes the change of electron density with the variation of the derivative of the universal functional and on this way, due to the change of the slope, of the universal interactions. The second equality of eq proves that this is the essence of reactivity. It describes an effect of the in/out flow of electrons between a system (e.g., an atom) and an external reservoir (e.g., other parts of a molecule), as to keep the chemical potential constant, when compensating the externally induced change Δν(r). The local softness is considerably richer as a potential source of information on reacting molecules.

The maxima of the local softness for hydrogenic orbitals (Figure ) indicate an interesting sequence: their heights and their distances from origin increase with the principal quantum number (n: 1, 2, and 3) and decrease with an increase in the orbital quantum number (l: s, p, d); this effect is reflected by the global softness of orbitals, as given in Table . Remarkably, the uniform relation S ∝ σ/Z2 between the atomic number and the global softness of hydrogenic orbitals has been proved, with σ constants characteristic for the orbital type. This property is of great importance for understanding the sequence in global softness of atoms in the periodic table.

Radial distributions of the Fukui function indices show another effect (Figure ): the heights of their maxima decrease with the main quantum number of an orbital. This supports the conclusion on the more informative character of the local softness over the Fukui function: large and distant maxima of s(r) indicate the regions where electron exchange might occur (Figure ), while the strongest maximum for the 1s orbital Fukui function (Figure ) does not seem to provide any relevant chemical information, given the internal character of this orbital. Small regions of the negative Fukui function (and local softness) for orbitals are not unexpected; the phenomenon has been discussed and theoretically explained in early DFT studies.[55,91−93] Actual values for the Fukui function calculated for atoms and ions by the presented method are all positive (eq –11c).

The global softness in ionized species and/or isoelectronic atoms/ions allows for observation of the effect of electronic interaction on softness; some instructive examples are presented in Table . Examples of atoms and ions in Table contain the 1s and 2s orbitals only. Softness S calculated by the pseudoanalytical method (cf. Section ) is systematically higher than the results from the practical numerical method by ca. 8% on average. Reasonable consistency between the global softness and a sum of orbital softness for electrons is found for one electron species (H and He+) and for some two-electron systems: He, Li+, and Be2+. Softness S for all anions H–, He–, and Li– significantly exceeds the expectation from summation of the orbital softness. Possible chemical explanation of this regularity is the role of electron repulsion; the sum of orbital softness values is calculated as if the contributions from each electron were simply additive. The same effect might be responsible for the regular tendency observed in 1s 2s atoms: softness S (numerical) is systematically larger than the sum of orbital contributions to global softness, and the effect is strong for Z = 1 and vanishes for Z = 4. The role of electron interactions is supported theoretically: for a set of isoelectronic atoms, the application of the virial theorem allows for estimation of the global hardness (inverse global softness) as η = 1/S ≅ Δ(J[ρ] + Exc[ρ])/ΔN. The global hardness provides a measure for a sensitivity of the electron–electron interactions to ΔN. This would explain the data in Table for anions with the highest softness (lowest hardness): their sensitivity to changing an electron number is indeed extra-low; the exceedingly high results of softness for anions by this method reflect the well-known weakness of this measure of softness for anions (last column in Table ), as compared to the method present hereby, equally applicable to atoms and any ions. Direct confrontation of the calculated softness data with the absolute hardness[8] adjusted by the improved correlation coefficient (Table ) for neutral atoms and cations proves how these experimental measures of softness match the sequence of softness S calculated by integration of s(r) within each group.

Table 4

Global Softness S Calculated for Isoelectronic Atoms and Ions by Three Methods: (i) S by Integration of the Radial Distribution of Local Softness s(r) by the Pseudoanalytical Method (Section ); (ii) Sum of Contributions from Each Electron and Orbital to the Global Softness (Table ); and (iii) I – A Values (refs (82) and (83)) Adjusted by the Correlation Parameter α = 2.10 (Table )

atom/ion	atomic number	number of electrons	orbitals (occupied)	S [a.u.] (pseudo-analytical)	S [a.u.] (numerical)	sum of orbital softness’ [a.u.]	S = α(I – A)−1 [a.u.]
H	1	1	1s	6.1209	5.83	6.00	4.45
He+	2	1	1s	1.5156	1.38	1.50	1.92
H–	1	2	1s2	20.5648	19.91	12.00	75.6
He	2	2	1s2	3.7635	3.48	3.00	2.28
Li+	3	2	1s2	1.5402	1.35	1.34	0.81
Be2+	4	2	1s2	0.8338	0.70	0.76	0.42
He–	2	3	1s22s	13.7064	12.45	9.00
Li	3	3	1s22s	6.5732	6.29	4.02	11.97
Be+	4	3	1s22s	3.0830	2.88	2.26	6.42
Li–	3	4	1s22s2	16.2529	14.09	6.67	92.4
Be	4	4	1s22s2	6.0468	5.77	3.75	5.82
B+	5	4	1s22s2	3.3763	3.16	2.40	3.39

Results for local and global softness for atoms in various periods of the periodic table of elements are much instructive (Figure ); clear systematic changes have been demonstrated. Minimum softness values for noble gas elements are in accordance with the maximum hardness principle.[94] As the elements in the second row are concerned, the calculated global softness is highest for Li and then falls to the lowest value for Ne, spanning the range of 6–3 a.u. Global softness in the third period is less diverse, 3.3–3.92 a.u., with the maximum at Al (Figure B). Following the analysis of orbital softness values presented above, this is reasonable for an atom with the sole 3p electron. An analogous effect is not observed in boron (Figure A), possibly due to the large difference of softness between neighbors: Be—5.77, B—5.15, and C—4.50.

The wide tail in the long distance from the nucleus (reaching as far as 10 a.u.) is responsible for the generally soft character of alkali metals (Figure ). When the electronic shell is filled with more electrons, this immediately affects softness: it rises systematically, becoming the dominant one; the internal shell peak is washed out when reaching the next noble gas configuration. This suggests that the valence shell is the one responsible for the observed effect of softness.

The calculated local softness of atoms allows some statements due to its behavior in isoelectronic entities—the noble gas configuration has been selected for this purpose (Figure ). With the number of electrons being constant, variation of local softness with the atomic number is exposed. Since only the stable species have been shown in Figure , building the wide tail of local softness in bivalent anions is an interesting observation. This effect is responsible for the sequence of global softness: S (O2–) > S (S2–). The dependence on Z and N for hardness has been first discussed by March;[95] the problem deserves to be revisited once the local and global softness had been endowed with numerical measures rooted in the electron density.

The overview of global softness for atoms presented in Figure has been supplemented with the marking of valence orbitals of each atom. The irregularities in S(Z) dependence are readily explained by the subsequent filling of the subshells in atoms; the general trend of decreasing softness with increasing atomic number illustrates the regular S ∝ Z–2 function for orbitals. Analogous data for cations and anions are presented in Figure ; global softness for anions is on average twice as large as that for atoms, while for cations, it is half as small; the ordinate scales on Figure have been adjusted accordingly.

The picture in Figure is coherent with the result for orbital softness and provides a valuable hint to the search for meaningful correlations between the new measure of global softness and the traditional measure of hardness introduced by Pearson and Parr,[8] with α = 2, leading to a barely satisfactory correlation. By limiting the correlation to periods (Table ), according to the finding for the variable global orbital softness (Section ), excellent linear correlations between (I – A) and hardness (inverse softness,η = 1/S) have been found with specific coefficients for each row. α = 2.10 holds roughly for atoms and cations in periods 1 and 2; better correlation is obtained for elements in period 2 only with average α = 1.66; in periods 3 and 4, the average coefficients for all atoms and cations are 1.05 and 0.86, respectively, with somewhat lower accuracy. Chemical intuition behind the concept of absolute hardness has been corroborated.[8] The reason why a general correlation should not be insisted on has now been unveiled in the discussion of orbital softness (Section ).

The proportional relation between the new measure of softness’ for atoms/ions and their diameters documented in Table confirms the early results in this matter. Various measures for the diameters have been tested in the past, including the atomic refractions understood as volumes of bonded atoms.[83,84,96] The linear correlations between atomic radii from three various sources strongly support calculated atomic softness S as the proper measures of this atomic property.

The quantified order of softness for ions has been originally pursued by Pearson.[30] The author introduced an arbitrary empirical hardness parameter Δ for this purpose; Δ is the difference of the dissociation energies of bonds to the chosen standards: F– and I– for the group of cationic acids and H+ and CH3+ for anionic bases . Results of this present work can be compared to those proposed by Pearson for a handful of ions as a sequence of increasing softness for the corresponding ions.

The sequence of increasing softness for alkali metal cations by the two measures coincides quite precisely: Li+ < Na+ < K+. Among the few atomic anions listed in Pearson’s book,[30] H– is the softest one as it is in this work. However, the sequence of increasing softness for halogen anions by the two methods is contradictory: F– < Cl– < Br– by Pearson and Br– < Cl– < F– by this work. By cross-checking this with electron affinity data (η ∝ A), the absolute hardness values are found to be rather close for all three anions, pointing out the chloride anion being only slightly harder than fluoride and bromide.[81] When the absolute hardness data (A) are combined with the correlation coefficients determined separately for rows 2, 3, and 4 in the periodic table (Table ), the order of global softness by this work has been confirmed.

A possible source of the discrepancy may be the factors determining the dissociation energy of bonds that served as chosen standards for anions in Pearson’s concept. The energy change due to the charge transfer is dependent both on electronegativity difference (Δχ) and the sum of hardness (η).[8,9]

Electronegativity and hardness of both reference cations H+ and CH3+ chosen as standards for anionic bases are considerably higher than electronegativity and hardness of anions; thus, by the conventional measure (eq ), the dissociation energies and for anionic bases will be dominated by the properties of the standards, masking the effect of very low hardness (and electronegativity) of changing partners in Δ. The Δ measure for cationic acids will not be affected since the chosen standards for this group (F– and Br–) are of much lower electronegativity (and hardness) than for the cations; hence, variable properties of cationic acids can be exposed.

Due to the high interest in the Fukui function index over the years, the complete overview of Fukui functions for atoms (and many of their ions) in 1–4 periods has been composed and attached to this paper as the presentation in the pdf format (Supporting Information). An example of the data available in this collection is shown in Figure for carbon. Section A gives the Fukui function together with global softness for an atom and its ions. By the general observation, the Fukui function (normalized to unity) shows a maximum for the cation, atom, and anion in the same place (distance to origin), with the height of the peak decreasing in the sequence of increasing global softness. The same is true for the Fukui functions of isoelectronic ions of an atom shown in Section B. Section C contains the new Fukui function together with the traditional ones (eqs and 11b); similarly, in Section D, the new Fukui function is confronted with the one given by eq and also with the normalized electron density function ρ(r)/N.

This comparison clearly shows the specialized and approximate character of the traditional Fukui functions (eq –11c), divided into characteristic measures for nucleophilic, electrophilic, and radical attacks separately. Typically, the new Fukui function for an atom is close to the one dominating among the traditional ones (period 2), hence providing the characteristics of an atom as such, with no need for a preliminary specification of the type of its reactivity. For atoms in periods 3 and 4, the new Fukui functions expose the role of internal shells that is typically washed out by the calculation method of the traditional FF (eq –11c). By the opinion of the authors, the future in chemical applications belongs more likely to the now available local softness, rather than to Fukui functions.

Conclusions

The key finding of this present work is eq and its consequences. The electron density gradient is related to the total electric field from nuclei of the molecule. For atoms, local softness is the coefficient in this linear dependence. For systems of lower symmetry, the local approximation for the softness kernel may not be held at the points of vanishing electric field. In such cases, one may need to introduce the softness tensor,[4] and the solution of eq will only provide its trace.

The interchange between the derivative over μ (local softness) and the derivative over r (gradient) in the analysis of electron density opens the route to evaluate the reactivity indices at the cDFT level. This could not have been achieved without the previous analysis leading to the electron density gradient theorem (eq ) and the nearsightedness property of the electronic matter as demonstrated by others.[68−70]

The working formula that leads to local softness (eq ) represents the rational approximation that allows for tractable computations of local softness. By application of the abovementioned finding to atoms, the authors demonstrated its usefulness in opening facile access to otherwise unavailable quantitative data of softness (local and global) for orbitals and atoms alike. The results properly reproduce the spherical symmetry of free atoms and ions, also providing unique quantitative information on the softness of individual orbitals, the key viewpoint in contemporary discussions on reactivity. The choice of atoms as the preliminary test group for the method of obtaining the Fukui function has been a standard approach.[37,39,40,45,46]

For atoms, ∇ρ(r) and (r)vectors are parallel according to the spherical symmetry of the entities. For molecules, there is no such topological constraint, and application of our method to molecules will require some modifications: the trace of the softness tensor provides the local softness function and the analytical integration will have to be replaced with numerical integration techniques. However, the idea remains the same and the numerical calculations for both global and local softness are feasible and will provide the Fukui function. The method also opens a new perspective: the electric field is clearly attributed to atoms in molecules and so may be the local softness and the Fukui function.

Results for the local and global softness and the Fukui functions demonstrated for orbitals and atoms indicate the potential for further use of the rational approximation in eqs and 7. By generalization of the method, higher derivatives of the electron density may be calculated, for example, the dual descriptor f(2)(r) = df(r)/dN,[97] the third derivative of energy over N (γ, hyperhardness[62,98,99]), and also the derivatives over Z.[95] The method is also applicable for molecules, thus opening a field for exploration of the potential measures of properties for atoms in molecules by the DFT indices that identify atoms, with no need for spatial divisions of the density into atoms: the indices rooted in the Hellmann–Feynman force. Their properties have already been described, but no appropriate computational tools have yet been elaborated for practical purposes: the nuclear reactivity ,[56] the nuclear stiffness ,[61] and the softening index for atoms and/or bonds .[2] The access to these potentially attractive tools for characterization of reactivity of atoms in molecules has been opened with the working procedure to the d/dN derivative of electron density as presented in this work.