A first (local) bridge between Kohn–Sham density functional theory and the quantum theory of atoms in molecules of Bader is built by means of a second order reduced density gradient expansion of the exchange-correlation energy density at a given bond critical point. This approach leads to the definition of new “mixed” descriptors that are particularly useful for the classification of the chemical interactions for which the traditional atoms in molecules characterization reveals insufficient, as for instance the distinction between hydrogen and agostic bonds.

## I. INTRODUCTION

Although all chemical information is in principle included in the Schrödinger equation and in the derived electron density (and wave function), chemists however still aim to develop interpretative tools that, even if they may contain less information, are very useful to straightforwardly account for (or even predict) a certain number of relevant physico-chemical properties. These tools can be either rooted in a mere pragmatic approach or in more fundamental theories. There is nevertheless no reason that these two frameworks should be separate. On the contrary, making them converse often constitutes a fruitful opportunity to validate and extend the experimental observations or a given theory, as epitomized by the conceptual density functional theory (DFT) that enables to physically ground and enrich important empirical chemical concepts and descriptors, such as hardness and electronegativity.^{1}

This appears all the most meaningful that DFT (Ref. 2) has become a choice tool in computational chemistry. As well known, this theory rests on the primary physical observable, the electron density, an ingredient that it shares with the quantum theory of atoms in molecules (QTAIM) developed by Bader and coworkers.^{3,4} QTAIM notably offers a powerful tool of analysis and rationalization, based on the density gradient field, this latter quantity stemming from theoretical calculations or from experimental analyses.

However, despite this common aspect, only few direct links have been evidenced between the two, mainly based on Thomas–Fermi related models^{5,6} or/and in a Bohmian mechanics perspective.^{5–7} Similarly, the use of a combined conceptual DFT/QTAIM approach has not been often exploited.^{8,9} Our purpose here is however slightly different because we will focus on Kohn–Sham DFT (KS-DFT) (certainly the most popular formulation of DFT), and we will give some hints for building a first (from the best of our knowledge) direct (local) bridge between KS-DFT and QTAIM.

## II. NEW LOCAL QTAIM DESCRIPTORS BASED ON EXACT PROPERTIES OF THE KOHN-SHAM EXCHANGE-CORRELATION FUNCTIONAL

In this communication, we will limit ourselves to a local point of view of QTAIM, focusing on the so-called bond critical points (BCPs) where the density gradient (but not the density itself) vanishes. The local properties at these points are in fact often used to discuss the nature of bonding,^{10} and we intend now evaluating the behavior of some KS-DFT quantities at these points.

To this end, within a sufficiently small sphere near a BCP, the KS exchange and correlation energy densities $ex,c$ can be supposed to be very close to their exact second order expansion for weak inhomogeneous electron densities,^{11–15}

where $s(r\u20d7)=As\Vert \u2207\u20d7\rho (r\u20d7)\Vert /\rho (r\u20d7)4/3$, $t(r\u20d7)=At\Vert \u2207\u20d7\rho (r\u20d7)\Vert /\rho (r\u20d7)7/6$, $As\u22121=2(3\pi 2)1/3$, $At=(\pi /3)1/6/4$, $\mu x=10/81$, $\mu c\u22480.067$.

Such an use of a truncated second order expansion is also a reminiscent of that used for the derivation of the electron localization function (ELF) (the expanded quantity being the spherically average conditional pair probability)^{16,17} and of the one used for the study of deformation of materials,^{18} the expanded quantity being the stress tensor.^{19,20}

Let now $\epsilon u\u20d7$ be an infinitesimal vector, with $\Vert u\u20d7\Vert =1$ defined by the azimuthal $\theta $ and zenith $\varphi $ angles (following Zwillinger’s convention^{21}). In the vicinity of a BCP ($H\u0347$ denoting the density Hessian and $\lambda i$ its eigenvalues),

For a local physical property $P$, we can define the variation

In the case $\delta P$ is of $\epsilon $ second order, we make the angular and the distance dependences disappear by performing the following spherical average,^{22}

If we now specify $P$ to be the previously mentioned DFT-related quantities, one gets ($Cx=\u2212(3/4)(3/\pi )1/3$, $Cc=1$, $[\u2009]\u2032\u2261d/d\rho $),

The structure of these equations is thus the same for exchange and correlation. The variation of the purely local energy density (due to the density inhomogeneity) is multiplied by the laplacian, whereas the variation stemming from the gradient inhomogeneity (controlled by $\mu x,c$) is multiplied by ◇. Defining

one obtains

Equation (10) is exact (within the limits of starting hypothesis) since only properties known to be respected by the exact KS exchange-correlation functional have been used. Obviously, the numerical applications will use approximate (and so inexact) electron densities.

Complementary, the following ratio (that locally compares the exchange and the correlation energy densities up to the second reduced gradient orders) can be evaluated,

