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.

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 models5,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.

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 

(1)
(2)

where s(r)=Asρ(r)/ρ(r)4/3, t(r)=Atρ(r)/ρ(r)7/6, As1=2(3π2)1/3, At=(π/3)1/6/4, μx=10/81, μc0.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 εu be an infinitesimal vector, with u=1 defined by the azimuthal θ and zenith ϕ angles (following Zwillinger’s convention21). In the vicinity of a BCP (H͇ denoting the density Hessian and λi its eigenvalues),

(3)

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

(4)

In the case δP is of ε second order, we make the angular and the distance dependences disappear by performing the following spherical average,22 

(5)

If we now specify P to be the previously mentioned DFT-related quantities, one gets (Cx=(3/4)(3/π)1/3, Cc=1, []d/dρ),

(6)
(7)
(8)

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 μx,c) is multiplied by ◇. Defining

(9)

one obtains

(10)
(11)

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,

(12)

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 ρ(rc) and 2ρ(rc) 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 β (Ti, V) and δ (Ti, Cr) agostic bonds, to which we will add a nickel(II) complex featuring an agostic Hδ that is nearly anagostic (or preagostic) (CδHδNi124°) 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 (ρ(rc),2ρ(rc)) 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.

FIG. 1.

The considered agostic and hydrogen-bonded systems.

FIG. 1.

The considered agostic and hydrogen-bonded systems.

Close modal
Table I.

Values of some density-related quantities at various BCPs (in atomic units).

 Ti(IV)V(III)Cr(II)Ni(II)H1H2aH3H4
TiHβagoTiHδagoVHβagoCrHδagoNiHδagoOHOH1OH2OHClH
ρ(rc) 0.040 0.030 0.047 0.030 0.014 0.053 0.042 0.050 0.028 0.013 
2ρ(rc) 0.141 0.114 0.202 0.120 0.046 0.147 0.141 0.168 0.083 0.036 
ρ(rc) 0.042 0.029 0.072 0.028 0.004 0.136 0.092 0.146 0.028 0.004 
Qxc02(rc) −3.297 −2.500 −3.672 −2.604 −1.946 −1.747 −1.684 −1.686 −1.698 −1.519 
Pxc02(rc) 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)H1H2aH3H4
TiHβagoTiHδagoVHβagoCrHδagoNiHδagoOHOH1OH2OHClH
ρ(rc) 0.040 0.030 0.047 0.030 0.014 0.053 0.042 0.050 0.028 0.013 
2ρ(rc) 0.141 0.114 0.202 0.120 0.046 0.147 0.141 0.168 0.083 0.036 
ρ(rc) 0.042 0.029 0.072 0.028 0.004 0.136 0.092 0.146 0.028 0.004 
Qxc02(rc) −3.297 −2.500 −3.672 −2.604 −1.946 −1.747 −1.684 −1.686 −1.698 −1.519 
Pxc02(rc) 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).

FIG. 2.

Hydrogen and agostic bonds represented by their ρ(rc) and 2(rc) values.

FIG. 2.

Hydrogen and agostic bonds represented by their ρ(rc) and 2(rc) values.

Close modal

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 OH1 bond, λ1=0.074, λ2=0.070, λ3=0.286, while for the titanium agostic one: λ1=0.047, λ2=0.011, λ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 (Qxc02(rc),Pxc02(rc)) representation is considered, two distinct and separate groups clearly appearing. This is notably due to the fact that the Qxc02(rc) 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.

FIG. 3.

Hydrogen and agostic bonds classified with respect to their Qxc02(rc) and Pxc02(rc) values.

FIG. 3.

Hydrogen and agostic bonds classified with respect to their Qxc02(rc) and Pxc02(rc) values.

Close modal

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 CO bonds in H1 and H2, and the TiCα bond, Pxc02 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 (rc) and ρ(rc) values for the studied bonds in the H3 and the chromium complexes; but, as Qxc02(rc) 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 yi for a BCP are not generally collinear to the eigenvectors for another BCP, the direct comparison of the λi has no direct meaning. So, in order not to lose information, the knowledge of the λ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: σ1=iλi, σ2=i<jλiλj, and σ3=iλi. Nevertheless, one can equivalently privilege the (σ12,) pair over (σ1,σ2), these two similarity invariants being related by σ2=(σ12)/2.

Such a choice is motivated by the fact that ◇ has a clear physical interpretation since δ2ρ(ε,u)rc=ε2(iλi2(u.yi)rc2+o(1)), so that: δ2ρrc(2)=13ρ(rc). Thus, ◇ measures the variation rate of ρ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 

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 (ρ(rc),2(rc)) 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.

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

