Matta, Boyd. The quantum theory of atoms in molecules

432 16 Topological Properties of the Electron Distribution in Hydrogen-bonded Systems

related to the local energy contribution of the Laplacian of rðrÞ by the local form of the virial theorem (Eq. 1 in a.u.):

	1	‘2rðrÞ ¼ 2GðrÞ þ VðrÞ	ð1Þ
4
Because GðrÞ > 0 and VðrÞ < 0, local depletion of the			charge density

ð‘2rðrÞ > 0Þ is related to a preponderance of twice the kinetic energy density over the potential energy density, whereas the local concentration ð‘2rðrÞ < 0Þ corresponds to the opposite situation. These local energy magnitudes cannot be derived from the electron density distribution alone; knowledge of the wavefunction is also necessary. Thus, the exact magnitudes of GðrÞ and VðrÞ (and therefore also that of HðrÞ) cannot be extracted from experimental analysis of rðrÞ in crystals. They can, however, be estimated by use of Eq. (1) and Abramov’s functional (Eq. 2 in a.u.), which links the topological properties of the electron distribution (rðrÞ, ‘rðrÞ and ‘2rðrÞ) to the local kinetic energy density [11, 24, 25]:

ð

Þ ¼

10 ð

Þ

ð

Þ

þ

72 j rðrÞ

þ

6

ð

Þ

ð Þ

G r

3

3p2

2=3r r

5=3

1

‘rðrÞj

1

‘2r r

2

At BCP positions the second term of Eq. (2) vanishes, because the first derivative becomes zero, and the local energy densities are functionals of rb and ‘2rb. It should be noted that in the original work of Abramov, good agreement with Hartree–Fock calculations of GðrÞ is obtained in the medium-range region only, i.e. for distances of approximately 0.5–2.1 A˚ from atomic nuclei.

Equations (1) and (2) were applied for the first time to a large experimental data set of 83 XaH O (X ¼ C, N, O) hydrogen-bonding interactions for which BCPs were experimentally observed between 0.5 and 1.2 A˚ from the hydrogen atom and between 0.5 and 1.6 A˚ from the oxygen atom [24]. For this data set the estimated values of Gb and Vb were exponentially dependent on the H O distance (Fig. 16.6), and there was an exponential interdependence between them. Later, similar features were also found for H F hydrogen-bonded systems by use of theoretical data [15].

Several theoretical articles have dealt with the validity of these local energy estimates and their associated limitations. Analysis along the bond paths of the HB complexes formed by (FH)2 and (OH2)2 in their minimum energy configuration show that the Abramov functional (Eq. 2) elegantly reproduces the magnitude of the electron kinetic energy density well in the intermediate H A region where the BCP is located but fails in regions close to the nuclei, especially around the hydrogen nucleus [26]. For a set of 37 H F HBs, comparison of the Gb magnitudes calculated by use of Eq. (2) with those obtained from the ab-initio wave function reveals an almost perfect match [27]. In the same study, estimated values of the potential component, Vb, and the total energy density Hb, derived

a0 3)

16.3 Energy Properties at the Bond Critical Point 433

