How fluorine substituents strengthen aryl C–H bonds

Abstract

We have investigated the nature and bond dissociation energies (BDEs) of the aromatic C–H bonds in fluorinated benzenes C₆R₅H (each R can be H or F) using quantitative Kohn–Sham molecular orbital theory and a matching energy decomposition analysis (EDA). The C–H bond becomes stronger as the number of fluorine atoms in the benzene ring increases. This increase in the calculated BDE is additive and most pronounced for ortho-substituted C–H bonds. Our analyses of the C–H bond between and H˙ reveal that a fluorine inductive effect is responsible for this. Fluorine polarizes the closed-shell molecular orbitals of away from the carbon radical center and in this way reduces Pauli repulsion between [C˙] and the H˙ radical, leading to a stronger C–H bond. The ortho effect can be accurately modelled by a combination of Pauli repulsion (main contribution) and orbital interactions. We extend our analysis to other substituents, including ones with the opposite effect on C–H bond strength.

---

Introduction

Fluorine alters the physical properties of chemical species. For instance, it modifies the lipophilicity and acid–base properties, leading to the adoption of fluorine-containing organic compounds, including numerous F- or CF₃-substituted aromatics, in pharmaceuticals, imaging, agrochemicals, and many other areas.^1a–g In materials science, fluorinated conducting materials and fluorinated anions are important, especially in rechargeable batteries.^1i–k Through its deep influence on the molecular dipole and quadrupole, it can serve as a conformational tool in organic and biological chemistry.^1l The formation of the halogen bond is strongly associated with fluorinated organic molecules such as iodopentafluorobenzene.^1m Through the modification of intermolecular interactions, it enhances charge mobilities in organic semiconductors.¹ⁿ Zooming in on some fundamental properties of fluorinated benzenes and polyaromatics, fluorine has been found to decrease the aromaticity of aromatic rings while enhancing π-stacking properties. However, the decrease in aromaticity of the ring is accompanied by increased thermostability and resistance to chemical attack, a phenomenon known as fluoroaromaticity.^1o–r The alternating stacking of benzene·hexafluorobenzene cocrystals and related cocrystals in a parallel but displaced geometry was originally attributed to interactions of the different quadrupoles of the two components. Although electrostatics make a major contribution to the interactions, the displaced geometry is suggested to result from Pauli repulsion and dispersion forces, illustrating the effectiveness of Energy Decomposition Analysis in these systems.^1s

Another key effect of fluorine substitution is its impact on the strength of aryl C–H bonds, typically quantified by the homolytic bond dissociation energy (BDE), which can be modulated by both the number and the position of fluorine atoms.² Several groups have investigated the correlations between C–H BDEs in (poly)fluorinated benzenes and C–metal bonds in corresponding metal aryl derivatives, demonstrating that the C–H bond strength increases upon fluorination.³ For instance, the calculated C–H bond dissociation energy (BDE) for pentafluorobenzene (C₆F₅H: 123.7 kcal mol⁻¹) is about 6 kcal mol⁻¹ higher than that for benzene (C₆H₆: 117.7 kcal mol⁻¹). Perutz, Eisenstein, Jones and their co-workers have also found that the most significant increases in bond strength occur when the fluorine atom is in the ortho position relative to the C–H bond analyzed.⁴ For those reactions that are thermodynamically driven, particularly reversible reactions, this can result in ortho regioselectivity as demonstrated in numerous experimental studies.³ These findings are intriguing, since many examples in the literature have shown that neighboring bulky groups tend to weaken the C–H bond through steric repulsion,⁵ but it is unclear how this will apply to fluorine whose van der Waals radius is only slightly larger than that of hydrogen. Although the bond dissociation energies of (poly)fluorinated benzenes have been reported, a detailed explanation, firmly grounded in quantum mechanics, of why fluorine substituents strengthen C–H bonds in fluorinated benzenes is still lacking in the literature.

In this work, we unravel the physical mechanism underlying the strengthening of C–H bonds in (poly)fluorinated benzenes. We also address the origin of the variation in strengthening of the C–H bonds that follows the order ortho ≫ para > meta.⁴ To this end, we investigate the C–H bonding nature in the systems C₆R₅–H (each R can be H or F; see Fig. 1) using quantitative Kohn–Sham molecular orbital theory (KS-MO) combined with a matching energy decomposition analysis (EDA).^6,7 Our results reveal that the main factor responsible for the strengthening of the C–H bonds upon fluorination is the inductive polarization by fluorine of the closed-shell molecular orbital density away from the aromatic ring. This lowers the Pauli repulsion in the bonding region and thus strengthens the C–H bond. The EDA analysis reveals that the regioselectivity of the C–H bond energies, especially the ortho/para ratio, can be reproduced in part by the Pauli repulsion but an accurate representation requires the sum of Pauli repulsion and orbital terms. We further include C₆RH₄–H systems (R = Cl, Br, I, Li) to assess the effect of different substituents beyond fluorine.

Fig. 1 Fluorinated benzene C₆R₅H (R = H, F) systems studied in this work.

Computational methods

Computational details

All calculations were carried out using the Amsterdam Density Functional (ADF) program (ADF2019.3 for Potential Energy Surface (PES) scans and ADF2024.1 for all other computations),⁸ which is part of the Amsterdam Modeling Suite (AMS2024.1).⁹ Geometries and energies were obtained with the BLYP level of the generalized gradient approximation (GGA).¹⁰ Dispersion interactions were accounted for using the DFT-D3(BJ) method developed by Grimme and coworkers,¹¹ which includes the damping function proposed by Becke and Johnson.¹² Scalar relativistic effects were treated using the zeroth-order regular approximation (ZORA).¹³ Molecular orbitals (MOs) were expanded in a large, uncontracted set of Slater-type orbitals (STOs), specifically the TZ2P basis set, which is of triple-ζ quality and includes two sets of polarization functions.¹⁴ All electrons were treated variationally. Radical fragments were treated using a spin-unrestricted formalism. The numerical accuracy¹⁵ was set to VERYGOOD. All optimized structures were confirmed as true minima by vibrational frequency analyses, showing no imaginary frequencies.¹⁶ The bond energies reported in this work correspond to electronic energies and do not include zero-point energy (ZPE) corrections. Multiple regression analyses were conducted using the statsmodels Python library,¹⁷ while graphical visualization of the results was performed with the Matplotlib library.¹⁸ The relative bond enthalpies computed at our final level of theory, ZORA-BLYP-D3(BJ)/TZ2P, are in good agreement with experimentally validated literature values,¹⁹ supporting the reliability of our methodological approach (see Table S1).

Activation strain model and energy decomposition analysis

The overall homolytic C–H bond enthalpy ΔH, corresponding to −ΔH_BDE, between and H˙ in C₆R₅–H (R = H, F), is shown in Scheme 1. The corresponding value of ΔE is decomposed into two major components using the activation strain model (ASM, eqn (1)):²⁰

ΔE = ΔE_strain + ΔE_int (1)

here, the strain energy ΔE_strain is the energy penalty required to deform the aryl fragment from its equilibrium structure to the geometry that it acquires in the final molecule. The interaction energy ΔE_int accounts for all chemical interactions between the geometrically deformed fragments in C₆R₅–H.

Scheme 1 Formation of the C–H bond in C₆R₅H (R = H, F).

The interaction energy ΔE_int is further analyzed within the framework of the quantitative Kohn–Sham molecular orbital (KS-MO)⁶ theory by partitioning it using our canonical energy decomposition analysis (EDA) scheme into electrostatic interactions, Pauli repulsion, (attractive) orbital interactions, dispersion corrections, and spin polarization (eqn (2)):⁷