2.
R. G.
Parr
and
W.
Yang
,
Density Functional Theory of Atoms and Molecules
(
Oxford University Press, NewYork
,
1989
).
3.
R. F. W.
Bader
,
Atoms in Molecules: A Quantum Theory
(
Oxford University Press
,
Oxford, U.K.
,
1990
).
4.
P. L. A.
Popelier
,
Atoms in Molecules An Introduction
(
Pearson Education
,
Harlow
,
2000
).
5.
6.
L.
Delle Site
,
Europhys. Lett.
57
,
20
(
2002
).
7.
H. J.
Bohórquez
and
R. J.
Boyd
,
J. Chem. Phys.
129
,
024110
(
2008
).
8.
F. A.
Bulat
,
E.
Chamarro
,
P.
Fuentealba
, and
A.
Toro-Labbé
,
J. Phys. Chem. A
108
,
342
(
2004
).
9.
S.
Liu
and
N.
Govin
,
J. Phys. Chem. A
112
,
6690
(
2008
).
10.
R. F. W.
Bader
and
H.
Essén
,
J. Chem. Phys.
80
,
1943
(
1984
);
notice that we will use the common term of “bond” even if its use has been criticized [see for instance
R. F. W.
Bader
,
J. Phys. Chem. A
113
,
10391
(
2009
)].
[PubMed]
11.
Atomic units will be used unless explicitly otherwise stated.
12.
P. S.
Svendsen
and
U.
von Barth
,
Phys. Rev. B
54
,
17402
(
1996
).
13.
S. -K.
Ma
and
K. A.
Brueckner
,
Phys. Rev.
165
,
18
(
1968
).
14.
Note that the exchange conventional gauge is used. Besides, it is also assumed that the inhomogeneity of the Laplacian near a BCP will only have a small impact on the energy densities variations. As, from the best of our knowledge, the role of the higher density derivatives (n4ρ) are unknown at least for correlation, it is anyway not possible to go behind this order; moreover, Eq. (2) is rigorously valid only in the high density limit, see
C. D.
Hu
and
D. C.
Langreth
,
Phys. Rev. B
33
,
943
(
1986
);
T.
Thonhauser
,
V. R.
Cooper
,
S.
Li
,
A.
Puzder
,
P.
Hyldgaard
, and
D. C.
Langreth
,
Phys. Rev. B
76
,
125112
(
2007
);
however, we will use it for every density regime, as done for instance by
Y.
Zhao
and
D. G.
Truhlar
,
J. Chem. Phys.
128
,
184109
(
2008
).
[PubMed]
15.
exUEG is the well known Slater’s functional, Slater,
The Self-Consistent Field for Molecules and Solids
(
McGraw-Hill
,
New York
,
1974
);
ecUEG is not exactly analytically known but the parameterization of
J. P.
Perdew
and
Y.
Wang
,
Phys. Rev. B
45
,
13244
(
1992
), can be considered a reference; as the BCP density values will typically belong to [0.007;0.300], the (spin-unpolarized) PW92 expression will be reasonably substituted by its second order diagonal Padé approximant calculated at ρ=0.06.
16.
A. D.
Becke
and
K. E.
Edgecombe
,
J. Chem. Phys.
92
,
5397
(
1990
).
17.
B.
Silvi
and
A.
Savin
,
Nature (London)
371
,
683
(
1994
).
18.
P. W.
Ayers
and
S.
Jenkins
,
J. Chem. Phys.
130
,
154104
(
2009
).
19.
R. F. W.
Bader
,
J. Chem. Phys.
73
,
2871
(
1980
).
20.
A.
Holas
and
N. H.
March
,
Int. J. Quantum Chem.
56
,
371
(
1995
).
21.
D.
Zwillinger
,
“Spherical Coordinates in Space” §4.9.3 in CRC Standard Mathematical Tables and Formulae
(
CRC
,
Boca Raton
,
1995
).
22.
V.
Tognetti
,
L.
Joubert
,
P.
Cortona
, and
C.
Adamo
,
J. Phys. Chem. A
113
,
12322
(
2009
).
23.
E.
Clot
and
O.
Eisenstein
,
Struct. Bonding (Berlin)
113
,
1
(
2004
).
24.
U.
Koch
and
P.
Popelier
,
J. Phys. Chem.
99
,
9747
(
1995
).
25.
P. L. A.
Popelier
and
G.
Logothetis
,
J. Organomet. Chem.
555
,
101
(
1998
).
26.
V.
Tognetti
,
L.
Joubert
,
R.
Raucoules
,
T.
De Bruin
, and
C.
Adamo
(unpublished).
27.
W.
Yao
,
O.
Eisenstein
, and
R. H.
Crabtree
,
Inorg. Chim. Acta
254
,
105
(
1997
).
28.
J. C.
Lewis
,
J.
Wu
,
R. G.
Bergman
, and
J. A.
Ellman
,
Organometallics
24
,
5737
(
2005
).
29.
The molecules were fully optimized using Gaussian 03 (Ref. 31) at the PBE0 level (Ref. 32), the nonmetallic atoms being described by the standard 6-31++G(d,p) Pople basis set and the metals by the Wachters+f one (Ref. 33). The topological analysis has been carried out with the MORPHY98 software (Refs. 34 and 35).
30.
Notice that the information carried by σ3 is lost, but is not relevant for our purpose. The ellipticity (which only focuses on the plane that is orthogonal to the bond path) would constitute an alternative. For the H bonds in Ref. 24, it lies in the [0.005–1.394] range, whereas the one for the agostic bonds in Ref. 25 is [0.703,1.534], so that overlap exists; thus ellipticity cannot be discriminative.
31.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
 et al., GAUSSIAN03, Revision C.02 (Gaussian, Inc., Wallingford, CT,
2004
).
32.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
33.
A. J. H.
Wachters
,
J. Chem. Phys.
52
,
1033
(
1970
).
34.
MORPHY98, a topological analysis program written by P. L. A. Popelier with a contribution from R. G. A. Bone (UMIST, Engl., EU).
35.
P. L. A.
Popelier
,
Comput. Phys. Commun.
93
,
212
(
1996
).