## III. APPLICATION TO SELECTED TEST CASES

In order to have a preliminary analysis of the interest of these new local descriptors, especially for discriminating bond families, we will focus on hydrogen and agostic interactions, whose nature is still controversial but whose importance in biochemistry and in organometallics is fundamental.^{23} For these bonds, the intervals of variations for $\rho (r\u20d7c)$ and $\u22072\rho (r\u20d7c)$ are thought to be entirely separated, suggesting the use of these descriptors as discriminating indexes.^{24,25} We have recently revised this last statement by a systematic study on 23 bis(imino)pyridyl complexes for a total of 40 agostic cases,^{26} leading to a partial reconsideration of this dichotomy. In order to illustrate this point and to motivate the induced need for new local descriptors, let us extract three molecules from this database (Fig. 1). A titanium(IV), a vanadium(III) and a chromium(II) complexes with $\beta $ (Ti, V) and $\delta $ (Ti, Cr) agostic bonds, to which we will add a nickel(II) complex featuring an agostic $H\delta $ that is nearly anagostic (or preagostic) $(\u2220C\delta H\delta Ni\u2248124\xb0)$ according to Yao’s and al.’s nomenclature.^{27,28} On the other hand, four hydrogen-bonded systems will be considered, two (H1 and H2) from the formic acid dimer family, and the water and HCl dimers (Fig. 1). As shown in Table I and Fig. 2,^{29} the $(\rho (r\u20d7c),\u22072\rho (r\u20d7c))$ pair is not sufficient to distinguish the agostic from the H bonds, that is to say one cannot unambiguously draw in this plane continuous zones corresponding to the two bond types.

. | Ti(IV) . | V(III) . | Cr(II) . | Ni(II) . | H1 . | H2^{a}
. | H3 . | H4 . | ||
---|---|---|---|---|---|---|---|---|---|---|

$Ti\ue5f8H\beta \u2002ago$ . | $Ti\ue5f8H\delta \u2002ago$ . | $V\ue5f8H\beta \u2002ago$ . | $Cr\ue5f8H\delta \u2002ago$ . | $Ni\ue5f8H\delta \u2002ago$ . | $O\u2026H$ . | $O\u2026H1$ . | $O\u2026H2$ . | $O\u2026H$ . | $Cl\u2026H$ . | |

$\rho (r\u20d7c)$ | 0.040 | 0.030 | 0.047 | 0.030 | 0.014 | 0.053 | 0.042 | 0.050 | 0.028 | 0.013 |

$\u22072\rho (r\u20d7c)$ | 0.141 | 0.114 | 0.202 | 0.120 | 0.046 | 0.147 | 0.141 | 0.168 | 0.083 | 0.036 |

$\u25c7\rho (r\u20d7c)$ | 0.042 | 0.029 | 0.072 | 0.028 | 0.004 | 0.136 | 0.092 | 0.146 | 0.028 | 0.004 |

$Qxc0\u20132(r\u20d7c)$ | −3.297 | −2.500 | −3.672 | −2.604 | −1.946 | −1.747 | −1.684 | −1.686 | −1.698 | −1.519 |

$Pxc0\u20132(r\u20d7c)$ | 5.399 | 5.134 | 5.567 | 5.120 | 4.490 | 5.698 | 5.443 | 5.628 | 5.053 | 4.438 |

. | Ti(IV) . | V(III) . | Cr(II) . | Ni(II) . | H1 . | H2^{a}
. | H3 . | H4 . | ||
---|---|---|---|---|---|---|---|---|---|---|

$Ti\ue5f8H\beta \u2002ago$ . | $Ti\ue5f8H\delta \u2002ago$ . | $V\ue5f8H\beta \u2002ago$ . | $Cr\ue5f8H\delta \u2002ago$ . | $Ni\ue5f8H\delta \u2002ago$ . | $O\u2026H$ . | $O\u2026H1$ . | $O\u2026H2$ . | $O\u2026H$ . | $Cl\u2026H$ . | |

$\rho (r\u20d7c)$ | 0.040 | 0.030 | 0.047 | 0.030 | 0.014 | 0.053 | 0.042 | 0.050 | 0.028 | 0.013 |

$\u22072\rho (r\u20d7c)$ | 0.141 | 0.114 | 0.202 | 0.120 | 0.046 | 0.147 | 0.141 | 0.168 | 0.083 | 0.036 |

$\u25c7\rho (r\u20d7c)$ | 0.042 | 0.029 | 0.072 | 0.028 | 0.004 | 0.136 | 0.092 | 0.146 | 0.028 | 0.004 |

$Qxc0\u20132(r\u20d7c)$ | −3.297 | −2.500 | −3.672 | −2.604 | −1.946 | −1.747 | −1.684 | −1.686 | −1.698 | −1.519 |

$Pxc0\u20132(r\u20d7c)$ | 5.399 | 5.134 | 5.567 | 5.120 | 4.490 | 5.698 | 5.443 | 5.628 | 5.053 | 4.438 |