ΔE_int = ΔV_elstat + ΔE_Pauli + ΔE_oi + ΔE_disp + ΔE_spinpol (2)

The electrostatic energy ΔV_elstat corresponds to the electrostatic interactions between the unperturbed charge distribution of the distorted fragments, which is usually attractive. The Pauli repulsion ΔE_Pauli comprises the destabilizing interactions between occupied orbitals (or, more precisely, same-spin electrons on either fragment) and is responsible for any steric repulsion. The orbital interactions ΔE_oi term, accounts for stabilizing orbital interactions between the fragments including both polarization and overlap effects and can be further decomposed into the electron pair-bond energy ΔE_pb and ΔE_rel (eqn (3)). ΔE_pb is defined as the energy change associated with the formation of a doubly occupied bonding combination of the two SOMOs while all other virtual orbitals are deleted. The relaxation energy ΔE_rel, results from full relaxation after including all virtual orbitals. The ΔE_rel term includes both charge transfer interactions (donor–acceptor interaction between an occupied orbital of one fragment with an empty orbital of the other fragment) and polarization effects (empty/occupied) orbital mixing on one fragment due to the presence of another fragment.⁷

ΔE_oi = ΔE_pb + ΔE_rel (3)

The dispersion energy ΔE_disp is added as a correction.¹¹ Finally, the ΔE_spinpol term refers to the spin polarization of the spin-α and spin-β electrons of the deformed unrestricted fragments and is destabilizing (i.e., the deformed unrestricted fragments without spin polarization lie consistently 2–4 kcal mol⁻¹ higher in energy and therefore have a too stabilizing ΔE_int).²¹ The open-source PyFrag2019 program was used to automate analyzing the bonding mechanism as a function of the C–H distance.²²

Results and discussion

General trends in bond strength

The bond enthalpies ΔH (see Scheme 1) under standard conditions (298.15 K and 1 atm) and the C–H bond lengths (r_C–H) of the C₆R₅–H model systems (R = H, F), obtained from our ZORA-BLYP-D3(BJ)/TZ2P computations, are collected in Table 1. The computed C–H bond strengths as represented by the homolytic bond enthalpies ΔH, are in excellent agreement with the previously reported results by Clot, Perutz, Eisenstein, and co-workers,^3c,4 confirming that the introduction of fluorine atoms into the benzene ring strengthens the C–H bond for all substitution patterns. The calculated bond enthalpies relative to benzene (B) vary from ΔΔH = −0.4 kcal mol⁻¹ for the meta-C–H bond in fluorobenzene (m) to ΔΔH = −6.4 kcal mol⁻¹ for 1,3,5-trifluorobenzene (2o-p). Fluorination of the benzene ring also leads to a shortening of the C–H bonds. The variations are very small but are approximately proportional to the bond strengthening, with stronger bonds generally associated with shorter bond lengths (Fig. S1). The largest contraction of Δr_C–H = −0.004 Å occurs in 2o-p, i.e., the system with the strongest C–H bond.

Table 1 Bond enthalpies and energies (ΔH and ΔE; in kcal mol⁻¹), activation strain model terms (in kcal mol⁻¹), and bond lengths (in Å) of the C–H bond in C₆R₅H (R = H, F)^a

System	rC–H	ΔH	ΔΔH	ΔE	ΔΔE	ΔEstrain	ΔΔEstrain	ΔEint	ΔΔEint
a Computed at ZORA-BLYP-D3(BJ)/TZ2P, enthalpies ΔH at 298.15 K and 1 atm. Energies relative to benzene (B) are given as ΔΔ.
B	1.088	−109.3	0.0	−115.9	0.0	1.8	0.0	−117.7	0.0
o	1.086	−111.8	−2.5	−118.4	−2.5	1.8	0.0	−120.2	−2.5
m	1.087	−109.7	−0.4	−116.3	−0.4	1.8	0.0	−118.1	−0.4
p	1.087	−110.6	−1.3	−117.2	−1.3	1.7	−0.1	−118.9	−1.2
o-m	1.087	−111.8	−2.5	−118.3	−2.4	1.9	0.1	−120.2	−2.5
m-p	1.087	−110.7	−1.4	−117.3	−1.4	1.8	0.0	−119.1	−1.4
2o	1.085	−114.8	−5.5	−121.1	−5.2	1.8	0.0	−122.9	−5.2
2m	1.087	−110.2	−0.9	−116.8	−0.9	1.9	0.1	−118.7	−1.0
o-p	1.086	−113.0	−3.7	−119.5	−3.6	1.7	−0.1	−121.2	−3.5
o-m′	1.086	−112.2	−2.9	−118.7	−2.8	1.8	0.0	−120.5	−2.8
o-m-p	1.086	−112.7	−3.4	−119.2	−3.3	1.9	0.1	−121.1	−3.4
2m-p	1.086	−110.9	−1.6	−117.5	−1.6	1.9	0.1	−119.4	−1.7
2o-m	1.085	−114.7	−5.4	−121.1	−5.2	1.8	0.0	−122.9	−5.2
o-2m	1.086	−112.4	−3.1	−118.9	−3.0	1.9	0.1	−120.8	−3.1
o-m-p′	1.085	−113.1	−3.8	−119.6	−3.7	1.8	0.0	−121.4	−3.7
2o-p	1.084	−115.7	−6.4	−122.1	−6.2	1.7	−0.1	−123.8	−6.1
o-2m-p	1.086	−113.0	−3.7	−119.5	−3.6	1.9	0.1	−121.4	−3.7
2o-2m	1.086	−115.0	−5.7	−121.3	−5.4	1.9	0.1	−123.2	−5.5
2o-m-p	1.085	−115.5	−6.2	−121.8	−5.9	1.8	0.0	−123.6	−5.9
2o-2m-p	1.085	−115.5	−6.2	−121.8	−5.9	1.9	0.1	−123.7	−6.0

---

As previously reported,^4,23 and further confirmed here via multiple regression analysis, the ΔΔH values can be expressed as a linear function (r² = 0.9890) of the number of ortho- (x_ortho), meta- (x_meta), and para-fluorine substituents (x_para), as shown in eqn (4) and Fig. 2a.

ΔΔH = a + bx_ortho + cx_meta + dx_para (4)

Fig. 2 Multiple linear regression relating the number of fluorine substituents at the ortho, meta, and para positions to relative (a) C–H bond enthalpies ΔΔH, (b) electrostatic interactions ΔΔV_elstat, (c) Pauli repulsions ΔΔE_Pauli, and (d) orbital interactions ΔΔE_oi; (e) energy decomposition analysis for all EDA terms from regression coefficients for all 20 fluorinated benzenes relative to benzene; and (f) energy decomposition analysis terms in fluorobenzene relative to benzene at a consistent geometry with a C–H distance of 1.088 Å.²⁴ Computed at ZORA-BLYP-D3(BJ)/TZ2P.

The regression coefficients b, c, and d, represent the energy variation associated with the addition of a fluorine atom at the ortho, meta, and para positions, respectively. Accordingly, fluorination at the ortho position strengthens the C–H bond by approximately 2.5 ± 0.1 kcal mol⁻¹, while meta- and para-fluorine substitution increase the bond strength by only 0.2 ± 0.1 and 0.9 ± 0.1 kcal mol⁻¹, respectively. It is worth noting how closely the regression coefficients calculated here match those reported in the literature^3c using the B3PW91 functional (a = −0.27 ± 0.12 vs. −0.17 ± 0.07, b = −2.49 ± 0.07 vs. −2.49 ± 0.05, c = −0.15 ± 0.07 vs. −0.07 ± 0.05, and d = −0.81 ± 0.07 vs. −0.87 ± 0.10 kcal mol⁻¹), indicating the robustness of the results with respect to the level of theory. Since these values match so closely with two different functionals, there is no reason to expect significant differences in the EDA analysis. We therefore carried out further analysis with the ZORA-BLYP-D3(BJ)/TZ2P level of theory only.

