The Distance between Minima of Electron Density and Electrostatic Potential as a Measure of Halogen Bond Strength.

In this study, we present results of a detailed topological analysis of electron density (ED) of 145 halogen-bonded complexes formed by various fluorine-, chlorine-, bromine-, and iodine-containing compounds with trimethylphosphine oxide, Me3PO. To characterize the halogen bond (XB) strength, we used the complexation enthalpy, the interatomic distance between oxygen and halogen, as well as the typical set of electron density properties at the bond critical points calculated at B3LYP/jorge-ATZP level of theory. We show for the first time that it is possible to predict the XB strength based on the distance between the minima of ED and molecular electrostatic potential (ESP) along the XB path. The gap between ED and ESP minima exponentially depends on local electronic kinetic energy density at the bond critical point and tends to be a common limiting value for the strongest halogen bond.

1. Introduction

Among various types of σ-hole interactions, halogen bonds (XBs) are among the most known and widely investigated [1,2,3,4]. XBs were shown to play a significant role in organocatalysis [5,6,7,8], crystal engineering and supramolecular chemistry [9,10,11,12,13], materials science [14,15,16,17], stabilization of explosives [18], drug design [19,20], etc. XBs R–X⋯A (X—halogen) are formed between an electron-depleted region on the continuation of the R–X bond and an electron-rich region of another atom or molecule A (Figure 1, bottom). Main geometric criteria of a XB formation are the short X⋯A distance, smaller than the sum of X and A van der Waals radii, and the proximity of the halogen bond angle to the linear [21]. As an electronic criterion, the presence of a critical point of type (3; −1) (bond critical point, BCP) for calculated [22,23,24] or experimentally measured electron density (ED) [25,26,27,28] along the X⋯A bond path in QTAIM analysis (Quantum Theory of Atoms in Molecules [29]) is often used.

Figure 1

Schematic representation of a halogen bond (X—halogen atom, A—nucleophilic site), ED (blue) and ESP (red) distributions along the bond path. The graphical definition of d(EDmin) (blue), d(ESPmin) (red), and Δd (black) are given.

