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.

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 liter- 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(ED min ) (blue), d(ESP min ) (red), and ∆d (black) are given.
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 halogenbonded complexes formed by various halogen donors with trimethylphosphine oxide, Me 3 PO ( Figure 2) at the B3LYP/jorge-ATZP level of theory. The Me 3 PO 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 Me 3 PO instead of Et 3 PO, 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 prop- erties as well as belonging to a certain class of chemical compounds) changes the properties of the halogen bond when the halogen bond acceptor Me 3 PO is fixed. The full list of halogen donors is given in Figure S1 in Supplementary Materials. 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 31 P 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 31 P NMR chemical shifts of Me3PO correlate with Δd. 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 31 P NMR chemical shift and the ν(P=O) stretching frequency [69] (similar correlations for hydrogen-bonded complexes with Me 3 PO 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(r BCP ), electron density ρ(r BCP ), Laplacian of electron density ∇ 2 ρ(r BCP ), local electron kinetic G(r BCP ) and potential V(r BCP ) energy densities and total electron energy density K(r BCP )). We also checked if 31 P NMR chemical shifts of Me 3 PO correlate with ∆d.

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 4a shows the dependence of distances from oxygen to ESP and ED minimad(ESP min ) and d(ED min ), respectively-along the XB path as a function of the local electron kinetic energy density at critical point of type (3; −1), G(r BCP ) (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· · · OPMe 3 complexes. The d(ESP min ) distance is among the largest ones for free Me 3 PO (see the black dot in Figure 4a) and falls upon the increase in the XB strength. The d(ED min ) is obviously infinite for free Me 3 PO, but rapidly-quasi-exponentially-decreases when the XB gets stronger. The d(ESP min ) and d(ED min ) values appear to be almost halogenindependent, except for complexes with F-donors, the d(ESP min ) data points for which deviate consistently from the other sets. The optimized geometries of some representative examples of halogen-bonded com-plexes 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).  However, a closer look at the ∆d = d(ESP min ) − d(ED min ) 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 Me 3 POX + , the optimized geometries of which together with the X· · · O bond parameters are shown in Figure S3. For these cations, the "asymptotic" G(r BCP ) 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 Me 3 POF + , as compared to non-covalent R-X· · · OPMe 3 complexes (i.e., the F-O bond in Me 3 POF + cannot be classified as a halogen bond). 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.  Table 1. Table 1. Fitting parameters X , Δ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.  (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).

Halogen
We have added Δd asymptotes as horizontal bold lines in Figure 4b. For F-donors, the Δd  Table 1. Table 1. Fitting parameters a X , ∆d 0 and b for Equation (2) and proportionality coefficient k between G(r BCP ) and complexation energy ∆E (Equation (4) [69]) for the data sets plotted in Figure 4b. The data sets in the ∆d plot versus G(r BCP ) for the studied sets of halogen-bonded complexes could be reasonably well fitted, assuming an overall exponential behavior:

Halogen Donor
Interestingly, only one fitting parameter in Equation (2), a X , is halogen-dependent (see Table 1). The rate of exponential fall b = 240 kJ/(mol·Å 3 ) and the limiting ∆d value for G(r BCP ) = 0 kJ/(mol·Å 3 ), ∆d 0 = 0.47 Å, seem to be virtually the same for all halogens. At the moment, it is not clear if these values for ∆d 0 and b reflect some property of Me 3 PO 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 Me 3 POX + at G(r BCP ) → ∞, though the relative position of the fitting asymptotes for F, Cl, Br and I is the same as for the corresponding cations Me 3 POX + . Solving Equation (2) for G(r BCP ) one obtains It should be noted that Equation (3) makes sense only when the logarithm is defined (∆d > a X ) and positive (∆d < ∆d 0 ). This means that Equation (3) might become inapplicable for extremely weak and extremely strong XBs. Subsequently, the G(r BCP ) 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(r BCP ). 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 Me 3 PO moiety ( Figure 5). For completeness of the subject, in Figures  S5-S8 we show also correlations of electron density ρ(r BCP ), Laplacian of electron density ∇ 2 ρ(r BCP ) and total electron energy density K(r BCP ) = G(r BCP ) + V(r BCP ) with d(ESP min ), d(ED min ), ∆d, and R norm , as well as G(r BCP ) correlation with V(r BCP ), ρ(r BCP ), ∇ 2 ρ(r BCP ) and ESP(r BCP ). In turn, Figure S9 shows ∆H, ρ(r BCP ), ∇ 2 ρ(r BCP ) and ESP(r BCP ) dependence on ∆d.
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 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 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. Finally, the ∆d values are plotted in Figure 6 as a function of R norm , defined as X· · · O interatomic distance, normalized by the sum of van der Waals radii (R vdW ) 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 R norm, 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 R norm on G(r BCP ) in Figure S10). The data sets plotted in Figure 6 indicate slightly non-linear dependencies of ∆d on R norm 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.
of direct interatomic distances XO, allows one-at least in principle-to compare X with participation of different halogens (see the dependence of absolute R(XO) distan and Rnorm on G(rBCP) in Figure S10). The data sets plotted in Figure 6 indicate slightly n linear dependencies of Δd on Rnorm for each halogen, but generally show the same tre as Figure 4b, because the terms halogen bond strength and shortness are almost in changeable: stronger/shorter XBs are characterized by smaller Δd values.  (5)).
To summarize the results so far, there is indeed a correlation between Δd and energy and length. The exponential fall of Δd upon strengthening/shortening of the X dependent on the type of halogen and most probably on the type of the electron-donat atom as well (oxygen in our case), but for a fixed pair of atoms the correlations seem to largely independent on the substituents, which makes Δd a promising tool in charact zation of XBs and-speculatively-other types of non-covalent interactions as well.
Now we turn attention to the spectroscopic manifestation of halogen bonding wit the studied series of complexes. In Ref. [69], some of us have demonstrated that isotro 31 P NMR chemical shift of Me3PO could be used as a spectroscopic marker for the halog donating ability of a probed molecule or, in other words, as a measure of the XB stren (length) in a R-XOPMe3 complex. For the data sets presented in this work, reasona correlations between the change of the 31 P NMR chemical shift upon complexation δ and the complexation enthalpy H exist for chlorine-and bromine-containing comple ( Figure S11a). For fluorine-and iodine-containing ones, there is a large scattering of data points, which in case of complexes with iodine is likely to be-at least partially-d to the presence of additional non-covalent interactions and, consequently, the influe of several competing factors (e.g., presence of additional interactions) on the elect  (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 31 P NMR chemical shift of Me 3 PO could be used as a spectroscopic marker for the halogendonating ability of a probed molecule or, in other words, as a measure of the XB strength (length) in a R-X· · · OPMe 3 complex. For the data sets presented in this work, reasonable correlations between the change of the 31 P NMR chemical shift upon complexation ∆δ 31 P 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 ∆δ 31 P correlation with G(r BCP ) values ( Figure S11b) is noticeably better, though with significant residual scattering still (it should be mentioned, that the shape of ∆δ 31 P(∆d) remains the same upon the change of a basis set). Figure 7 shows the correlation of ∆δ 31 P 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 31 P 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. 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 31 P 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.

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.

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 Me 3 PO 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 ∆δ 31 P = σ free − σ, where σ free is the 31 P nuclear shielding constant for an isolated Me 3 PO 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(ESP min ) < d(ED min ), 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).

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· · · OPMe 3 (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(r BCP ) 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.