Fig. 16.6 Exponential relationships between Gb and Vb (kJ mol 1 and the H O distance (A˚ ). (Values taken from Ref. [24]). The fitted curves are Gb ¼ 12ð2Þ 103 exp½ 2:73ð9ÞdðO HÞ& and

Vb ¼ 50:0ð1:1Þ 103 exp½ 3:6dðO HÞ&.

from application of the virial theorem were, again, both in very good agreement with the corresponding theoretical magnitudes. Study of NaH covalent and N H hydrogen-bonding interactions shows that the Abramov functional works well for HBs longer than 1.9 A˚ , overestimating the ab-initio results for closer interaction distances and providing random results for covalent NaH bonds [28]. For dihydrogen bonded complexes it has been found that Eq. (2) underestimates Gb magnitudes at long distances and overestimates them at short distances by up to 10% [21].

In addition to the exponential dependencies on H O distance, excellent linear relationships between V and the sum of the perpendicular curvatures, l1 þ l2, and between G and the curvature along the bond path direction, l3, have been found for both experimental H O and the theoretical H F data sets. It is particularly worthy of note that the linear regressions performed with these data sets have equivalent fitting data (Fig. 16.7). The topological and energy dependencies on H A distance observed for the reported experimental data have enabled interpretation of the strengthening of the HB interaction in terms of those rðrÞ properties [14, 24]. Local energy densities are dimensionally equivalent to a pressure. Hence VðrÞ is interpreted as the pressure exerted on the rðrÞ distribution to concentrate it and GðrÞ is interpreted as the counterpart pressure exerted by rðrÞ against the environment, as a reaction to the former because of the electron– electron repulsion. Thus, the greater the potential energy density the stronger the repulsion, leading to a greater kinetic energy.

434 16 Topological Properties of the Electron Distribution in Hydrogen-bonded Systems

Fig. 16.7 Linear relationships between local energy densities Gb and Vb (kJ mol 1 a0 3) and electron density curvatures (e A˚ 5) for experimental O H and calculated F H hydrogen bonds. (Values taken from Refs [14] and [15]). The linear fitting corresponds to the equations Vb ¼ 35:1ð7Þðl1 þ l2Þ and Gb ¼ 15:3ð1Þl3 for the (O H) HB.

As a consequence of bond formation, a maximum distribution of rðrÞ is created along the bond path; this is the indication of the established interaction. The curvatures of the electron distribution in the perpendicular plane at any point of the bond path, and in particular at the BCP, increase with increasing concentration of rðrÞ. Along the bond path rðrÞ decreases from the nuclei towards the BCP, where it has a local minimum, and the curvature along this direction increases with charge depletion. The increase of the negative magnitude l1 þ l2 is accompanied by an increase of the positive amplitude, l3, indicating that a sharper electron concentration in the plane where rðrÞ is maximum is necessarily followed by larger rðrÞ depletion along the direction of the bond path. Thus, at BCP, Vb and Gb are related to the charge concentration of rðrÞ in the perpendicular plane to the bond path and to its charge depletion along the path direction, respectively. The linear relationships observed between the local energy densities and the topological curvatures reflect this situation. In this way, when the H A distance shortens and the HB interaction becomes stronger, accumulation of charge in the internuclear region increases in magnitude with all topological and energetic properties at the BCP. Indeed, for a pure closed shell interaction, and as a consequence of Pauli’s principle, the increase of rb is accompanied by a more important rðrÞ depletion ð‘2rb > 0Þ that follows from l3 > jl1 þ l2j within the range of distances considered. According to the linear correlations Gb z l3 and Vb z l1 þ l2, this is related to the observed ratio of pressures jVbj=Gb < 1 for this type of interaction.

16.4 Topological Properties and Interaction Energy 435

16.4

Topological Properties and Interaction Energy

The interaction energy, which is the stabilization occurring as a result of complex formation, provides a measure of the strength of the interaction. It is calculated as the di erence between the energy of the complex and the energies of the isolated monomers, and corresponds to the negative value of the dissociation energy.

With the objective of gaining insight into the connection between the local behavior of the electron distribution in hydrogen-bonding regions and the integrated properties of complexes, the dependence of the interaction energy on the BCP properties has been explored in parallel with the distance dependencies of the latter. This approach has been undertaken to characterize the strength of hydrogen bonds in molecular crystals, thus providing a link between local rðrÞ quantities associated with intermolecular interactions and crystal properties.

With the initially proposed linear dependencies of rb and ‘2rb on hydrogen bonding distance, linear correlations of these topological properties with the interaction energy have also been reported for XCN HCl and XCN HF complexes (Fig. 16.8), and for carbenes as HB acceptors [8–10]. Similar results have been also recently been obtained for a wide variety of HB systems with a wide range of interaction strengths, from weakly bonded complexes in the van der Waals limit, for example CH4 Ar, to the almost covalent interaction in the H3Oþ H2O complex [29]. In this last study, linear regression was applied to the rb data represented against the interaction energy in the full interaction range

Fig. 16.8 Dependence of rb (a.u.) on the interaction energy (kJ mol 1) [9]. The linear regressions are: Ei ¼ 13ð1Þ þ 2133ð66Þrb and

Ei ¼ 21:8ð6Þ þ 2558ð30Þrb for the XCN HCl and XCN HF complexes, respectively.

436 16 Topological Properties of the Electron Distribution in Hydrogen-bonded Systems

taken into account. As for ‘2rb, however, although the proportionality is valid for weak HBs the regression should be regarded as an unpolished approximation for strong interactions. According to these analyses, the accumulation of electrons between the acceptor and the hydrogen, represented by rb, reflects the interaction strength, characterization that has also been used for other type of interaction, for example SiaO bonds in silicates [30]. In this last example the interaction energy calculated from rb confirms predictions of bond strength by empirical models.

On the basis of experimental analysis of electron density for 83 XaH O hydrogen-bonding interactions in crystals, a di erent approach to the characterization of the interaction energy has been proposed in terms of the local energy densities at BCP [31]. By comparing the dependence of Vb on the HB distance observed for this set of experimental H O interactions with that for the dissociation energies, De, theoretically calculated for several dimers, an approximately linear relationship between Vb and the interaction energy, Ei ¼ De A1=2Vb was found.

By following this phenomenological relationship, derivation of an interaction potential for H O hydrogen bonds was undertaken [31] on the basis of the topological analysis of rðrÞ for this experimentally characterized interaction. The interaction potential function is defined as:

U ¼ nHb

ð3Þ

where n is a constant in volume units (0.982 a0 3) which was calculated by using the force constant of H O hydrogen bonds in ice VIII (k ¼ 22:7 N m 1) experimentally determined from Raman spectroscopy [32]. Here Hb is expressed as the sum of the exponential dependencies fitted for Gb and Vb (Fig. 16.6). From Eq. (3) it should be noted that the potential well depth at equilibrium geometry corresponds to ðHbÞmax, which represents the greatest excess of kinetic energy that r can a ord at BCP, leading to the most e cient r depletion at the interatomic surface. Thus, the internuclear distance associated with ðHbÞmax (i.e. equilibrium geometry) is the best compromise between the quantity of bonding charge and its degree of depletion (reflected by rb and ‘2rb, respectively) furnishing the deepest stabilization of a pure closed-shell interaction [15]. According to the terms involved in Eq. (3), the interaction potential is expressed as:

U ¼ 49100 expð 3:6rÞ 11800 expð 2:73rÞ

ð4Þ

where U is in kJ mol 1 and r is the H O interatomic distance in A˚ . It is worthy of note that when the interaction energy function Ei ¼ De ¼ 1=2Vb has been corrected for the polarization energy jU0 polj at their equilibrium distance, which was theoretically calculated for ice VIII [33], the new function De þ jU0 polj crosses the U ¼ nHb potential at its minimum. As far as DeðrÞ represents the interaction energy at any equilibrium distance r, the crossing betweenDeðrÞ þ jU0 polj and U ¼ nHb at Umin elegantly indicates the internal agreement between both descriptions (Fig. 16.9).

16.4 Topological Properties and Interaction Energy 437

Fig. 16.9 Interaction potential UðrÞ ¼ nHb and polarization-corrected dissociation energy DeðrÞ þ jU0 polj ¼ 1=2Vb (kJ mol 1) along the interatomic distance (A˚ ), with the Morse ðUMorseÞ and Buckingham ðUBuckgÞ potentials for the hydrogen bond [31].

Both U ¼ nHb and De ¼ 1=2Vb were checked against experimentally derived thermodynamic properties of ice VIII, leading to a good estimate of the heat of sublimation at 0 K and excellent agreement with the H O expansion distance from 0 to 273 K calculated from the linear coe cient of expansion of ice. Several energy properties, for example polarization, cohesion, and binding energies, calculated by use of those functions were also found to be in good agreement with theoretical calculations for ice VIII and ice Ih for experimental and calculated geometries. The U ¼ nHb potential was also successfully compared with Buckingham and Morse-type potentials, which are used in semi-empirical atom–atom potential methods and in spectroscopy, respectively, for representing hydrogen bond energy. Indeed, when they were constrained to have the same features (i.e. the same position, potential, and curvature) as those describing the potential curve U ¼ nHb at its minimum, the comparisons revealed almost perfect matching over the complete range of distances considered.

In a later theoretical study involving H F hydrogen-bonding interactions [15] the proportionality factor obtained (@0.42) was similar to that derived between EI and Vb for H O interactions (@0.5). In this work, study of the (FH)2 system at di erent intermolecular distances, also showed that the distance of the minimum interaction energy at equilibrium geometry also corresponds to the maximum of Hb [15]. Both theoretical findings support the U ¼ nHb and De A 1=2Vb correspondences that were derived from experimental data. In the same theoretical study, the degree of softening term (SD), defined as Hb=rb for Hb > 0, was also used to estimate Ei magnitudes. For the data set for neutral complexes XaH FaY, plots of the theoretically calculated energies Ei against the corre-

438 16 Topological Properties of the Electron Distribution in Hydrogen-bonded Systems

sponding values of Vb and of SD gave a narrower distribution of data for the latter when both were fitted by use of linear regression.

16.5

Electron Localization Function, h(r)

The electron localization function (ELF), hðrÞ, was originally defined by Becke and Edgecombe as [34]:

hðrÞ ¼ 1=ð1 þ q2Þ

ð5Þ

where q ¼ ðT TwÞ=TTF, and T is the kinetic energy, and Tw and TTF are the von Weizsa¨cker and the Thomas–Fermi kinetic energy functionals, respectively. TTF gives the kinetic energy of a homogeneous electron gas having the same density as the point in the electron density which is under consideration and Tw is an inhomogeneity correction. ELF takes values in the range 0 < hðrÞ < 1. (See Chapter 5 of this book for a review on the ELF).

In the regions where the electrons are localized Pauli repulsion has little e ect and, therefore, T ATw and hðrÞ A1. If, however, the Pauli repulsion is strong and electrons are delocalized, ðT TwÞ is large and hðrÞ A0. The value hðrÞ ¼ 0:5 corresponds to the electron localization of a homogenous electron gas. ELF local maxima are known as attractors and correspond to local charge concentrations. For these attractors, the topological partition of the ‘hðrÞ gradient vector field yields basins that can be associated either with core electrons (core basins) or with bonds and lone pairs (valence basins). The valence basin is characterized by the synaptic order, which is the number of participating atomic valence shells. The synaptic order is given by the number of core basins sharing a boundary with the valence basin, plus the number of nuclei of hydrogen atoms it contains. Depending on their synaptic order, valence basins are classified as monosynaptic, disynaptic, trisynaptic, etc. Monosynaptic basins correspond to the lone pairs of the Lewis model, and polysynaptic basins to the shared pairs of the Lewis model. Disynaptic basins therefore correspond to two-center bonds, trisynaptic basins to three-center bonds, etc. The basins are delimited by zero-flux surfaces S (‘hðrÞ nðrÞ ¼ 0, Er A S), where ð3; 1Þ critical points topologically analogous to BCPs are also found.

For hydrogen bonds two valence basins appear in the bonding region. One, which corresponds to the electron pair that belongs to the DaH bond, is a disynaptic basin in contact with the core basin of the donor and contains the hydrogen. The other is a monosynaptic basin that corresponds to the acceptor lone pair. Three ð3; 1Þ critical points are observed in DaH A hydrogen bonds, two between the core and the valence basins, the other in the surface separating both valence basins. The behavior of ELF in the hydrogen-bonding region is illustrated by the ELF profiles depicted in Fig. 16.10.

16.5 Electron Localization Function, hðrÞ	439

Fig. 16.10 hðrÞ along the N HF hydrogen bond of the N2 HF and H3N HF complexes. The N F distance has been normalized to the same value in both traces [35].

Topological analysis of the ELF has been used to characterize the strength of the HB interactions on the basis of the values of hðrÞ at the ð3; 1Þ critical points [36]. This is done by use of the core-valence bifurcation index, which is defined as:

Q ¼ hvv hcv

ð6Þ

where hvv is the value of hðrÞ at the ð3; 1Þ point between both valence basins and hcv is the largest of the hðrÞ values at the two ð3; 1Þ points between the core and valence basins.

The Q index has been interpreted from the localization domains, which are the volumes enclosed by isosurfaces of the ELF (hðrÞ ¼ f , where f is a constant with 0 < f < 1). Hence, if Q is negative, in the range hvv < f < hcv, there are separated localization domains for the donor and the acceptor. If Q is positive, however, a localization domain containing both valence attractors is observed. Q < 0 corresponds to weak complexes, for example N2 HF, whereas Q > 0 is observed for stronger complexes, for example H3N HF (both are depicted in Fig. 16.10). In Fig. 16.10, the di erent sign of Q can be deduced from the values of ELF at the ð3; 1Þ critical points, which are the minima in the plot. For weak complexes the value of hðrÞ at the critical point between both valence basins ðhvvÞ, which is the central minimum in the plot, is below the value of hðrÞ for the other two minima; the opposite situation is observed for the strong complexes. It has been shown that the more negative the value of Q, the stronger the interaction; there is an almost linear correlation between this index and the dissociation energy of the complexes [35] for a given donor (Fig. 16.11).

440 16 Topological Properties of the Electron Distribution in Hydrogen-bonded Systems

Fig. 16.11 Relationship between dissociation energy (kJ mol 1) and core-valence bifurcation index (Q) for FaH A complexes [36]. The fitted line corresponds to the equation De ¼ 334ð22ÞQ þ 26ð1Þ.

The changes both in the distribution of the ELF magnitudes and in the observed domains during the proton-transfer process in strong HB complexes suggest that hydrogen bonding is accompanied by strong electron localization in the intermolecular regions. This has been shown for proton transfer between an imidazole and a carboxyl group [37]. At the middle stage of the N H O transfer, a single localization domain with an ELF value larger than 0.5 bridges both molecules through the hydrogen, indicating strong electron localization throughout the hydrogen-bonding region. This situation contrast with the initial and the final stages of proton transfer, in which the donor and acceptor valence shells seem to be separated by a region of depleted electron distribution, represented by ELF values of @0.2. Three ELF attractors are, moreover, observed in the middle stage of the proton transfer, because the hydrogen is no longer included in the valence basins of either the nitrogen or the oxygen but appears as a sharp ELF peak at the hydrogen position.

16.6

Complete Interaction Range

16.6.1

Dependence of Topological and Energy Properties on the Interaction Distance

From weak van der Waals to strong covalent interactions, hydrogen atoms interact and bind with other atoms in very di erent forms. This can be easily observed

16.6 Complete Interaction Range 441

by making use of the topological properties of rðrÞ. For example, for H X hydrogen-bonding interactions ‘2rb is positive whereas for XaH covalent bonds ‘2rb is negative. Thus, according to topological characterization of rðrÞ, ‘2rb < 0 and ‘2rb > 0 are, respectively, indicative of closed-shell and shared-shell interactions. Despite these di erences, the existence of very short HB complexes has enabled almost continuous consideration of the binding properties of the pair of atoms X and H from covalent to very weak interactions.

The first attempt to find a unique dependence of rb on the interaction distance, including covalent and HB complexes, corresponds to the CaH and C H interactions. For these a semi-logarithmic relationship was proposed for 33 systems with internuclear distances between 1.0 and 2.8 A˚ (correlation coe cient r2 ¼ 0:996) [38]:

dðCaH; C HÞ ¼ 0:52ð2Þ 0:43ð3Þ lnðrbÞ

ð7Þ

This equation can be also expressed as an exponential dependence of rb on the distance. Double logarithm functions, which are equivalent to power functions rb ¼ AdB have also been successfully used to fit similar dependencies of rb in covalent bonds, however [12]. Statistical analysis of the best fitting of rb data for 16 types of bond, including covalent and hydrogen bonds, and using either a double logarithm or single logarithm model, has been reported [39]. According to the results, the second model, which corresponds to exponential dependence, provides slightly better statistical results for all kind of interaction.

Although a single exponential enables good fitting for the complete range of distances, statistical analysis conducted later for a large set of NaH N complexes shows that a significantly better fit results when two separate exponential regression functions of the same type are used, one for the covalent NaH and the other for the H N hydrogen-bonding interactions [40]. In the same way, in a theoretical study of F H hydrogen bonded complexes [15], in which covalent FaH bonds were also included, it was proposed that data corresponding to the shared-shell and to the closed-shell regions are better fitted with independent exponentials. Then, to describe the whole interaction range with a single continuous dependence, a joint function (Eq. 8) linking both exponentials was derived. By following this method the joint function is obtained by dividing the complete range of distances into three parts, which correspond to the covalent, transition, and noncovalent interactions. After fitting of the data with single exponentials within the covalent and noncovalent regimes independently, the data included in the transition regime are fitted to a three-variable function defined as:

JðxÞ ¼	f ðxÞ	1	þ gðxÞ	1	ð8Þ
		1 þ eEðx aÞ=c		1 þ e Eðx bÞ=c

where a, b, and c are the fitting terms, E is a normalization factor, and f ðxÞ and gðxÞ are the fitted exponentials within the covalent and noncovalent ranges, respectively. Further studies of covalent bonds and hydrogen-bonding interactions