The philicity of halogen bonding participants can be determined using geometric criterion (the angle around electrophilic site is close to 180°) [21], spectral manifestations (for instance, in ultraviolet–visible (UV–vis) spectra [30,31,32,33] or using calculated parameters (electron density deformation (EDD) [34,35,36,37,38,39,40], electron localization function (ELF) distribution [41,42,43,44], natural bond orbital charge transfer [45,46], sums of atomic charges [47], and molecular electrostatic potential (ESP) distribution on the van der Waals surface with ρ(r) = 0.001 e/Bohr3, in representative planes or along the bond path [48,49,50,51,52,53,54,55]).

In [54], it was mentioned for the first time for a series of hydrogen-bonded complexes that the relative arrangement of ED and ESP minima positions along the hydrogen bond path is the same for all complexes. Moreover, the superposition of gradient fields of ED and ESP as discussed for different intermolecular interactions is solids, including the studies based on experimental charge density [56,57]. These observations formed the basis of electronic criterion proposed in Ref. [24], which makes it possible to unambiguously determine the type of electrostatically driven noncovalent bonds. In accordance with electronic criterion, the minimum of ESP along the bond path is located closer to the atom that donates electrons, whereas the minimum of ED is located closer to the atom that delivers its electrophilic site for noncovalent bonding (Figure 1). The latter atom prescribes the name of “-ogen bonds” (e.g., hydrogen, halogen, chalcogen, pnictogen, etc.).

While for R–X⋯A XBs it is often clear which atom acts as an electrophile and which one acts as a nucleophile, the information value of the relative positions of ED and ESP minima seems to hold for less trivial cases of halogen-halogen contacts [58,59,60] and also beyond the halogen bonding in general [61,62].

Up to now, it is not known if there is an information value in the distance between ED and ESP minima, Δd, defined as (see also the graphical definition in Figure 1). A question arises: does a robust correlation exist between Δd and such XB properties as its strength or the X⋯A distance? To the best of our knowledge, it is an original question, not yet discussed in the QTAIM-related literature. t Of fundamental interest, we expect that this sort of information could be useful for studying XBs in solids, where ED and ESP distributions could in principle be measured experimentally [63], whereas the direct experimental evaluation of XB energy is very difficult or even not possible.

In this study, we checked the information value of Δd for a series of 145 halogen-bonded complexes formed by various halogen donors with trimethylphosphine oxide, Me3PO (Figure 2) at the B3LYP/jorge-ATZP level of theory. The Me3PO molecule was chosen as a “standard” halogen acceptor, in a sense following the approach started many years ago by Gutmann and Beckett [64,65] and explored in recent years by a number of authors [66,67,68]. In order to reduce the number of internal degrees of freedom in the electron donor, we have selected Me3PO instead of Et3PO, which was originally used in the Gutmann–Beckett method. As halogen donors the molecules belonging to the following classes of F-, Cl-, Br- and I-containing compounds were considered: halogens, interhalides, oxohalides, pseudohalides, halogenated methane, ethylene, acetylene, benzene, phosgene and their derivatives, as well as thionyl- and sulfurylhalides, sulfur halides and sulfur hypohalites, and several halogenated nitrogen-containing inorganic compounds and some others. This choice of model systems makes it possible to track how a change in the halogen bond acceptor (i.e., its electronic properties as well as belonging to a certain class of chemical compounds) changes the properties of the halogen bond when the halogen bond acceptor Me3PO is fixed. The full list of halogen donors is given in Figure S1 in Supplementary Materials.

Figure 2

Schematic representation of halogen-bonded complexes formed between R–X (X = F, Cl, Br and I) and Me3PO. In blue are given geometric (interatomic distance R(X⋯O)), energetic (complexation enthalpy ΔH), and electronic (ρ stands for electron density, ∇2ρ for its Laplacian, V, G and K stand for local electron potential, kinetic and total energy densities at the bond critical point, BCP) parameters that were correlated with Δd in this work (see Figure 1 for the definition of Δd).

Previously, we have considered a similar but somewhat smaller set of complexes in order to build a correlation between the XB strength/geometry and the changes of two spectral parameters upon complexation: the 31P NMR chemical shift and the ν(P=O) stretching frequency [69] (similar correlations for hydrogen-bonded complexes with Me3PO were recently published as well [70]). We have shown that decent correlations of this kind do exist, and they could be fitted by simple analytical functions. Now, we turn our attention to electronic properties, such as Δd, which is experimentally accessible—if the conditions are right—for single crystal samples. Several XB parameters were considered in this work and correlated with Δd: the complexation enthalpy ΔH, the X⋯O interatomic distance R(X⋯O), and various parameters at XB electron density critical point of type (3; −1): (molecular electrostatic potential ESP(rBCP), electron density ρ(rBCP), Laplacian of electron density ∇2ρ(rBCP), local electron kinetic G(rBCP) and potential V(rBCP) energy densities and total electron energy density K(rBCP)). We also checked if 31P NMR chemical shifts of Me3PO correlate with Δd.

2. Results and Discussion

The optimized geometries of some representative examples of halogen-bonded complexes belonging to different classes of inorganic and organic compounds are shown in Figure 3 and all 145 halogen-bonded complexes are shown in Supplementary Materials (Figure S1). The majority of XBs are linear (or close to linear) and formed along the direction of expected oxygen lone pair localization. The numerical values of relevant parameters are listed in Supplementary Materials: geometric, energetic and spectroscopic parameters (Table S1), and QTAIM electronic parameters (Table S2).

Figure 3

Some examples of optimized geometries of halogen-bonded complexes formed between (a) F2; (b) I2; (c) ClF; (d) BrF; (e) IF; (f) FCN; (g) ClCN; (h) BrCN; (i) ICN; (j) CF3Cl; (k) CF3Br; (l) CF3I; (m) C6F5Cl; (n) C6F5Br; (o) C6F5I and Me3PO considered in this work. Blue dot marks the position of the X⋯O bond critical point.

Figure 4a shows the dependence of distances from oxygen to ESP and ED minima—d(ESPmin) and d(EDmin), respectively—along the XB path as a function of the local electron kinetic energy density at critical point of type (3; −1), G(rBCP) (in kJ/(mol⋅Å3)), taken as a measure of XB strength. In Figure S2 (Supplementary Materials) we show several examples of ED and ESP profiles along the bond path for weak, medium, and strong R–Cl⋯OPMe3 complexes. The d(ESPmin) distance is among the largest ones for free Me3PO (see the black dot in Figure 4a) and falls upon the increase in the XB strength. The d(EDmin) is obviously infinite for free Me3PO, but rapidly—quasi-exponentially—decreases when the XB gets stronger. The d(ESPmin) and d(EDmin) values appear to be almost halogen-independent, except for complexes with F-donors, the d(ESPmin) data points for which deviate consistently from the other sets.

Figure 4

(a) Distances from oxygen atom to minima of molecular electrostatic potential d(ESPmin) (lighter colors) and electron density d(EDmin) (darker colors) along the X⋯O (X = F, Cl, Br and I) bond path; (b) distances between ED and ESP minima Δd as a function of G(rBCP) for a series of halogen-bonded complexes formed between R–X and Me3PO. The solid curves correspond to Equation (2) with fitted parameters listed in Table 1.

However, a closer look at the Δd = d(ESPmin) − d(EDmin) values (Figure 4b) reveals that there are some systematic differences between complexes with Cl-, Br-, and I-donors as well. The general trend is the same for all halogens: the stronger is the XB, the closer are the ESP and ED minima. The Δd values within each series seem to approach a limiting (asymptotic) value for the strongest XBs. A question arises, what are the limiting Δd values for hypothetical strongest possible bonds? As models for such X⋯O bonds, we took cations Me3POX+, the optimized geometries of which together with the X⋯O bond parameters are shown in Figure S3. For these cations, the “asymptotic” G(rBCP) values are the largest and Δd values are the smallest within the respective series of complexes (except for complexes with iodine, where there is a significant scattering of data points for strong XBs). We have added Δd asymptotes as horizontal bold lines in Figure 4b. For F-donors, the Δd even gets slightly negative, suggesting that the electrophile/nucleophile roles are practically swapped for the F–O bond in Me3POF+, as compared to non-covalent R–X⋯OPMe3 complexes (i.e., the F–O bond in Me3POF+ cannot be classified as a halogen bond).

The data sets in the Δd plot versus G(rBCP) for the studied sets of halogen-bonded complexes could be reasonably well fitted, assuming an overall exponential behavior:

Interestingly, only one fitting parameter in Equation (2), , is halogen-dependent (see Table 1). The rate of exponential fall b = 240 kJ/(mol⋅Å3) and the limiting Δd value for G(rBCP) = 0 kJ/(mol⋅Å3), Δd0 = 0.47 Å, seem to be virtually the same for all halogens. At the moment, it is not clear if these values for Δd0 and b reflect some property of Me3PO or a property of XBs in general. The result of the fitting is added to Figure 4b as solid lines. Another noteworthy thing is that the fitting with Equation (2) does not reproduce the asymptotic values for Me3POX+ at G(rBCP) → ∞, though the relative position of the fitting asymptotes for F, Cl, Br and I is the same as for the corresponding cations Me3POX+. Solving Equation (2) for G(rBCP) one obtains

Table 1

Fitting parameters , Δd0 and b for Equation (2) and proportionality coefficient k between G(rBCP) and complexation energy ΔE (Equation (4) [69]) for the data sets plotted in Figure 4b.

Halogen Donor	a_X, Å	Δd0, Å	b, kJ/(mol⋅Å3)	k, Å3
F	0.09	0.47	240	0.18
Cl	0.15	0.47	240	0.47
Br	0.18	0.47	240	0.57
I	0.19	0.47	240	0.74

It should be noted that Equation (3) makes sense only when the logarithm is defined (Δd > aX) and positive (Δd < Δd0). This means that Equation (3) might become inapplicable for extremely weak and extremely strong XBs. Subsequently, the G(rBCP) values could be converted into the XB “complexation energies” ΔE (the difference of total electronic energies of the complex and its isolated relaxed constituents) using previously published halogen-dependent proportionality coefficients k (see the last column in Table 1; the values were taken from Table 2 in Ref. [69]),

At this point, it is worth making a brief comment concerning various quantitative models based on the properties of electron density at the bond critical points, which were previously successfully used for characterization of the energy for various types of noncovalent bonds [70,71,72,73,74,75]. Figure S4 shows the correlation between the complexation enthalpy ΔH and G(rBCP). Note that there is a significant scattering of the data points due to the fact that a number of complexes within the studied series are held not only by an XB, but also by other noncovalent interactions, most prominently weak hydrogen bonds between electronegative atoms of the halogen donor (including the halogen itself) and the methyl protons of the Me3PO moiety (Figure 5). For completeness of the subject, in Figures S5–S8 we show also correlations of electron density ρ(rBCP), Laplacian of electron density ∇2ρ(rBCP) and total electron energy density K(rBCP) = G(rBCP) + V(rBCP) with d(ESPmin), d(EDmin), Δd, and Rnorm, as well as G(rBCP) correlation with V(rBCP), ρ(rBCP), ∇2ρ(rBCP) and ESP(rBCP). In turn, Figure S9 shows ΔH, ρ(rBCP), ∇2ρ(rBCP) and ESP(rBCP) dependence on Δd.

Figure 5

Examples of additional non-covalent interactions (red dashed lines) present in some halogen-bonded (black dashed line) R–X⋯OPMe3 complexes: hydrogen bonds between methyl protons of the Me3PO moiety and electronegative atoms of the halogen donor; (а) two H-bonds with two proton acceptors; (b) two H-bonds with one proton acceptor; (c) two H-bonds with the halogen donating atom.

Finally, the Δd values are plotted in Figure 6 as a function of Rnorm, defined as X⋯O interatomic distance, normalized by the sum of van der Waals radii (RvdW) of X (X = F, Cl, Br and I) and O: where the following values of the van der Waals radii were used: 1.52 Å (O), 1.47 Å (F), 1.75 Å (Cl), 1.85 Å (Br), 1.98 Å (I) [76]. The usage of normalized and unitless Rnorm, instead of direct interatomic distances X⋯O, allows one—at least in principle—to compare XBs with participation of different halogens (see the dependence of absolute R(X⋯O) distances and Rnorm on G(rBCP) in Figure S10). The data sets plotted in Figure 6 indicate slightly non-linear dependencies of Δd on Rnorm for each halogen, but generally show the same trends as Figure 4b, because the terms halogen bond strength and shortness are almost interchangeable: stronger/shorter XBs are characterized by smaller Δd values.

Figure 6

Distances between ED and ESP minima Δd as a function of normalized interatomic distances Rnorm between X (X = F, Cl, Br and I) and O (see definition in Equation (5)).

To summarize the results so far, there is indeed a correlation between Δd and XB energy and length. The exponential fall of Δd upon strengthening/shortening of the XB is dependent on the type of halogen and most probably on the type of the electron-donating atom as well (oxygen in our case), but for a fixed pair of atoms the correlations seem to be largely independent on the substituents, which makes Δd a promising tool in characterization of XBs and—speculatively—other types of non-covalent interactions as well.

Now we turn attention to the spectroscopic manifestation of halogen bonding within the studied series of complexes. In Ref. [69], some of us have demonstrated that isotropic 31P NMR chemical shift of Me3PO could be used as a spectroscopic marker for the halogen-donating ability of a probed molecule or, in other words, as a measure of the XB strength (length) in a R–X⋯OPMe3 complex. For the data sets presented in this work, reasonable correlations between the change of the 31P NMR chemical shift upon complexation Δδ31P and the complexation enthalpy ΔH exist for chlorine- and bromine-containing complexes (Figure S11a). For fluorine- and iodine-containing ones, there is a large scattering of the data points, which in case of complexes with iodine is likely to be—at least partially—due to the presence of additional non-covalent interactions and, consequently, the influence of several competing factors (e.g., presence of additional interactions) on the electron shells of the phosphorus atom. Indeed, the Δδ31P correlation with G(rBCP) values (Figure S11b) is noticeably better, though with significant residual scattering still (it should be mentioned, that the shape of Δδ31P(Δd) remains the same upon the change of a basis set).

Figure 7 shows the correlation of Δδ31P with Δd, which has a similar degree of data point scattering as Figure S11b. Because of that, we have added to Figure 7 two trend curves for chlorine- and bromine-containing complexes only. These curves serve only as guides for the eye, as it seems premature to propose a functional fit. Still, one could confirm that 31P NMR chemical shift sensitively reflect the changes in the electronic structure of the complexes—including the parameter in focus of this work, Δd—and could serve for the characterization of XBs with phosphine oxides.

Figure 7

The correlation of Δδ31P on Δd for a series of R–X⋯OPMe3 complexes. The solid lines are guides for the eye added for Cl- and Br-bonded complexes.

3. Computational Methods

The full geometry optimization, harmonic vibrational frequencies and NMR calculations were performed for studied complexes in vacuum using B3LYP functional [77] and nonrelativistic all-electron augmented triple-zeta valence quality basis set with polarization functions jorge-ATZP [78,79,80] (adopted from the Basis Set Exchange site [81]) using Gaussian16 software package [82]. Geometry optimization was performed using the default for Gaussian16 Berny algorithm without any geometry restrictions and standard convergence criteria (maximum/RMS forces and displacements of atoms smaller than 0.000450 a.u./0.000300 a.u. and 0.001800 a.u./0.001200 a.u., respectively). The optimized geometries were checked for the absence of imaginary vibrational frequencies. The presence of halogen bonds in the calculated complexes was confirmed by the criteria for halogen bond formation, according to the IUPAC recommendation [21]. For each halogen donor, only one optimized structure of its complex with phosphine oxide was considered. For a subset of complexes, we have checked that variations in the initial geometries lead to the same optimized structures. Because of this, we believe that they are true global minima. However, even if it is not the case, any local true minimum containing a halogen bond and satisfying the Virial theorem would be legitimate structure to perform topological analysis of electron density and use the resulting data to construct correlations.

The complexation enthalpy ΔH was calculated at 298.15 K as the enthalpy required to separate the interacting molecules at infinite distance (including the relaxation of monomers). Note that in this way the ΔH values include contributions not only from the XB, but also from any other non-covalent interactions which might be present between monomers of Me3PO and XB donor.

Isotropic NMR shielding constants σ were calculated using the gauge-independent atomic orbital (GIAO) approach [83] and converted into changes of chemical shifts upon complexation as Δδ31P = σfree − σ, where σfree is the 31P nuclear shielding constant for an isolated Me3PO molecule.

The topological analysis of ED and ESP along the XB path was carried out within the framework of QTAIM methodology from the wave function files using MultiWFN software (http://sobereva.com/multiwfn/ accessed on 22 June 2022; Beijing Kein Research Center for Natural Sciences; Beijing, China; version 3.3.8) [84]. The ED and ESP minima positions from the oxygen atom along the XB path were determined with the following path searching parameters: maximum number of points of a path is 100000, with stepsize of 0.0005 Bohr and generation stop threshold of 0.005 Bohr distance to any critical point. This takes into account that for all studied complexes d(ESPmin) < d(EDmin), the definition given in Equation (1) makes all Δd values positive.

Visualization of the studied complexes was performed using the Chemcraft software (available online at www.chemcraftprog.com accessed on 22 June 2022; version 1.8) [85].

All of the proposed correlation functions were fitted by Levenberg–Marquardt algorithm and visualized using Origin software (OriginLab Corporation, Northampton, MA, USA) [86]. The complexes with X⋯H (X = N, O, F, Cl, Br or I; H are protons of methyl groups) contacts shorter than the sum of Bondi’s van der Waals radii of corresponding atoms are marked in red in the Supplementary Materials (Figure S1 and Tables S1 and S2).

4. Conclusions

The summary and the main conclusion of this work are rather concise. For a homologous series of 145 halogen-bonded complexes with the general formula R–X⋯OPMe3 (X = F, Cl, Br, and I), we showed that the XB strength/energy correlate well with the distance between ED and ESP minima along the X⋯O bond path, Δd (see Equation (2); Figure 4b). The maximum Δd value (for weakest XBs) as well as the exponent of the fall of the Δd dependence on G(rBCP) are halogen-independent, whereas the limiting values for strongest XBs are halogen-dependent. One could expect that there is also a dependence on the type of electron-donating atom (here, oxygen), though this and other limits of applicability of the proposed correlations remain to be studied. Quite likely, the most robust conclusion one could make is that for a pair of homologous non-covalently bound complexes the one with larger Δd is weaker. The numerical value of Δd might appear to be useful in analysis of experimental high-resolution X-ray data on ED of halogen-bonded single crystals. We also propose Δd as a new tool for routine QTAIM analysis of the energies of XBs.