As shown previously, the ΔΔH values cluster into three groups corresponding to a progressive increase in C–H bond strength as the number of ortho-F substituents increases (0 → 1 → 2). Because the para-coefficient (d = −0.87 kcal mol⁻¹) is significant, a secondary subdivision is observed within each region, distinguishing systems that contain a para-fluorine from those without a para-fluorine. Overall, the linear correlation demonstrates that substituting a fluorine atom for hydrogen at a specific position has an additive effect on the C–H bond strength, with this effect being most pronounced at the ortho position, followed by the para and meta positions, respectively. Here we only analyze the changes in the C–H bond dissociation energy. In contrast,²³ the C–F bond dissociation energy also changes but in the opposite direction, such that C–F BDEs decrease markedly with the number of ortho fluorine substituents, with smaller effects from meta- and para-fluorine substitution.

The overall trends in bond enthalpies ΔH are well reproduced by the electronic bond energies ΔE, as shown by the relative ΔΔ values in Table 1 and the multiple regression analyses shown in Fig. S2. Therefore, to elucidate the origin of the C–H bond strengthening in fluorinated benzenes, we analyzed the electronic bond energy ΔE using the Activation Strain Model (ASM).²⁰ Within this framework, ΔE is decomposed into a strain energy ΔE_strain and an interaction energy ΔE_int (eqn (1); see Computational methods). The corresponding ASM terms are also listed in Table 1. The only significant geometrical deformation observed upon bond formation is the in-plane bending of the neighboring substituents (H or F) in the aryl fragment. This deformation leads to a small strain energy ΔE_strain, which remains essentially constant across different substitution patterns. Thus, the trends in bond strength (i.e., in both ΔH and ΔE) are entirely governed by the interaction energy ΔE_int. We further investigated the underlying C–H bonding mechanism using Kohn–Sham Molecular Orbital theory (KS-MO)⁶ and a matching Energy Decomposition Analysis (EDA).⁷ This method decomposes ΔE_int into distinct physically meaningful components: electrostatic interactions ΔV_elstat, Pauli repulsion ΔE_Pauli, orbital interactions ΔE_oi, among others (see eqn (2) in the Computational methods section). The EDA terms for the C–H bonds in fluorinated benzenes, relative to benzene, can also be expressed as linear functions of the number of ortho- (x_ortho), meta- (x_meta), and para-fluorine substituents (x_para), as depicted in Fig. 2b–d. The corresponding absolute values are provided in Table S1. The graphs in Fig. 2a–c exhibit the characteristic patterns corresponding to the dominance of the ortho contribution (x_ortho = 0, 1, or 2) as shown by the rings around the points. Within each subset, we have all possible combinations of meta- and para-Fs (6 species for 0 and 2 ortho-Fs, 7 species for 1 ortho-F). The linear correlations indicate that each para substitution has a similar effect (see above). The contribution from meta substituents is negligible in Fig. 2a and c (ΔΔH and ΔΔE_Pauli) but is comparable to the para contribution for ΔΔV_elstat. The spans of the electrostatic, Pauli, and orbital contributions are 8.5, 13.6, and 3.6 kcal mol⁻¹, respectively (Fig. 2b–d).

The trends in the interaction energy, and therefore in the C–H bond enthalpies, are primarily dictated by the Pauli repulsion ΔE_Pauli, as shown in Fig. 2c and e and through linear correlation between ΔΔH and ΔΔE_Pauli (Fig. 3a). Relative to benzene, ΔE_Pauli is the most stabilizing contribution, becoming less repulsive as the number of ortho-F substituents increases and exhibiting the same clustering pattern observed for the bond enthalpies (Fig. 2c). Fluorination at the ortho position decreases the Pauli repulsion by approximately 4.8 ± 0.2 kcal mol⁻¹, while para-F substitution reduces ΔE_Pauli by 2.5 ± 0.3 kcal mol⁻¹. The stronger effect of ortho compared to para substitution is quantified by the ratio b/d. The value of b/d for ΔΔH is 2.86 ± 0.33. The corresponding value for ΔΔE_Pauli is 1.91 ± 0.26 thus accounting for approximately 2/3 of the ortho preference.²⁵ Similar to the bond enthalpy, fluorination at the meta position has almost no effect on the Pauli repulsion.

Fig. 3 Linear correlation between ΔΔH and (a) ΔΔE_Pauli, (b) ΔΔE_Pauli + ΔΔE_oi; (c) multiple linear regression relating the number of fluorine substituents at the ortho, meta, and para positions to ΔΔE_Pauli + ΔΔE_oi. Computed at ZORA-BLYP-D3(BJ)/TZ2P.

In contrast to ΔE_Pauli, the electrostatic interaction term ΔV_elstat follows the opposite trend to ΔH, becoming less attractive as the number of ortho-F substituents increases, as indicated by the positive regression coefficients (Fig. 2b). The magnitude of the electrostatic effects follows the pattern ortho ≫ meta > para, which may be ascribed to simple distance effects, 1/r, arising because of the 9 electrons of fluorine compared to 1 of hydrogen. Finally, unlike the other EDA terms, the orbital interaction term ΔE_oi does not exhibit the same clustering pattern based on the number of ortho-F substituents and therefore follows a different trend from the bond strength (Fig. 2d). In fact, the stabilization upon fluorine substitution is larger for the meta position than for the ortho position (−0.9 ± 0.1 versus −0.5 ± 0.1 kcal mol⁻¹, respectively). In contrast, para-F substitution leads to a destabilization of 0.8 ± 0.2 kcal mol⁻¹ in the orbital interaction term. This destabilizing contribution at the para position accounts for the smaller, remaining difference between the ortho- and para-C–H bond strengths, as evidenced by the improved correlation coefficient and the closer match between the b/d ratios obtained from the multiple regressions of ΔΔE_oi + ΔΔE_Pauli and of ΔΔH (2.99 ± 0.28 and 2.86 ± 0.33, respectively, Fig. 3b and c).

A systematic analysis of the ortho-, meta- and para-fluorine substituent effect on the EDA terms, using different systems as initial references, reveals trends similar to those obtained from the regression coefficients when all 20 C–H bonds are considered (Fig. 2e, f and S5).²⁶ Additionally, the same conclusions obtained from the analysis at equilibrium and consistent geometries can be drawn from the analysis of the EDA terms as a function of the C–H bond distance (see Fig. S7 and S8).

Our results discussed so far demonstrate that the Pauli repulsion provides a unified explanation for both the overall C–H bond strengthening in (poly)fluorobenzenes and the major part of the selective increase in bond strength for the C–H bond ortho to fluorine, which is fully captured by the combined Pauli and orbital interaction terms. In this way, we address two important open questions of which the underlying causes had remained unknown until now.^2b A comprehensive analysis of the origin of the Pauli repulsion trends is provided in the following section.

Origin of the Pauli repulsion

As discussed above, our quantitative MO and EDA analyses reveal that Pauli repulsion ΔE_Pauli plays a key role in the strengthening of the C–H bonds upon fluorine substitution on the benzene ring. As will become clear in the following, the inductive effect of fluorine polarizes the closed-shell molecular orbitals of the aryl fragment towards the fluorine atom, which reduces their amplitude in the C–H bonding region of the radical fragment. This, in turn, leads to a decrease in Pauli repulsion with the singly-occupied hydrogen 1s orbital, resulting in a stronger C–H bond.