^{a}

Contrary to H1 for which the two H bonds are strictly equivalent ($C2h$ symmetry), these bonds are not equivalent in the H2 optimized geometry (symmetry group: $Cs$).

Such a failure is particularly evident for the titanium and the H2 complexes. However, while the laplacian values are confused, the eigenvalues of the density Hessian are quite different at the corresponding BCPs. Indeed, for the $O\u2026H1$ bond, $\lambda 1=\u22120.074$, $\lambda 2=\u22120.070$, $\lambda 3=0.286$, while for the titanium agostic one: $\lambda 1=\u22120.047$, $\lambda 2=\u22120.011$, $\lambda 3=0.200$. Yet, their sums are equal due to numerical compensations between positive and negative values, so that the laplacian cannot discriminate them.

As shown in Fig. 3, differentiation becomes possible when the $(Qxc0\u20132(r\u20d7c),Pxc0\u20132(r\u20d7c))$ representation is considered, two distinct and separate groups clearly appearing. This is notably due to the fact that the $Qxc0\u20132(r\u20d7c)$ values for the agostic bonds can be almost twice (in absolute value) those for the H bonds. Notably, the almost preagostic complex lies in an intermediate position in this representation, between pure agostic and H bonded systems.

However, exchange and correlation always vary following opposite trends, the variation rates of the exchange energy density being higher in absolute value, even if it is less pronounced for H bonds. Let us also recall that the ratio of the total (integrated) exchange energy over the total correlation one is about nine for atoms or molecules. Obviously, this result is a kind of spatial average, and such a ratio is not verified at any point. However, one can show that for the $C\ue5fbO$ bonds in H1 and H2, and the $Ti\ue5f8C\alpha $ bond, $Pxc0\u20132$ is close to this value. This is not any longer the case for the H and agostic bonds for which the weight of the correlation energy is larger.

Interestingly, it can be noticed that we have the same $\u25c7(r\u20d7c)$ and $\rho (r\u20d7c)$ values for the studied bonds in the H3 and the chromium complexes; but, as $Qxc0\u20132(r\u20d7c)$ also includes a laplacian contribution, our new descriptors succeed in our typology endeavor, proving that the laplacian and the diamond are complementary. In order to show how this works, we must compare the information carried by the corresponding Hessians and, in particular, their eigenvalues.

Nevertheless, as the eigenvectors $y\u20d7i$ for a BCP are not generally collinear to the eigenvectors for another BCP, the direct comparison of the $\lambda i$ has no direct meaning. So, in order not to lose information, the knowledge of the $\lambda i$ must be replaced by the one of three (independent) scalar combinations of them that are invariant by rotation. A natural choice would be to use the three symmetric functions that arise from the expansion of the Hessian characteristic polynomial: $\sigma 1=\u2211i\lambda i$, $\sigma 2=\u2211i<j\lambda i\lambda j$, and $\sigma 3=\u220fi\lambda i$. Nevertheless, one can equivalently privilege the $(\sigma 1\u2261\u22072,\u25c7)$ pair over $(\sigma 1,\sigma 2)$, these two similarity invariants being related by $\sigma 2=(\sigma 12\u2212\u25c7)/2$.

Such a choice is motivated by the fact that ◇ has a clear physical interpretation since $\delta \Vert \u22072\rho \Vert (\epsilon ,u\u20d7)r\u20d7c=\epsilon 2(\u2211i\lambda i2(u\u20d7.y\u20d7i)r\u20d7c2+o(1))$, so that: $\u27e8\delta \Vert \u22072\rho \Vert \u27e9r\u20d7c(2)=13\u25c7\rho (r\u20d7c)$. Thus, ◇ measures the variation rate of $\Vert \u2207\u20d7\rho \Vert 2$ near a given BCP. As ◇ only involves the eigenvalues squares, there cannot be compensations between negative and positive values; besides small eigenvalues will give a negligible contribution, so that the considered agostic and H bonds may have sufficiently distinguishable ◇ values.^{30}

## IV. CONCLUSIONS

In this communication, we have calculated the average variations of the exchange and the correlation energy densities in the vicinity of a BCP, thanks to second order reduced gradient expansions. This approach led to the definition of new descriptors that were then shown to be useful to classify bonds, since they provide crucial supplementary information to that provided by the usual $(\rho (r\u20d7c),\u22072(r\u20d7c))$ pair that can be not sufficient to univocally characterize chemical bonding. Obviously, extended systematic tests are in progress to confirm this outlined typology based on KS-DFT energetics. These results constitute, from our point of view, incentive arguments to foster a deeper study of the relationships that can exist between the related KS-DFT and QTAIM, and to address the new questions this first local bridge raises.

## ACKNOWLEDGMENTS

The authors thank the referees for their shrewd comments and suggestions.