Strength of the [Z–I···Hal]− and [Z–Hal···I]− Halogen Bonds: Electron Density Properties and Halogen Bond Length as Estimators of Interaction Energy

Bond energy is the main characteristic of chemical bonds in general and of non-covalent interactions in particular. Simple methods of express estimates of the interaction energy, Eint, using relationships between Eint and a property which is easily accessible from experiment is of great importance for the characterization of non-covalent interactions. In this work, practically important relationships between Eint and electron density, its Laplacian, curvature, potential, kinetic, and total energy densities at the bond critical point as well as bond length were derived for the structures of the [Z–I···Hal]− and [Z–Hal···I]− types bearing halogen bonds and involving iodine as interacting atom(s) (totally 412 structures). The mean absolute deviations for the correlations found were 2.06–4.76 kcal/mol.

Properties of functional materials and supramolecular structures formed by halogen bonding are controlled by its strength. Therefore, halogen bond energy is the most important characteristic of these interactions. However, determination of the bond energy (E b ) is a difficult task. Experimental methods (e.g., mass spectrometry, vibrational or NMR spectroscopy, calorimetry, vapor density measurements) are usually associated with complex technical procedures and can be applied to a very limited number of structures. The direct theoretical calculations of E b for the A···B bond as the energy difference between the structure A···B and the isolated molecules A and B may be successfully applied to intermolecular non-covalent interactions in the gas phase. However, in the condensed phase, molecules are usually bound with each other by a network of several non-covalent interactions, and an adequate fragmentation which affects only the bond of interest is often impossible. In such a situation, an approximation of E b through other parameters easily accessible from experiment becomes very important [51][52][53][54][55][56][57][58][59][60]. Among these parameters, bond distance [61][62][63][64][65][66][67][68][69][70][71][72][73][74][75][76][77][78], electron density properties [61,, and force constants [57,73,114,115] are the most common. In fact, the application of E b -property relationships is often the only possibility to estimate the energy of non-covalent interactions in the condensed phase.
A broad practical utilization of the E b -property relationships for the estimates of E b started with a publication by Espinosa, Molins, and Lecomte [82], where a simple correlation between the bond energy and the potential energy density (V b ) at the bond critical point (BCP) was obtained for the X-H···O hydrogen bonds (X = C, N, O): where E int is the interaction energy of the A···B bond defined as E int = E A···B − E A − E B (A and B are fragments with unrelaxed geometry corresponding to the structure A···B). Later, several other equations for the approximation of E int through the electron density properties were published for hydrogen and halogen bonds ( Table 1 and Table S1). Meanwhile, it was shown that these relationships are not universal and may be applied only to the type of interaction for which they were derived [116][117][118]. Table 1. Relationships between interaction energy (E int , kcal/mol) and potential or kinetic energy density (V b or G b , kcal/(mol•bohr 3 )).
Recently, the author started a project to establish practically useful correlations between the interaction energy and properties which could be easily determined from experiment for halogen bonds of various types, including those formed by a halide anion. In the previous work [116], such relationships were found for structures of the [(A) n Z-Y···X] − type (X, Y = F, Cl, Br; Z = F, Cl, Br, C, N, O, H, S, P, Si, B; 441 structures in total). The aim of this work was to include in the consideration the structures bearing the interacting iodine atom(s), thus completing the whole series of halogen bonds with a halide anion. In this work, structures of the [(A) n Z-I···X] − and [(A) n Z-Y···I] − types (X = F, Cl, Br, I; Y = Cl, Br, I; Z = F, Cl, Br, I, C, N, O, H, S, P, Si, B) were calculated. The correlations between E int and the electron density (ρ b ), its Laplacian (∇ 2 ρ b ), the curvature of ρ(r) which is parallel to the bond path direction (positive) (λ ||,b ), the potential, kinetic, and total energy densities (V b , G b , and H b ) at BCP, and the halogen bond length (d Y···X ) were established. Two points should be mentioned here. First, the interactions with the iodine atom as a halogen bond donor are the most important among all halogen bonds since the existence of a prominent σ-hole at the iodine atom provides particularly high stability of these interactions. Second, reliable relationships can be obtained only for sufficiently large sets of the structures [(A) n Z-Y···X] − which include different types of second-order atoms Z and third-order groups A, otherwise, the relationships found may suffer a significant overfitting. Therefore, a large statistically significant set of 412 structures bearing twelve different types of the Z atom was used in this work.

Computational Details
Full geometry optimization of the main set of structures was carried out at the density functional theory (DFT) level by using the M06-2X functional [121] with the help of the Gaussian 09 [122] program package applying tight optimization criteria and an ultrafine integration grid. Cartesian d and f basis functions (6d, 10f) were used in all calculations. The ADZP basis set taken from the Basis Set Exchange library [123] was applied for geometry optimization. The energies of the equilibrium structures were refined by the second-order scalar relativistic correction using the Douglas−Kroll−Hess method (DKH) and single-point calculations with the ADZP-DKH basis set [123,124]. The ADZP-DKH basis set was constructed from the DZP-DKH basis set [125][126][127] by addition of diffuse functions taken from the ADZP basis set.
Similar level of theory was applied in the previous work for the analysis of correlations for the [(A) n Z-Y···X] − structures (X, Y = F, Cl, Br) [116]. It was shown by Kozuch and Martin [128] for a set of 51 structures bearing halogen bonds that the interaction energies calculated with the M06-2X functional correlate well with the reference CCSD(T) level providing the root-mean-square deviation, the mean signed error, and the maximum error of 0.43, 0.01, and 1.58 kcal/mol, respectively. It was also shown [116] that the E int (V b ) correlations for the [(A) n Z-Y···X] − structures (X, Y = F, Cl, Br) obtained at the M06-2X level of theory have similar parameters as those for the MP4, CCSD, and CCSD(T) methods. On the other hand, it was found [116] that the simple double-zeta quality basis set 6-31+G* which is similar to ADZP provides the interaction energy values and parameters of the E int (V b ) correlations close to those obtained at the much more extended triple-zeta basis set 6-311++G(3df,3pd).
The applicability of the M06-2X/ADZP-DKH level of theory for the reproduction of experimental interaction energies and electron densities is additionally discussed below in the Results. For this analysis, the interaction energies and electron density values were calculated at the CCSD, CCSD(T), and PBE0-D3BJ levels of theory. The ADZP-DKH, 6-31+G*, 6-311+G**, and aug-cc-pVTZ basis sets were applied.
The Hessian matrix was calculated for all optimized structures to prove the location of correct minima. No symmetry operations were applied during the calculations. The stability test was performed, and stable solutions were achieved for all structures using the keyword STABLE(OPT).
The interaction energies for the Y . . . X bond in [(A) n Z-Y···X] − were calculated using Equation (1): where BSSE is a basis set superposition error estimated using the counterpoise (CP) method [129,130] with the ADZP basis set. Geometries of the (A) n Z-Y fragments were unrelaxed and corresponded to those in [(A) n Z-Y···X] − . The topological analysis of the electron density distribution was performed with the help of the Atoms in Molecules (AIM) method developed by Bader [131] using the AIMAll program [132].  For Z = C and N, different orbital hybridizations of these atoms were considered. For Z = S and P, different oxidation states of these atoms were considered. The groups A (totally 66 groups) vary from the electron donor to the electron acceptor ones providing a broad span of interaction energies for the halogen bond. The complete list of the calculated structures is given in Table S2.

Computational Models
The

Test of the Computational Method
Since the main goal of this work is to establish relationships between the interaction energy and the experimentally accessible properties such as electron density and derivatives and the length of the halogen bond, it is essential, first, to test the ability of the computational method used to reasonably reproduce these experimental properties. Results of such analysis are presented in this section.
Interaction energy. It is a usual practice to test reliability of a computational method for the estimates of interaction energies against a reference level of theory, such as coupled cluster (CC) approaches. Meanwhile, even the CCSD(T) level often considered as the "gold standard" may be associated with some degree of uncertainty, for instance, due to issues related with the basis set superposition error and its correction [133]. Therefore, when it is possible, experimental data should always be considered as the ultimate reference for method benchmarking.
Here, interaction energies were calculated for several structures bearing the hydrogen or halogen bonds with a halide anion for which experimental data are available (i.e., H 3 C-H···F − , [F-H-F] − , HO-H···Cl − , HO-H···I − , H 3 C-H···I − , and [I-I-I] − ). As can be seen in Table 2, the M06-2X functional with the moderate double-zeta quality basis sets 6-31+G* or ADZP-DKH and the CP correction provide good agreement with the experimental bond energies with the mean absolute deviation (MAD) of 1.2 kcal/mol for the whole set of these structures and of 0.6 kcal/mol excluding [F-H-F] − . The experimental error is in the range of 0.1-1.6 kcal/mol, and theoretical M06-2X values fall within the 3σ interval for four of six structures. The interaction energies calculated at the CCSD and CCSD(T) levels were strongly underestimated when the CP correction was applied. The CCSD method performed worse than M06-2X for H 3 C-H···F − and HO-H···Cl − even with the aug-cc-pVTZ basis set, and only the CCSD(T) level with the triple-zeta basis set and without the CP correction provided better results than the DFT approach (MAD = 0.3 kcal/mol). Thus, considering that typical MAD values for the E int property relationships found for the structures [(A) n Z-Hal 1 ···Hal 2 ] − are 2-4 kcal/mol (see [116] and results of this work below), the M06-2X/6-31+G*(or ADZP-DKH) level of theory is reliable for the analysis of interaction energies. Table 2. Enthalpies of formation of hydrogen and halogen bonds (in kcal/mol).
. k Adiabatic interaction energy in terms of E int . l Ref. [137]. m Enthalpy at 0 K.

Electron density.
Recently, it has been shown that the M06-2X functional is not sufficiently good to reproduce the experimental electron density properties for isolated atoms and atomic cations [138]. Here, the ability of this functional to reproduce the experimental ρ b values was tested for the BCPs corresponding to the halogen bonds in eleven structures ( Table 3). The calculated structures corresponded to the experimental X-ray ones.
The results indicate that the M06-2X/ADZP-DKH method adequately reproduces the experimental electron density at BCPs in these structures with the MAD and RMSD values of 0.006 and 0.009 e/A 3 , respectively. The slope of the correlation ρ b,theor against the experimental values ρ b,exp was quite close to unity (1.07, Figure S1). The PBE0 functional found by Medvedev et al. [138] as one of the best to reproduce the electron density properties in isolated atoms and atomic cations demonstrates the performance similar to M06-2X for halogen bonds (MAD = 0.006 e/A 3 , RMSD = 0.009 e/A 3 , slope, 1.12). Correspondingly, the E int (V b ) relationships obtained for the small set of ten structures of the [(A) n Z-I···F] − type were similar for both M06-2X and PBE0-D3BJ functionals ( Figure S2). The deviation of E int estimated by these two methods was 2.27 kcal/mol at V b = 20 kcal/(mol•bohr 3 ) and it was 3.96 kcal/mol at V b = 100 kcal/(mol•bohr 3 ). These deviations were lower than MAD for the whole set of seventy [(A) n Z-I···F] − structures (4.76 kcal/mol, see below). All these results indicate that the relationships found in this work at the M06-2X/ADZP-DKH level are reliable for the determination of E int using experimental electron density properties.   An interesting relationship distinct from all other bonds with halide anions was obtained for the [(A) n Z-Cl···I] − series ( Figure 2). These structures could not be treated by a single linear or quadratic E int (V b ), E int (G b ), or E int (ρ b ) function. Instead, the horizontal branch is clearly visible in the range of E int of 10-20 kcal/mol. The interaction energy could be reasonably approximated by polynomial fourth-order functions, but such an approximation is usually associated with significant overfitting and cannot be recommended for practical use. The reasons of this peculiar situation are discussed below. An interesting relationship distinct from all other bonds with halide anions was obtained for the [(A)nZ-Cl···I] − series ( Figure 2). These structures could not be treated by a single linear or quadratic Eint(Vb), Eint(Gb), or Eint(ρb) function. Instead, the horizontal branch is clearly visible in the range of Eint of 10-20 kcal/mol. The interaction energy could be reasonably approximated by polynomial fourth-order functions, but such an approximation is usually associated with significant overfitting and cannot be recommended for practical use. The reasons of this peculiar situation are discussed below. Analysis of the relationships obtained for the structures [(A)nZ-Y···X] − in this and previous [116] works allows making the following additional conclusions and generalizations. First, all three estimators, Vb, Gb, and ρb, behave similarly (Table 4)  Third, the slope of the relationships is strongly affected by the nature of the Y and X atoms. For instance, for the V b estimator, the lowest slope (0.69-0.84) and the highest negative intercept (−3.68 ÷ −6.50) were found for X = F. The highest slopes were detected for [(A) n Z-I···I] − (2.34) as well as for the [(A) n Z-Y···X] − (Y = Br, X = I; Y = I, X = Br, X = Y = Br) series (1.55-2.12). The series with (X or/and Y) = Cl but (X and Y) = F exhibited an intermediate slope (1.26-1.46). Such a behavior can be qualitatively interpreted by the different properties of the halogen atoms forming the Y···X bond. First of all, the electronegativity of the X atom decreases along the series F > Cl > Br > I. There is a clear trend between the slope of the E int (V b ) dependence and the difference of the electronegativities of X and Y ( Figure S3). Thus, the higher the electronegativity difference, the lower the E int at a given energy density at BCP. This conclusion becomes understandable if we recall that ionic bonds have high interaction energies but low electron densities at BCP. The second factor is an enhancement of polarizability, diffuse character of electron shells, and the role of dispersion forces from X = F to I. Electronic effects from the (A) n Z group on the Y···X bond are more pronounced if the electron shells of X are more diffuse, i.e., from X = F to I. Since the higher dispersion contribution and more diffuse electron shells are associated with lower concentration of electron or energy density at BCP, E int values are higher for most electron acceptor groups (A) n Z at a given value of E int in the case of heavier atom X.  Fourth, some series form distinct groups which can be approximated by a single relationship of a reasonable quality. There are three such groups when considering the E int (V b ) dependence, i.e., dependences are characterized by a higher dispersion between series. They are more sensitive to the nature of the Y and X atoms and, therefore, such a grouping is not efficient for those relationships.

Series
Fifth, the relationships for any series of the [(A) n Z-Y···X] − structures were significantly different from those obtained previously for the hydrogen bonds X-H . . . O and FH . . . FR [82,83] and for the neutral halogen bonds [117,119,120] (Table 1).

The E int (∇ 2 ρ b ) Relationship
The Laplacian of electron density may be used as an estimator of E int only for the

The E int (λ ||,b ) Relationship
The curvature of the electron density distribution can serve for the estimates of E int only for three series calculated in this work, i.

The E int (H b ) Relationship
The quality of the total energy density as an estimator of E int depends on the nature of the Y and X atoms. It cannot be used for the [(A) n Z-I···F] − series due to high dispersion of the E int (H b ) function at lower E int and its low sensitivity at higher E int ( Figure S6). For other series bearing the iodine Y and/or X atom, the E int (H b ) relationship is not sensitive for the weakest interactions with E int < 10 kcal/mol. However, for stronger interactions, good linear dependences were observed (R 2 = 0. 95 (the fittings were performed for points corresponding to −Eint > 10 kcal/mol, data for the non-iodine structure were taken from [116]).

The Eint(dY…X) Relationship
All series except for [(A)nZ-Cl···I] − demonstrate reasonable exponential relationships between interaction energy and the length of the halogen bond ( Figure 6). The quality of these relationships is similar to those with Vb, Gb, and ρb. There are only two couples of series which can be approximated by single correlations, i.

The E int (d Y . . . X ) Relationship
All series except for [(A) n Z-Cl···I] − demonstrate reasonable exponential relationships between interaction energy and the length of the halogen bond ( Figure 6). The quality of these relationships is similar to those with V b , G b , and ρ b . There are only two couples of series which can be approximated by single correlations, i.

The [(A) n Z-Cl···I] − Series
As was mentioned above, the E int -property relationships for the [(A) n Z-Cl···I] − series have a peculiar character. The more detailed analysis shows that such behavior is due to the tremendous effect of the second-order atom Z and the third-order groups A.

The [(A)nZ-Cl···I] − Series
As was mentioned above, the Eint-property relationships for the [(A)nZ-Cl···I] − series have a peculiar character. The more detailed analysis shows that such behavior is due to the tremendous effect of the second-order atom Z and the third-order groups A.  Third, the [(A) n N-Cl···I] − structures could be divided into two distinct groups (Figure 7), i.e., structures with a weaker Cl···I interaction (E int < 10 kcal/mol, aliphatic and monoaryl amines, methylidene imine and nitroso compounds, Group I) and those with a stronger Cl···I interaction (E int > 10 kcal/mol, diarylamines, difluoro(methylidene) imine and nitro compounds, Group II). The trend slope for these two groups is similar (2.1) but the intercept is very different (−5.7 and −30.4, respectively). The different behavior of these two groups is apparently associated with a different nature of the (A) n N part of the molecule. In Group I, (A) n N exhibits either electron donor or weak electron acceptor properties. In Group II, (A) n N is a strong electron acceptor. Furthermore, the structures of Group II have more extended conjugation systems compared to those of Group I. All these relationships considered together provide such a peculiar trend for the whole [(A) n Z-Cl···I] − series. Interestingly, this effect of Z and A is much smaller for other series making possible the existence of linear correlations. Similar behavior was also found for the G b and ρ b estimators.
The great effect of the Z atom and A group on E int does not permit an approximation of the latter via

Discussion
The estimate of E int for the condensed phases is a much more difficult task than that for the gas phase. Both periodic solid state and cluster approaches often do not allow direct estimates of E int as an energy difference because fragments may have multiple mutual weak intermolecular interactions. When it is possible, the adequate computational model is usually exceedingly large. For instance, accurate calculations of E int for the A···B bond in a solid-state structure within the cluster approach requires explicit consideration of all intermolecular interactions which molecules A and B are involved in. An example of an adequate cluster for the direct estimates of E int between the simple fragments I-C 4 F 8 -I and Cl − in the X-ray structure ACIPOU is shown in Figure S7 and it has 386 atoms (this model is composed of the binuclear cluster I-C 4 F 8 -I···Cl − surrounded by all molecules which form short contacts with the I-C 4 F 8 -I or Cl − fragments).
Considering this situation, the practice of utilization of small bimolecular clusters which include only two molecules forming a contact of interest becomes increasingly popular. In other words, there are multiple attempts to approximate the properties of solid-state structures by gas phase calculations.
There are two points which invalidate this practice. First, such a bimolecular model ignores other intermolecular interactions in which both considered molecules participate. These secondary interactions may significantly affect the energy of the contact under study. For instance, E int calculated for the I-CF 2 CF 2 -I···Cl − halogen bond in the X-ray structure ACIPIO using the bimolecular approach is 20.07 kcal/mol. However, when four molecules are considered (three I-CF 2 CF 2 -I and one Cl − ion, Figure S8), E int for the I···Cl bond becomes 14.75 kcal/mol (see caption to Figure S8 for details) (this structure was taken as one of the examples bearing a short contact between the I atom and the Cl − ion and where a relatively small cluster may be selected for reliable calculations of the interaction energy).
The second point is that geometry of bimolecular clusters is usually not optimized but taken from the X-ray structures since optimization often leads to the collapse of the structure compared to X-ray. In such calculations, the geometry considered is not equilibrium and does not belong to a minimum energy path for either adiabatic or vertical bond cleavage, and, hence, the resulting E int has little physical meaning.
Application of the correlations established in this work for the rapid estimates of E int in the condensed phase is free from the difficulties discussed above since it requires a single parameter (electron density or halogen bond length) directly obtained from experiment.

Final Remarks
In this work, correlations between the halogen bond interaction energy and experimentally accessible properties, i.e., electron density and bond length, were established for the series of structures [(A) n Z-I···Hal] − and [(A) n Z-Hal···I] − bearing the iodine atom (412 structures in total). Thus, together with the results published previously [116], this work finalizes the analysis of halogen bonds formed by a halide anion. Obtained correlations are of practical use for rapid estimates of the halogen bond energy in systems with a network of multiple weak interactions. Several concluding remarks and generalizations can be made.
First, among all structures [(A) n Z-Y···X] − , those with X = F demonstrate quite poor approximations of E int with MAD of 3.29-4.76 kcal/mol. For other structures, except [(A) n Z-Cl···I] − , more reasonable correlations were obtained (MAD = 2.06-3.01 kcal/mol). Meanwhile, the approximations of E int for halogen bonds formed by halide anions are worse than those obtained for homohalogen bonds formed by two neutral fragments [117].
Second, the most difficult case is the structures of the [(A) n Z-Cl···I] − type. In fact, only two parameters were found to be able to adequately estimate E int for this bond, i.e., H b (for strong bonds with negative H b values) and d Y . . . X .
Third, there is no unique equation to approximate E int for all structures [(A) n Z-Y···X] − . For the electron density-based properties, the highest slopes (linear correlations) or curvatures (non-linear correlations) were found for the bonds formed by heavy halogens (I and Br) while the lowest slopes and curvatures were observed for X = F. Fourth, the correlations derived for halogen bonds with halide anions are very different from those established for some hydrogen bonds and halogen bonds formed by two neutral fragments (Table 1). This and the previous [116] works indicate that each type of bonding requires the application of its own relationship.
Fifth, the ρ b , V b , G b , and d Y···X estimators behave similarly. In contrast, the ∇ 2 ρ b , λ ||,b , and H b parameters have limited significance as E int estimators and they may be used only for some series or a certain range of the E int values.
The relationships recommended for practical use to estimate E int of the structures [(A) n Z-Y···I] − and [(A) n Z-I···X] − are given in Table 5. Table 5. Relationships recommended for the estimates of E int (in kcal/mol) for the structures [(A) n Z-Y···I] − and [(A) n Z-I···X] − at the M06-2X/ADZP-DKH level of theory (V b , G b , and H b in kcal/(mol•bohr 3 ), ρ b in e/Å 3 , λ ||,b in e/Å 5 , d Y···X in Å).

Series
Estimator Equation

Supplementary Materials:
The following are available online: Figure S1: Correlations between theoretical and experimental values of electron density; Figure S2: Correlations between −E int and −V b calculated at the M06-2X/ADZP-DKH and PBE0-D3BJ/ADZP-DKH levels of theory; Figure S3: Plot of the slope of the E int (V b ) dependence against the difference of Pauling's electronegativities of the atoms X and Y; Figures S4-S6: Plots of −E int against ∇ 2 ρb , λ ||,b , and H b ; Figure S7: Example of an adequate cluster for the direct estimates of E int between the fragments I-C 4 F 8 -I and Cl − in the X-ray structure ACIPOU; Figure S8: Four-molecule computational model of the structure ACIPIO; Table S1: Relationships between interaction energy and electron density-based estimators, Table S2: Calculated structures.

Data Availability Statement:
The data presented in this study are available in supplementary material and from the author upon request.

Conflicts of Interest:
The author declares no conflict of interest.