In the EDA framework, the ΔE_Pauli term arises from repulsive interactions between same-spin occupied orbitals on the respective fragments.⁷ Therefore, all ΔE_Pauli originates from occupied–occupied orbital interactions between the singly occupied hydrogen 1s orbital and the same-spin component of the occupied molecular orbitals in the σ-electron system of the aryl fragment (see Fig. 4a). The greater the orbital overlap between the hydrogen 1s and the aryl σ-MOs, the stronger the Pauli repulsion. Table S3 lists key overlaps between the hydrogen 1s AO and the σ-MOs of the aryl fragment. The occupied orbital most affected by the presence of a fluorine substituent is the σ_HOMO−6 of the aryl radical (Fig. 4b). The orbital overlap between the hydrogen 1s orbital and the aryl σ_HOMO−6 orbital amounts to 〈H_1s|σ_HOMO−6〉 = 0.18 for benzene (Fig. 4c). This orbital overlap significantly decreases for the ortho- and para-C–H bonds, namely to 〈H_1s|σ_HOMO−6〉 = 0.12 and 0.09, respectively. This leads to the lower Pauli repulsion for the latter two bonds. Note that for the meta-C–H bond, the overlap remains nearly unchanged: 〈H_1s|σ_HOMO−6〉 = 0.19; hence the ΔE_Pauli value remains similar to that in benzene.

Fig. 4 (a) Schematic MO diagram of the C–H bond formation in benzene, (b) σ_HOMO−6 occupied FMO ψ_HOMO−6 (isovalue = 0.03 au) of the aryl fragment, and (c) overlapping FMOs, 〈H_1s|σ_HOMO−6〉 orbital overlap (above arrows), and overlap density given by the product of fragment orbitals ψ_1s·ψ_HOMO−6 (isovalue = 0.002 au). Computed at ZORA-BLYP-D3(BJ)/TZ2P.

A similar, though less pronounced, effect is observed in the 〈H_1s|σ_HOMO−2〉 orbital overlap. Fig. S9 presents the 〈H_1s|σ_HOMO−6〉 and 〈H_1s|σ_HOMO−2〉 overlaps as a function of the C–H bond distance. Interestingly, although the largest reduction in Pauli repulsion occurs for the ortho-C–H bond, the largest decrease in 〈H_1s|σ_HOMO−6〉 overlap occurs at the para position. However, when considering the combined effect of both the 〈H_1s|σ_HOMO−6〉 and 〈H_1s|σ_HOMO−2〉, a more significant reduction is observed for the ortho-C–H bond, whereas the para position shows only a modest reduction in the second overlap. An approximate way to assess the overall contribution of all orbitals to the Pauli repulsion is by considering the sum of the squared overlaps (S²) between the H_1s orbital and all same-spin occupied σ-orbitals of the aryl fragment.²⁷ Accordingly, the ∑S² values for the C–H bond in benzene and for the ortho-, meta-, and para-C–H bonds in fluorobenzene are 0.244, 0.227, 0.242, and 0.234, respectively. These values follow the same trend observed for both the Pauli repulsion and the C–H bond strength.

The essential question remaining is: why does fluorination reduce the Pauli-repulsive overlap between the occupied orbitals on the fragment and the radical electron of the hydrogen atom and thus strengthen the aryl C–H bond? To address this question, we analyze the formation of the aryl fragment in benzene (R = H) and fluorobenzene (R = F) from the biradical and R˙ (see Fig. 5a). This approach enables us to understand how the relevant orbitals arise from the same biradical, and how their nature becomes different in the case of a fluorine substituent R = F in fluorobenzene versus a hydrogen substituent R = H in benzene. The resulting MO diagram is shown qualitatively in Fig. 5b (for details, see Fig. S10–S12 and Tables S4–S6). From this point onward, we focus on the ortho position, where the reduction of Pauli repulsion is most pronounced and where the proximity of the substituent would intuitively suggest a repulsive effect.

Fig. 5 (a) Fragmentation of the C–R bond (R = F versus H). (b) Qualitative MO diagram showing the formation of the σ_SOMO and σ_HOMO−6 of the aryl radical from the interaction of SFO gross Mulliken contributions (in gray) are shown only for σ_SOMO and σ_HOMO−6; for details on the other orbitals, see Fig. S8 and Table S4. Computed at ZORA-BLYP-D3(BJ)/TZ2P.

The σ_HOMO−6 orbitals of the and radicals arise from a σ-bonding interaction between a doubly occupied σ-orbital of the biradical and the singly occupied hydrogen 1s (H_1s) or fluorine 2p (F_2p) orbitals, respectively (see Fig. 5b). The F_2p orbital lies lower in energy than the H_1s orbital because fluorine is more electronegative than hydrogen. As a result, the F_2p orbital is closer in energy to the doubly occupied σ-orbital of the fragment than the H_1s orbital, which leads to a larger fluorine contribution to the formation of the bonding combination, which is the σ_HOMO−6 orbital (see gross Mulliken contributions in Fig. 5b). This greater contribution of fluorine polarizes the σ_HOMO−6 orbital towards the fluorine atom and away from the carbon radical. Consequently, the overlap between σ_HOMO−6 and the hydrogen 1s orbital across the aromatic C–H bond between the carbon-radical center of and H˙ is reduced. This is a clear manifestation of fluorine's inductive effect, which predominantly affects the σ-system in fluorobenzene, as previously discussed elsewhere.²⁸ Furthermore, note that neither the F_2p nor the H_1s contributes significantly to the formation of the σ_SOMO in the aryl radical. As a result, the σ_SOMO is not polarized by the substituent, leading to similar 〈H_1s|σ_SOMO〉 orbital overlaps and, consequently, similar electron-pair bond ΔE_pb and orbital interaction ΔE_oi energies for benzene and fluorobenzene (see Fig. S9 for overlaps and Table S7 for ΔE_pb values).²⁹

Predictive power of the model

Finally, we wish to extend our model to related systems. To this end, we have extended our MO and EDA analyses of the C–H bonds to a broader set of monosubstituted benzenes C₆H₅R with R still including H and F, but now also Cl, Br, I, and Li. Proceeding from the model we have established, one would expect that electronegative substituents R strengthen C–H bonds, unless, in cases of adjacent C–H bonds, R becomes so big that it destabilizes this C–H bond through direct R⋯H Pauli repulsion.⁵ On the other hand, electropositive substituents R are then expected to weaken the C–H bonds because they polarize orbital amplitude toward the pertinent carbon atoms. The results confirm this expectation.

The bond dissociation enthalpies ΔH, bond dissociation energies ΔE, together with the ASM and EDA terms for this series of monosubstituted benzenes, are given in Table S7. The C–H bond dissociation enthalpy is larger than that in benzene for C₆FH₅ (−111.8 kcal mol⁻¹) and C₆ClH₅ (−110.9 kcal mol⁻¹) but very close to that in benzene (−109.3 kcal mol⁻¹) for the other halogens (Br, −110.3; I, −109.4 kcal mol⁻¹). However, the C–H bond dissociation enthalpy is significantly smaller in C₆LiH₅ (−91.9 kcal mol⁻¹). A comparison with experimental BDE data would be valuable; however, such data are not listed in ref. 19 for the systems with R = Cl, Br, I, and Li. As expected, the ΔH values follow the same trends as the bond energies ΔE, which are in turn dominated by the interaction energy ΔE_int. The relative EDA terms for the C–H bonds in the C₆RH₅ (R = F, Cl, Br, I, Li) systems are depicted in Fig. 6a. As shown in Fig. 6b, the inductive effect responsible for the reduction in Pauli repulsion in fluorobenzene is also observed in the other halobenzenes. In contrast to the case of fluorine, the orbital term is of the same order of magnitude as the Pauli repulsion term for Cl, Br, and I. Thus, the increase in ΔE_oi compensates for the variation in the ΔE_Pauli and results in values of ΔE_int for Cl, Br, and I within 1.4 kcal mol⁻¹ of one another (Fig. 6a and Table S7).

Fig. 6 (a) Energy decomposition analysis of the ortho-C–H bond in mono-substituted benzenes C₆RH₅ (R = F, Cl, Br, I, and Li) relative to benzene at a consistent geometry with a C–H distance of 1.088 Å, and (b) the associated repulsive overlap 〈H_1s|σ_HOMO−6〉. Computed at ZORA-BLYP-D3(BJ)/TZ2P.

As in the case of fluorine, strong contributions to the decrease of the Pauli repulsion are found in σ_HOMO−6. The halogen substituent polarizes the σ_HOMO−6 orbital toward itself and away from the C–H bonding region, leading to a progressive decrease in 〈H_1s|σ_HOMO−6〉 overlap as the electronegativity of the substituent increases from I to F. However, the EDA results in Fig. 6a (the complete dataset is provided in Table S7) reveal that inductive effects alone cannot fully account for the observed trends in Pauli repulsion across the halobenzene series. Although bromine and iodine are more electronegative than hydrogen and thus exert a higher inductive effect, the Pauli repulsion is actually higher in bromobenzene and iodobenzene compared to benzene. This apparent discrepancy arises from the increase in size of the substituent atom. As the size of R increases, the spatial extension of the valence AOs and the number of subvalence shells increase, causing a stronger Pauli repulsion. Thus, the observed trends in Pauli repulsion reflect a balance between two opposing effects, namely, the electron-withdrawing (inductive) effect, which reduces Pauli repulsion, and the steric (size-related) effect, which enhances Pauli repulsion.

For fluorine and chlorine, the most electronegative atoms in the series, the atom size effect is offset by the electron-withdrawing (inductive) effect, resulting in reduced 〈H_1s|σ_HOMO−6〉 overlap (see Fig. 6b), negative ΔΔE_Pauli, and stronger C–H bonds compared to benzene. In the case of bromine, the atom size effect begins to outweigh the inductive effect, leading to a modest increase in Pauli repulsion (ΔΔE_Pauli = +0.9 kcal mol⁻¹). Nevertheless, this destabilization is compensated by an enhanced orbital interaction,³⁰ still resulting in a stronger ortho-C–H bond compared to benzene. For iodobenzene, the atom size effect becomes dominant, resulting in greater Pauli repulsion. However, even in this case, the gain in attractive orbital interactions counterbalances the increase in ΔE_Pauli, yielding a C–H bond with nearly the same strength as in benzene.

Finally, we have also investigated phenyllithium (C₆H₅Li), a system in which the substituent R is less electronegative than hydrogen. The 2s orbital of lithium lies significantly higher in energy than the np orbitals of the halogens and does not contribute meaningfully to the formation of the σ_HOMO−6 orbital in the aryl radical. As a result, the 〈H_1s|σ_HOMO−6〉 orbital overlap in phenyllithium is nearly identical to that in benzene (see Fig. 7a). However, lithium does affect higher-energy orbitals of the aryl fragment. Notably, Fig. 7b shows a substantial increase in the 〈H_1s|σ_HOMO−1〉 orbital overlap in C₆H₄Li–H relative to both benzene and fluorobenzene. Unlike the halogens, lithium raises the orbital amplitude in the C–H bonding region rather than depleting it, effectively donating electron density into the carbon scaffold of C–C σ bonds through an effect opposite to that of fluorine. This electron-donating character leads to increased occupied–occupied overlap and thus increased Pauli repulsion across the aromatic C–H bond. This is reminiscent of the lone-pair shielded radical effect induced by the substituent on the C–C bond strength.³¹ This effect outweighs all stabilizing contributions, resulting in a substantially weaker C–H bond.

Fig. 7 For the model systems with R = H, F, and Li: (a) H_1s AO and C₆H₄R˙ σ_HOMO−6 (isovalue = 0.03 au) at their position in C₆H₄R–H, 〈H_1s|σ_HOMO−6〉 orbital overlap, and overlap density given by the product of fragment orbitals ψ_1s·ψ_HOMO−6 (isovalue = 0.002 au). (b) H_1s AO and C₆H₄R˙ σ_HOMO−1 (isovalue = 0.03 au) at their position in C₆H₄R–H, 〈H_1s|σ_HOMO−1〉 orbital overlap, and overlap density given by the product of fragment orbitals ψ_1s·ψ_HOMO−1 (isovalue = 0.002 au). Computed at ZORA-BLYP-D3(BJ)/TZ2P.

To disentangle size and inductive effects, we also performed EDA analyses for the same substituents at the para position. In this configuration, direct R⋯H repulsion is eliminated, allowing the intrinsic electronic nature of the substituent (electron-withdrawing or electron-donating) to dominate the Pauli repulsion trends. As shown in Fig. S13, all substituents more electron-withdrawing than hydrogen lead to a decrease in Pauli repulsion, following the electronegativity trend. In contrast, for phenyllithium, the only electron-donating substituent, Pauli repulsion remains significantly higher than in benzene, even in the absence of direct steric interactions.

Conclusion

Fluorine substitution strengthens the C–H bond in fluorinated benzenes C₆R₅H (R = H, F), as indicated by the increase in the homolytic bond dissociation energies (BDE). Our quantum chemical bonding analyses, based on dispersion-corrected, (scalar) relativistic density functional theory, uncover that fluorination strengthens the C–H bond in large part by reducing steric (Pauli) repulsion between filled orbitals. As demonstrated previously, the effects of multiple fluorine substitution are additive and markedly greater for ortho than for para substitution while the effect of meta substitution is barely significant. The ortho effect on BDE is calculated to be ca. 2.9 times the effect of para substitution. The EDA analysis shows that the Pauli repulsion is an important contributor to the bond strengthening but that this ortho/para ratio can be reproduced by the sum of the Pauli repulsion and orbital terms.

Our quantitative MO and energy decomposition analyses reveal how the strong inductive effect of the fluorine atom withdraws electron density, especially from the ipso carbon of the aryl fragment and polarizes its closed-shell molecular orbitals towards the fluorine atom. This polarization reduces the spatial extension of the occupied orbitals on the aryl fragment towards the hydrogen 1s orbital. Therefore, the same-spin orbital overlap across the C–H bond shrinks, and Pauli repulsion becomes weaker, resulting in stronger C–H bonds.

Extension of the analysis to other monosubstituted benzenes C₆RH₅ (R = Cl, Br, I, Li) shows that substituents R modulate the strength of the adjacent C–H bond by reducing or enhancing Pauli repulsion, either via an indirect, inductive effect or via direct steric repulsion between C–R and C–H bonds. For R = F, the C–H BDE is notably bigger than for other halogens which, in turn, are slightly bigger than that of benzene. For R = F, the strong inductive effects decrease the overall Pauli repulsion, resulting in a stronger C–H bond. For R = Cl, the overall Pauli repulsion remains smaller than for benzene. In contrast, for R = Br and I, weaker inductive effects are offset by larger atomic size, which increases direct R⋯H steric repulsion. For all halogens except fluorine, this repulsion is balanced by enhanced orbital attraction, yielding C–H bonds of similar strength for R = Cl, Br, I, and H. Thus, fluorine remains as a special case among electron-withdrawing substituents that we have studied so far. Finally, for R = Li, the inductive effect is inverted: lithium pushes electron density to the aryl fragment, increasing Pauli repulsion and leading to the weakest C–H bond in the series.

Author contributions

RNP, OE, and FMB conceived the project, which was supervised by DRS, AACB, and CFG. DASO carried out the quantum-chemical computations and bonding analyses and drafted the manuscript. All authors discussed the results and reviewed the manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

All data of this study are available in the main text and supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6sc01846g.

Acknowledgements

The authors thank the Dutch Research Council (NWO) for financial support. This work was carried out on the ADA (VU) Cluster and the Dutch national e-infrastructure with the support of SURF Cooperative. Additionally, D. A. S. O. acknowledges the São Paulo Research Foundation (FAPESP) (Grants #2022/01685-6 and #2024/03646-3) and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (Grant #162366/2021-3) for funding their doctoral scholarship. A. A. C. B. thanks the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Brazil) for financial support (Grant #313720/2023-1). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES), Finance Code 001, Program 33002010191P0. OE acknowledges the support of the Research Council of Norway through the Centre of Excellence Hylleraas Centre for Quantum Molecular Sciences (No. 262695).

References

1. (a) D. O'Hagan, Understanding organofluorine chemistry. An Introduction to the C−F Bond, Chem. Soc. Rev., 2008, 37, 308–319; (b) M. Inoue, Y. Sumii and N. Shibata, Contribution of Organofluorine Compounds to Pharmaceuticals”, ACS Omega, 2020, 5, 10633–10640; (c) Y. Zhou, J. Wang, Z. Gu, S. Wang, W. Zhu, J. L. Aceña, V. A. Soloshonok, K. Izawa and H. Liu, Next Generation of Fluorine-Containing Pharmaceuticals: Compounds Currently in Phase II–III Clinical Trials of Major Pharmaceutical Companies: New Structural Trends and Therapeutic Areas, Chem. Rev., 2016, 116, 422–518; (d) M. Bernús, G. D. Núñez, W. C. Hartley, M. Guasch, J. Mestre, M. Besora, J. J. Carbó and O. Boutureira, Impact of Fluorine Pattern on Lipophilicity and Acid–Base Properties of 2-(Thiofluoroalkyl)pyridines: Insights from Experiments and Statistical Modeling, J. Med. Chem., 2025, 68, 4787–4800; (e) Y. Du, Y. Bian, D. Baecker, G. Dhawan, A. Semghouli, L. Kiss, W. Zhang, A. E. Sorochinsky, V. A. Soloshonok and J. Han, Fluorine in the Pharmaceutical Industry: FDA-Approved Fluorine-Containing Drugs in 2024, Chem.–Eur. J., 2025, 31, e202500662; (f) E. Henary, S. Casa, T. L. Dost, J. C. Sloop and M. Henary, The Role of Small Molecules Containing Fluorine Atoms in Medicine and Imaging Applications, Pharmaceuticals, 2024, 17, 281; (g) S. Preshlock, M. Tredwell and V. Gouverneur, ¹⁸F-Labeling of Arenes and Heteroarenes for Applications in Positron Emission Tomography, Chem. Rev., 2016, 116, 719–766; (h) P. Jeschke, Recent Developments in Fluorine-Containing Pesticides, Pest Manage. Sci., 2024, 80, 3065–3087; (i) T. Zhang, Z. Chen, W. Zhang, L. Wang and G. Yu, Recent Progress of Fluorinated Conjugated Polymers, Adv. Mater., 2024, 36, 2403961; (j) H. Hyrondelle, A. Terry, J. Lhoste, S. Tencé, K. Lemoine, J. Olchowka, D. Dambournet, C. Tassel, J. Gamon and A. Demourgues, Fluorine as a Key Element in Solid-State Chemistry of Mixed Anions 3d Transition Metal-Based Materials for Electronic Properties and Energy, Chem. Rev., 2025, 125, 4287–4358; (k) Y. Wang, X. Yang, Y. Meng, Z. Wen, R. Han, X. Hu, B. Sun, F. Kang, B. Li, C. Wang and G. Wang, Fluorine Chemistry in Rechargeable Batteries: Challenges, Progress, and Perspectives, Chem. Rev., 2024, 124, 3494–3589; (l) L. Hunter, The C–F Bond as a Conformational Tool in Organic and Biological Chemistry, Beilstein J. Org. Chem., 2010, 6, 38; (m) G. Cavallo, P. Metrangolo, R. Milani, T. Pilati, A. Priimagi, G. Resnati and G. Terraneo, The halogen bond, Chem. Rev., 2016, 116, 2478–2601; (n) B. Maiti, A. Schubert, S. Sarkar, S. Bhandari, K. Wang, Z. Li, E. Geva, R. J. Twieg and B. D. Dunietz, Enhancing Charge Mobilities in Organic Semiconductors by Selective Fluorination: A Design Approach Based on a Quantum Mechanical Perspective, Chem. Sci., 2017, 8, 6947–6953; (o) N. M. P. Rosa, M. Máximo-Canadas and I. Borges Jr., Assessing the Aromaticity of Fluorinated Benzene Derivatives Using New Descriptors Based on the Distributed Multipole Analysis (DMA) Partition of the Electron Density, J. Phys. Org. Chem., 2025, 38, e70024; (p) J. J. Torres-Vega, A. Vásquez-Espinal, L. Ruiz, M. A. Fernández-Herrera, L. Alvarez-Thon, G. Merino and W. Tiznado, ChemistryOpen, 2015, 4, 302–307; (q) S. Balser and O. V. Boltalina, Electron Deficient Fluorous Polycyclic Aromatic and Heteroaromatic Hydrocarbons, Coord. Chem. Rev., 2026, 554, 217592; (r) T. J. Fuhrer, M. Houck and S. T. Iacono, Fluoromaticity: The Molecular Orbital Contributions of Fluorine Substituents to the π-Systems of Aromatic Rings, ACS Omega, 2021, 6, 32607–32617; (s) K. Carter-Fenk and J. M. Herbert, Electrostatics does not dictate the slip-stacked arrangement of aromatic π–π interactions, Chem. Sci., 2020, 11, 6758–6765.
2. (a) O. Eisenstein, J. Milani and R. N. Perutz, Selectivity of C−H activation and competition between C−H and C−F bond activation at fluorocarbons, Chem. Rev., 2017, 117, 8710–8753; (b) Y. P. Budiman, R. N. Perutz, P. G. Steel, U. Radius and T. B. Marder, Applications of Transition Metal-Catalyzed ortho-Fluorine-Directed C−H Functionalization of (Poly)fluoroarenes in Organic Synthesis, Chem. Rev., 2024, 124, 4822–4862; (c) T. Rogge, N. Kaplaneris, N. Chatani, J. Kim, S. Chang, B. Punji, L. L. Schafer, D. G. Musaev, J. Wencel-Delord, C. A. Roberts, R. Sarpong, Z. E. Wilson, M. A. Brimble, M. J. Johansson and L. Ackermann, C–H activation, Nat. Rev. Methods Primers, 2021, 1, 43; (d) T. Dalton, T. Faber and F. Glorius, C−H Activation: Toward Sustainability and Applications, ACS Cent. Sci., 2021, 7, 245–261; (e) T. P. Bast, J. V. Obligacion, E. Rochette, I. Pappas and P. J. Chirik, Cobalt-Catalyzed Borylation of Fluorinated Arenes: Thermodynamic Control of C(sp2)−H Oxidative Addition Results in ortho-to-Fluorine Selectivity, J. Am. Chem. Soc., 2019, 141, 15378–15389.
3. (a) W. D. Jones and E. T. Hessel, Photolysis of Tp’Rh(CN-neopentyl)(η²-PhN=C=N-neopentyl) in Alkanes and Arenes: Kinetic and Thermodynamic Selectivity of [Tp’Rh(CN-neopentyl)] for Various Types of C−H Bonds, J. Am. Chem. Soc., 1993, 115, 554–562; (b) B. A. Vanchura, S. M. Preshlock, P. C. Roosen, V. A. Kallepalli, R. J. Staples, R. E. Maleczka Jr., D. A. Singleton and M. R. Smith, Electronic Effects in Iridium C−H Borylations: Insights from Unencumbered Substrates and Variation of Boryl Ligand Substituents, Chem. Commun., 2010, 46, 7724–7726; (c) E. Clot, C. Mégret, O. Eisenstein and R. N. Perutz, Validation of the M−C/H−C Bond Enthalpy Relationship through Application of Density Functional Theory, J. Am. Chem. Soc., 2006, 128, 8350–8357; (d) P. J. Chirik and T. P. Pabst, A Tutorial on Selectivity Determination in C(sp2)−H Oxidative Addition of Arenes by Transition Metal Complexes, Organometallics, 2021, 40, 813–831.
4. (a) E. Clot, M. Besora, F. Maseras, C. Mégret, O. Eisenstein, B. Oelckers and R. N. Perutz, Bond energy M–C/H–C correlations: dual theoretical and experimental approach to the sensitivity of M–C bond strength to substituents, Chem. Commun., 2003, 490–491; (b) E. Clot, C. Mégret, O. Eisenstein and R. N. Perutz, Exceptional Sensitivity of Metal-Aryl Bond Energies to ortho-Fluorine Substituents: Influence of the Metal, the Coordination Sphere, and the Spectator Ligands on M−C/H−C Bond Energy Correlations, J. Am. Chem. Soc., 2009, 131, 7817–7827; (c) M. E. Evans, C. L. Burke, S. Yaibuathes, E. Clot, O. Eisenstein and W. D. Jones, Energetics of C−H Bond Activation of Fluorinated Aromatic Hydrocarbons Using a [Tp’Rh(CNneopentyl)] Complex, J. Am. Chem. Soc., 2009, 131, 13464–13473; (d) T. Tanabe, W. W. Brennessel, E. Clot, O. Eisenstein and W. D. Jones, Synthesis, structure, and reductive elimination in the series Tp’Rh(PR₃)(Ar^F)H; Determination of rhodium–carbon bond energies of fluoroaryl substituents, Dalton Trans., 2010, 39, 10495–10509; (e) E. Clot, B. Oelckers, A. H. Klahn, O. Eisenstein and R. N. Perutz, cis–trans Isomerisation of CpRe(CO)₂(H)(Ar^F) (Ar^F = C₆F_nH_5-n; n = 0–5) is the rate determining step in C–H activation of fluoroarenes: a DFT study, Dalton Trans., 2003, 4065–4072.
5. (a) E. Blokker, W.-J. van Zeist, X. Sun, J. Poater, J. M. van der Schuur, T. A. Hamlin and F. M. Bickelhaupt, Methyl Substitution Destabilizes Alkyl Radicals, Angew. Chem. Int. Ed., 2022, 61, e202207477 (Angew. Chem., 2022, 134, e202207477); (b) P. Vermeeren, W.-J. van Zeist, T. A. Hamlin, C. Fonseca Guerra and F. M. Bickelhaupt, Not Carbon s–p Hybridization, but Coordination Number Determines C–H and C–C Bond Length, Chem.–Eur. J., 2021, 27, 7074–7079; (c) J. Poater, R. Visser, M. Solà and F. M. Bickelhaupt, Polycyclic Benzenoids: Why Kinked is More Stable than Straight, J. Org. Chem., 2007, 72, 1134–1142; (d) J. Poater, J. Paauwe, S. Pand, G. Merino, C. Fonseca Guerra and F. M. Bickelhaupt, Kekulene: Structure, stability and nature of H•••H interactions in large PAHs, Mol. Astrophys., 2017, 8, 19–26; (e) D. Rodrigues Silva, P. Vermeeren, J. M. van der Schuur and F. M. Bickelhaupt, Radical Stability Paradox: Substituent Effects versus Heats of Formation, Chem.–Eur. J., 2025, e03600.
6. (a) R. van Meer, O. V. Gritsenko and E. J. Baerends, Physical Meaning of Virtual Kohn−Sham Orbitals and Orbital Energies: An Ideal Basis for the Description of Molecular Excitations, J. Chem. Theory Comput., 2014, 10, 4432–4441; (b) T. A. Albright, J. K. Burdett and W. H. Whangbo, Orbital Interactions in Chemistry, Wiley, New York, 2nd edn, 2013, ISBN 978-0-471-08039-8.
7. F. M. Bickelhaupt and E. J. Baerends, in Reviews in Computational Chemistry, ed. K. B. Lipkowitz and D. B. Boyd, Wiley–VCH, Hoboken, 2000, pp. 1–86.
8. (a) ADF2019.3/2024.1, SCM: Theoretical Chemistry, Vrije Universiteit Amsterdam, The Netherlands, 2019, http://www.scm.com; (b) ADF2024.1, SCM, Theoretical Chemistry, Vrije Universiteit Amsterdam, The Netherlands, 2024, http://www.scm.com; (c) G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders and T. Ziegler, Chemistry with ADF, J. Comput. Chem., 2001, 22, 931–967; (d) C. Fonseca Guerra, J. G. Snijders, G. te Velde and E. J. Baerends, Towards an order-N DFT method, Theor. Chem. Acc., 1998, 99, 391–403.
9. (a) E. J. Baerends, N. F. Aguirre, N. D. Austin, J. Autschbach, F. M. Bickelhaupt, R. Bulo, C. Cappelli, A. C. T. van Duin, F. Egidi, C. Fonseca Guerra, A. Förster, M. Franchini, T. P. M. Goumans, T. Heine, M. Hellström, C. R. Jacob, L. Jensen, M. Krykunov, E. van Lenthe, A. Michalak, M. M. Mitoraj, J. Neugebauer, V. P. Nicu, P. Philipsen, H. Ramanantoanina, R. Rüger, G. Schreckenbach, M. Stener, M. Swart, J. M. Thijssen, T. Trnka, L. Visscher, A. Yakovlev and S. van Gisbergen, The Amsterdam Modeling Suite, J. Chem. Phys., 2025, 162, 162501; (b) AMS 2024.1, SCM: Theoretical Chemistry, Vrije Universiteit Amsterdam, The Netherlands, 2024, http://www.scm.com.
10. (a) A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A, 1988, 38, 3098–3100; (b) C. T. Lee, W. T. Yang and R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B, 1988, 37, 785–789.
11. (a) S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132, 154104; (b) S. Grimme, S. Ehrlich and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, J. Comput. Chem., 2011, 32, 1456–1465.
12. E. R. Johnson and A. D. Becke, A post-Hartree–Fock model of intermolecular interactions, J. Chem. Phys., 2005, 123, 024101.
13. E. van Lenthe, E. J. Baerends and J. G. Snijders, Relativistic total energy using regular approximations, J. Chem. Phys., 1994, 101, 9783–9792.
14. E. van Lenthe and E. J. Baerends, Optimized Slater-type basis sets for the elements 1–118, J. Comput. Chem., 2003, 24, 1142–1156.
15. (a) M. Franchini, P. H. T. Philipsen, E. van Lenthe and L. Visscher, Accurate Coulomb Potentials for Periodic and Molecular Systems through Density Fitting, J. Chem. Theory Comput., 2014, 10, 1994–2004; (b) M. Franchini, P. H. T. Philipsen and L. Visscher, The Becke Fuzzy Cells Integration Scheme in the Amsterdam Density Functional Program Suite, J. Comput. Chem., 2013, 34, 1819–1827.
16. (a) A. Bérces, R. M. Dickson, L. Fan, H. Jacobsen, D. Swerhone and T. Ziegler, An implementation of the coupled perturbed Kohn-Sham equations: perturbation due to nuclear displacements, Comput. Phys. Commun., 1997, 100, 247–262; (b) H. Jacobsen, A. Bérces, D. P. Swerhone and T. Ziegler, Analytic second derivatives of molecular energies: a density functional implementation, Comput. Phys. Commun., 1997, 100, 263–276; (c) S. K. Wolff, Analytical second derivatives in the Amsterdam density functional package, Int. J. Quantum Chem., 2005, 104, 645–659.
17. S. Seabold and J. Perktold, Statsmodels: Econometric and statistical modeling with Python, Proceedings of the 9th Python in Science Conference, 2010.
18. J. D. Hunter, Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 2007, 3, 90–95.
19. Y. R. Luo, Comprehensive handbook of chemical bond energies. CRC press, Boca Raton, 2007.
20. (a) P. Vermeeren, T. A. Hamlin and F. M. Bickelhaupt, Chemical reactivity from an activation strain perspective, Chem. Commun., 2021, 57, 5880–5896; (b) P. Vermeeren, S. C. C. van der Lubbe, C. Fonseca Guerra, F. M. Bickelhaupt and T. A. Hamlin, Understanding chemical reactivity using the activation strain model, Nat. Protoc., 2020, 15, 649–667; (c) F. M. Bickelhaupt and K. N. Houk, Analyzing Reaction Rates with the Distortion/Interaction-Activation Strain Model, Angew. Chem. Int. Ed., 2017, 56, 10070–10086 (Angew. Chem., 2017, 129, 10204–10221); (d) L. P. Wolters and F. M. Bickelhaupt, The activation strain model and molecular orbital theory, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2015, 5, 324–343; (e) I. Fernández and F. M. Bickelhaupt, The activation strain model and molecular orbital theory: understanding and designing chemical reactions, Chem. Soc. Rev., 2014, 43, 4953–4967; (f) W.-J. van Zeist and F. M. Bickelhaupt, The activation strain model of chemical reactivity, Org. Biomol. Chem., 2010, 8, 3118–3127; (g) F. M. Bickelhaupt, Understanding reactivity with Kohn–Sham molecular orbital theory: E2–S_N2 mechanistic spectrum and other concepts, J. Comput. Chem., 1999, 20, 114–128.
21. (a) F. M. Bickelhaupt, M. Solà and C. Fonseca Guerra, Highly polar bonds and the meaning of covalency and ionicity—structure and bonding of alkali metal hydride oligomers, Faraday Discuss., 2007, 135, 451–468; (b) D. Rodrigues Silva, E. Blokker, J. M. van der Schuur, T. A. Hamlin and F. M. Bickelhaupt, Nature and strength of group-14 A–A’ bonds, Chem. Sci., 2024, 15, 1648–1656; (c) E. Blokker, X. Sun, J. Poater, J. M. van der Schuur, T. A. Hamlin and F. M. Bickelhaupt, The Chemical Bond: When Atom Size Instead of Electronegativity Difference Determines Trend in Bond Strength, Chem.–Eur. J., 2021, 27, 15616–15622.
22. (a) W. J. van Zeist, C. Fonseca Guerra and F. M. Bickelhaupt, PyFrag—Streamlining Your Reaction Path Analysis, J. Comput. Chem., 2008, 29, 312–315; (b) X. Sun, T. M. Soini, J. Poater, T. A. Hamlin and F. M. Bickelhaupt, PyFrag 2019—Automating the Exploration and Analysis of Reaction Mechanisms, J. Comput. Chem., 2019, 40, 2227–2233; (c) X. Sun, T. Soini, L. P. Wolters, W.-J. van Zeist, C. Fonseca Guerra, T. A. Hamlin and F. M. Bickelhaupt, PyFrag 2019, Vrije Universiteit Amsterdam, The Netherlands, 2019.
23. Concomitant changes are also observed in the C–F bond dissociation energies. S. A. Macgregor, D. McKay, J. A. Panetier and M. K. Whittlesey, Computational study of the hydrodefluorination of fluoroarenes at [Ru(NHC)(PR₃)₂(CO)(H)₂]: predicted scope and regioselectivities, Dalton Trans., 2013, 42, 7386–7396.
24. To ensure a consistent comparison across systems, all analyses for the monosubstituted system were carried out at a fixed C–H bond length of 1.088 Å, which corresponds to the equilibrium bond length in benzene. In this way, the original variations in bonding capabilities along the various model systems are separated from the resulting geometrical relaxations. Table S2 and Fig. S4 report the EDA terms at this consistent geometry for a representative subset of systems (Fig. S3). Comparison of the data in Tables S1 and S4 shows that the same conclusions are obtained for both equilibrium and fixed geometries.
25. The b/d ratio can be interpreted as a measure of how much stronger the ortho-C–H bond becomes for each 1 kcal mol⁻¹ increase in the para-C–H bond strength. As such, it provides a meaningful quantitative descriptor for analyzing the regioselectivity of C–H bond strength between the ortho and para positions.
26. An approximately additive effect for the polyfluorinated benzene is also observed for both the bond enthalpies and Pauli repulsion when the changes in their values are compared with those of the monosubstituted system (Fig. S6).
27. The use of the sum of squared overlaps (∑S²) provides an approximate, although representative, measure of the overall Pauli repulsive interactions. A more rigorous descriptor would require weighting each S² contribution by the corresponding kinetic energy of the interacting orbitals, which more accurately reflects their Pauli repulsive character. However, such an approach would substantially increase the complexity of the analysis and obscure the clear trends discussed here. For theoretical details see reference 7.
28. C. Nieuwland, P. Vermeeren, F. M. Bickelhaupt and C. Fonseca Guerra, Understanding chemistry with the symmetry-decomposed Voronoi deformation density charge analysis, J. Comput. Chem., 2023, 44, 2108–2119.
29. As the electron pair-bond energy is approximately proportional to Δε + c • S²/Δε (c is a constant) it is not possible to deduce the variation of this term from the variation in Δε (see also Ref. 6b).
30. The increase in atomic size contributes not only to greater Pauli repulsion but also to enhanced orbital interaction stabilization of the C–H bond. In the case of halogens, this stabilization arises primarily from an increase in the C–H electron-pair bond (ΔE_pb), which becomes progressively more significant from fluorine to iodine (see values in Table S7). For R = Li, however, the observed increase in orbital attraction cannot be attributed solely to ΔE_pb. Instead, other orbital interaction components, such as, donor–acceptor interactions and polarization, also play a substantial role in the overall stabilization.
31. E. Blokker, M. ten Brink, J. M. van der Schuur, T. A. Hamlin and F. M. Bickelhaupt, Origin of the Captodative Effect: The Lone-Pair Shielded Radical, ChemistryEurope, 2023, 1, e202300006.

---