Tuning halogen bond donor ability: benchmarking substituent effects in fluoroiodobenzenes

Abstract

In high-resolution spectroscopy, electric nuclear quadrupole coupling can provide insights into a molecule's electronic structure. In a systematic study, we investigated the mesomeric and inductive effects exhibited by the fluorine substituents in various fluoroiodobenzenes. We characterised all 19 possible substitution patterns of fluorinated iodobenzenes using cavity resonator Fourier transform microwave jet spectroscopy, which included measuring the aforementioned coupling constants. By systematically implementing the extended Townes–Dailey model, we quantified the polarisation of the iodine atom to investigate its ability to act as a halogen bond donor, using the iodine atom as a probe. In this way it can be determined how halogen bond donor abilities are influenced by mesomeric and inductive effects. The results were compared to those obtained from an intrinsic basis bonding analysis (IBBA). Additionally, we used a large quantity of experimental data to rigorously benchmark numerous levels of theory in terms of their structure predictions and nuclear quadrupole coupling constants. To gauge the influence of zero-point effects on the rotational constants, vibrational perturbation theory of second order (VPT2) was also applied.

---

1 Introduction

Halogen atoms are commonly known to act as electron donor sites in inorganic salts or haloorganic molecules due to their electronegativity, which leads to an excess of electron density in comparison to the isolated atom. However, a halogen atom that is covalently bonded exhibits an anisotropy in its electron density. This effect is particularly pronounced in the heavier homologues due to their increased polarisability. On the side of the halogen atom opposite the covalent bond, an area of lower electron density and positive electrostatic potential forms, known as the σ-hole. This formation is facilitated by electron-withdrawing substituents. The area of electron depletion can interact with electron-rich binding sites, in a way that differs fundamentally from the non-covalent interactions commonly seen with halogen atoms acting as acceptors in hydrogen bonding. In a halogen bond, the halogen atom that is covalently bonded acts as an electrophile, forming an attractive, directed interaction with a nucleophile or nucleophilic area within the acceptor. This interaction can occur both inter- and intramolecularly, and is primarily driven by electrostatic forces. However, charge transfer and dispersion can also play a role.^1,2

The first halogen-bonded complex, I₂⋯NH₃, was reported over two centuries ago,^3,4 followed by dimers in which the lighter dihalogens act as donors. Other examples of known halogen bonds in the literature include halocarbons as donors, both aromatic and aliphatic, as well as interhalogens combined with various amines or other acceptors.¹ There are also cases where halide anions act as acceptors, forming halogen bonds with other halogen atoms. This has been observed, for example, in combination with iodopyridines.⁵ Systems that have been studied using rotational spectroscopy include H₂O⋯ICl,⁶ H₃N⋯ICl,⁷ H₂O⋯ICF₃,⁸ and H₃N⋯ICF₃.⁹

The concept of halogen bonding has numerous applications. Halogen bonds are widely used in crystal engineering in order to design and modify the desired properties of solids.¹⁰ In soft materials, they are found in polymers¹¹ and liquid crystals.¹² Since halogen atoms are present in many drug molecules, halogen bonding has many applications in the pharmaceutical field, e.g. in drug–target binding.¹³ Furthermore, halogen bonding is used in organic catalysis, where the directionality of the halogen bond, due to the position of the σ-hole, is exploited to design highly nucleophile selective ligands.¹⁴

As a group of structurally similar halogen bond donors for this study, the fluoroiodobenzenes, with all their substitution patterns resulting from the exchange of hydrogen atoms for fluorine atoms while retaining one iodine atom, were selected. As the title suggests, the motivation for this work was to systematically study the mesomeric and inductive substituent effects in relation to halogen bond donor ability. Microwave spectroscopy was selected as the experimental technique due to its ability to determine the molecular structure and resolve the hyperfine splitting caused by the quadrupolar iodine nucleus. The nuclear electric quadrupole coupling contains information of a nucleus's electronic environment, specifically the electric field gradient experienced by the nucleus. The resulting splitting allows for the application of the extended Townes–Dailey model^15,16 in order to gain more detailed information of the electronic structure of the aromatic molecules. Ionic characters and valence p-orbital occupations are derived to gauge the extent to which the electron-withdrawing effect of the aromatic substituent is influenced by certain substitution patterns. From the ionic characters and occupations we can then derive their potential as a halogen bond donors. A related analysis was conducted by Lv. et al.,¹⁷ in which the p_z-orbital anisotropy, e.g. the difference in valence p_z-orbital occupation compared to the average valence p-orbital occupation, was employed to measure the effects of full fluorination on the σ-hole strength in chloro- and bromobenzene. To aid the analysis and provide additional details that are not accessible through experiments, an intrinsic basis bonding analysis (IBBA) was carried out.

In addition, rigorous benchmarking against theory is performed using a large quantity of monomer data from 19 iodine-containing molecules that have not yet been studied by rotational spectroscopy. This benchmarking involves investigating the agreement between the optimised structures obtained experimentally and theoretically. Second-order vibrational perturbation theory (VPT2) was incorporated into the benchmark to determine the impact of zero-point effects on rotational constants. The selected group of molecules enables a systematic investigation of zero-point effects, given the increased rigidity of the systems with an increased degree of substitution. Furthermore, the electric field gradient calculations required for predicting the nuclear quadrupole coupling of the iodine nucleus are investigated by comparing the results with and without relativistic effect treatment against the experimentally determined nuclear quadrupole coupling constants. To enhance theoretical predictions based on experimental data, parametrisations in the form of semi-empirical scaling approaches are presented for both structure optimisation and electric field gradient calculation.

2 Methods

2.1 Nomenclature

Given the large number of molecules characterised, a nomenclature was defined to abbreviate the IUPAC names of all the fluoroiodobenzenes. Fig. 1 shows the Lewis structure of all fluoroiodobenzenes, which differ in the number of fluorine atoms and their position on the aromatic ring. This includes the non-substituted iodobenzene. Numbers indicate the positions of the fluorine atoms, counting from the ipso position 1, which is the carbon atom connected to the iodine atom. The only exception is iodobenzene, where 0 indicates the absence of fluorine atoms. For example, 2,3,5,6-tetrafuoroiodobenzene would be abbreviated with 2356. Iodobenzene (0) represents the only molecule that has already been studied by rotational spectroscopy.^18–21

Fig. 1 Lewis structures for all the different fluoroiodobenzenes and iodobenzene, using a shorthand notation. Numbers indicate the positions of the fluorine atoms, counting from the ipso position 1, which is the carbon atom connected to the iodine atom. The only exception is iodobenzene, where 0 indicates the absence of fluorine atoms. A training set (red) was created in order to scale the equilibrium rotational constants and to derive a semi-empirical scaling approach for the calculation of the nuclear quadrupole coupling constants, as described in Section 2.3.

2.2 Experimental methods

The chemicals in this study were bought and used without further purification, with the exception of 2,3,4,5-tetrafluoroiodobenzene and 2,3,4,6-tetrafluoroiodobenzene which were synthesised and purified in the lab of Prof. Manuel Alcarazo. The syntheses were conducted with a modified procedure of Perry et al.²² and the successful synthesis verified using literature ¹H and ¹⁹F NMR values^23,24 as well as HR-MS-EI spectrometry. Detailed information on the manufacturers, the substances' purity, the experimental and synthetic details can be found in Section S1 of the SI.

All experimental data were obtained from a pulsed-jet Balle–Flygare²⁵ type cavity resonator Fourier transform spectrometer, QCUMBER (Q-factor cavity utilising molecular beam electric resonator). The original setup from 1994²⁶ was updated on multiple occasions with the most recent version being described in ref. 27. Rotational spectra were recorded in the range between 7 and 16 GHz with the carrier gas varying between helium, neon, and argon at stagnation pressures between 1.1 and 1.6 bar. The resulting free induction decays (FIDs) were recorded for 51.2 µs at a maximum of 20 000 averages. Experimental details are summarised in Section 1.5 of the SI.

Grabow's FTMW++ software²⁸ was used for data extraction. For all microwave spectra measured in this work, Pickett's^29–31 SPFIT/SPCAT (32-bit version) programs were used paired with Kisiel's graphical interface AABS^32,33 to fit the Hamiltonian alongside the auxiliary programs found on Kisiel's Prospe³⁴ website. Specifically, PMIFST (version 25a.VII.2017) was used for structure calculations like isotope scaling/visualisation, QDIAG (version 12.II.2023) was used for the diagonalisation of the inertial quadrupole tensor with errors and PIFORM (version 31.XII.2016) was used for reformatting the fit output file. The Watson S^35–37 reduction of the Hamiltonian in the I^r representation was used throughout this work.

2.3 Theoretical methods

For a detailed analysis, numerous quantum chemical calculations on the monomers were performed. Density functional theory (DFT) based calculations were employed using the ORCA 6.0.1^38–40 program package. MOLPRO 2024.1,^41–43 and Gaussian 16 (Rev. C.01)⁴⁴ were used for wavefunction (WF) theory and DFT based methods. The specific use cases and the respective levels of theory for the different program packages are discussed below. A summary alongside exemplary input files can be found in Section S5 of the SI. Where applicable, the frozen core approximation was used. All the molecules analysed were treated equally.

Geometry optimisations of the fluorinated iodobenzenes and iodobenzene were carried out with Gaussian in order to use its implementation of vibrational perturbation theory of second order (VPT2)^45,46 to account for zero-point effects in the calculation of the rotational constants and to enable benchmarking against different levels of theory. Note that the geometry optimisation and the VPT2 calculation were performed separately. The choice of functionals was inspired by ref. 47. The def2-TZVP⁴⁸ basis set was used to limit computational cost. To compare with the ubiquitously used B3LYP functional,^49–52 a second hybrid functional, PBE0,^53,54 was used. To represent the range-separated hybrid functionals,⁵⁵ where the added exchange is no longer constant, CAM-B3LYP⁵⁶ and LC-ωPBE⁵⁷ were chosen. The empirically fitted meta-hybrid M06-2x⁵⁸ functional, as well as two different double-hybrid functionals, B2PLYP^59,60 and DSD-PBEP86,^61,62 were also added to the benchmark. Double-hybrid functionals use parts of the correlation energy, obtained from second-order Møller–Plesset perturbation theory (MP2).^63–68 DSD-PBEP86 uses spin-component scaled MP2 as an improvement upon B2PLYP.⁶⁹ All DFT methods were dispersion corrected using Grimme's D3(BJ) scheme.^70,71 Structures were also optimised with MP2, in order to represent the post Hartree–Fock methods.

Additionally, geometry optimisations and frequency calculations, within the double harmonic approximation, were carried out with ORCA at the B3LYP level of theory. Unlike in the Gaussian calculations, the minimally augmented version of the def2-TZVP basis set (i.e. ma-def2-TZVP) of Truhlar and co-workers⁷² was employed as well as three body terms^73,74 being added to the dispersion correction (D3(BJ,abc)).

Due to the similarity of the fluoroiodobenzenes structures, the training set highlighted in red in Fig. 1 was used to obtain the scaling parameters p_A, p_B, and p_C connecting the calculated equilibrium rotational constants A_theo,e, B_theo,e and C_theo,e to the experimentally determined rotational constants. This was done to provide a point of reference for the performance of the VPT2 calculations. The training set included one commercially available representative molecule for each degree of substitution with the largest dipole moment. For this, the results of the geometry optimisation at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory were employed. To obtain the scaling parameters, the ratios of the experimental and theoretical values were calculated, A_exp/A_theo,e, B_exp/B_theo,e and C_exp/C_theo,e. The averages of these ratios for the six molecules then yield the scaling parameters p_A, p_B, and p_C.

Experimental nuclear electric quadrupole coupling constant (NQCC) tensors were benchmarked against a variety of different levels of theory for a large set of experimental results. To achieve this, EFG calculations were performed in ORCA using the numerous different DFT functionals listed below. An attempt was made to parametrise an effective nuclear quadrupole moment, Q_eff, in order to replace the literature value Q = −688.22 mb.⁷⁵ This procedure extends the work of Bailey,⁷⁶ who found a linear correlation between the measured and predicted NQCC values for iodine at the B1LYP/6-311G(df,p) level of theory. For iodine, an f-polarisation function with a coefficient of 0.4 was added, by Bailey, to the 6-311G(d,p)⁷⁷ basis set, resulting in a linear NQCC correlation for a training set of 16 molecules. Only the diagonal components of the NQCC tensors were considered within the inertial principal axis system. As an inexpensive entry level-method for benchmarking, the GGA functional PBE⁷⁸ was paired with the def2-SVP^48,79 basis set. For all other functionals, the x2c-TZVPPall⁸⁰ basis set was used. The B3LYP, B2PLYP, BP86,^49,81,82 B97,⁸³ and B97M-V⁸⁴ functionals completed the set. To take relativistic effects into account, the second order Douglas–Kroll–Hess transformation (DKH2)^85–87 was employed, in combination with the finite-nucleus model⁸⁸ and picture change effects.⁸⁹ For comparison, the exact two-component theory (X2C) implementation^90,91 was used as well, combined with the aforementioned finite-nucleus model and picture change effects. The choice of functionals was inspired by ref. 92, in which chlorine NQCCs were benchmarked. Sun et al.⁹³ recently adopted Bailey's scaling approach, where relativistic effects were not explicitly treated, revealing that B3LYP can be parametrised well as described above.

An intrinsic basis bonding analysis (IBBA)⁹⁴ was performed in MOLPRO for all monomers, employing the structures optimised in ORCA at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level. Herein, the exact wavefunction from a self-consistent field calculation is represented by atomic valence and core orbitals, which are polarised by the molecular environment. From this, the occupations of atomic orbitals, partial charges, and the contributions of these orbitals to the partial charges of atoms can be deduced. This analysis can be compared to the extended Townes–Dailey model presented in the following subsection. B3LYP-D3(BJ,abc) was paired with the cc-pVTZ-DK^95,96 basis set along with the so-called intrinsic minimal basis (IMB) MINAO/MINAO-PP, from which the molecular Lewis structure emerges.⁹⁴ Relativistic effects were again treated with DKH2, employing the finite-nucleus model. Picture-change effects cannot be individually defined in MOLPRO. The chosen basis set has been specifically contracted for use with DKH2, enabling analysis of the importance of relativistic effects in terms of iodine electron occupations.

2.4 Extended Townes–Dailey model

2.4.1 Theory. In 1949, Townes and Dailey¹⁵ formulated a model to obtain information about the electronic structure of a molecule from its nuclear quadrupole coupling constants. This model was later applied to estimate the ionic characters of diatomics.⁹⁷ In 2011, Novick¹⁶ proposed an extension to the Townes–Dailey model, expanding it to the three dimensional space for polyatomic molecules while also taking hybridisation into account. In general, given its assumptions and restraints discussed below, the model aims to strike a balance between achieving a pure chemically intuitive understanding, and determining the true nature of a bond, which would require the precise measurement of valence electron occupations. In the following, a brief description of the extended Townes–Dailey model will be given, as many of the following results are based on or derived from Novick's extension of the original model. Unless stated otherwise, all information was taken from ref. 16.

In both, the original and the extended cases, the electronic structure derived from NQCCs only considers valence p-electrons. This strict cut-off may no longer be a good approximation for heavy nuclei, such as iodine, as relativistic effects become more significant. The nuclear quadrupole coupling requires a non-zero electric field gradient, which results from a break in spherical symmetry. This is not the case for s-orbitals, which are spherically symmetric and therefore have no EFG. According to Novick, 3d-orbitals contribute 47 times less than the 2p-orbitals. The relationship between the NQCCs and the p-orbital populations n_x, n_y and n_z is shown in eqn (1)–(3). For the systems studied here, the NQCC resulting from the one unpaired 5p electron in an isolated iodine atom is denoted by χ₀ (eQq₅₁₀). In the case of ¹²⁷I, this value is equal to +2292.712 MHz.⁹⁸χ_xx, χ_yy and χ_zz are the NQCCs, defined in the nuclear principal axis system, with |χ_xx| ≤ |χ_yy| ≤ |χ_zz|.

Experimentally, the NQCCs in the inertial principal axis system are obtained. Depending on the molecular geometry, a rotation to the nuclear axis system, and therefore knowledge of the whole tensor, including the off-diagonal components is required.⁹⁹ This kind of rotation is qualitatively illustrated in Fig. 2 for 2-fluoroiodobenzene. In this case, the b-axis in the inertial principal axis system now becomes the y-axis in the nuclear case, while the a-axis becomes the z-axis. The previously defined n_x, n_y and n_z correspond to the p-orbitals in the nuclear framework. The Townes–Dailey model assumes the z-axis in this example to be oriented exactly along the C–I bond, which is not always true, where the NQCC has non-zero off-diagonal components.

Fig. 2 Illustration of the structure of 2-fluoroiodobenzene in the ab-projection at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. The inertial principal axes a and b are shown in black. Qualitative nuclear principal axes y and z are blue coloured, showing the rotation from inertial to the nuclear principal axis system. The c and x axes positioned at the crossing points of the two remaining axes for the respective axes system are pointing out of the plane of projection and are not labelled.

Since the NQCC tensor is traceless, all the matrices following eqn (1)–(3) must also be traceless. This leaves only two independent parameters, resulting in the equation being underdetermined when solving for the three p-orbital populations, leading to an infinite number of sets of solutions. In the original formulation of the model, n_x and n_y are set to 2, but in the extended model, a restraint is added so that the populations are within 0 and 2, with the populations of the non-bonding orbitals being as close as possible to 2 in order to identify chemically sensible solutions. This allows for more flexibility in the extended version.

2.4.2 Systematic implementation. Novick's¹⁶ extended model Townes–Dailey model has been implemented in a systematic way. First, the Townes–Dailey model shown in eqn (1)–(3) is defined with added penalty functions to fulfil the first restraint defined by Novick, whereby the non-bonding orbital populations should stay as close to two as possible. To account for uncertainties in the experimental results for the NQCCs, all possible combinations are considered for each tensor component: without uncertainty and with uncertainty added. This results in nine combinations in total. To take into account different starting guesses, a three-dimensional linear space has been created, spanning four different starting guesses, 0, 2/3, 4/3 and 2, resulting in 64 total starting combinations. The Python package Scipy¹⁰⁰ was used to solve eqn (1)–(3) with the added penalty functions. Bounds were set for the p-orbital occupations to be between zero and two, thus fulfilling Novick's second restraint for chemically reasonable results. For each solution, the resulting objective value is calculated resembling how close the results from n_x, n_y and n_z match the NQCCs and how well the restraints are fulfilled. The objective value is defined in eqn (4). Additionally, the sum of the p-orbital occupations n_tot was calculated. As a preliminary uncertainty, the absolute difference in n_tot between the lowest and highest objective value was calculated.After repeating the procedure for all possible combinations and starting guesses, a total of 576 least-square problems per molecule, the final best result was found to be the one with the lowest objective value, o, of all the results. To calculate the final uncertainty, the maximum difference in n_tot between this best result and any other combination, each including its preliminary error bar, was determined. This procedure was repeated for each molecule analysed. The ionic character I_TD of the I–C bond was calculated according to eqn (5) as the difference from the sum of the valence p-orbital populations to a population of five electrons for the isolated iodine atom ([Kr]5s²4d¹⁰5p⁵). As all three valence p-orbitals are considered, this definition of the ionic characters encompasses both σ- and π-characters. A negative I_TD indicates a gain in valence electrons, whereas a positive value corresponds to an electron depletion. The uncertainties of I_TD depend linearly on n_tot and were calculated directly using its uncertainty.

I_TD = (5 − n_tot) × 100% (5)

3 Results and discussion

Due to the large number of systems studied, the experimental results regarding rotational constants, distortion constants, and nuclear quadrupole coupling constants are not discussed in detail, but are summarised in the SI (see Section 1.3). For the sake of completion, results for the non-fluorinated iodobenzene (i.e. 0) have been taken from ref. 21 and added to our own data set. As the original fit uses the Watson A reduction in the I^r representation, the data provided in the supplementary information was refitted employing the Watson S reduction, to ensure consistency with the other fits in this work.

3.1 Parametrisation of rotational constants

Following the procedure described in the previous section, the scaling parameters p_A, p_B, and p_C were obtained from the training set, highlighted in red in Fig. 1, as shown in Fig. 3. As stated above, the training set was constructed by selecting, from commercially available sources, one representative molecule for each total number of substituted hydrogen atoms with the largest dipole moment to ensure good measurability. This was done to potentially improve the B3LYP-D3(BJ,abc)/ma-def2-TZVP calculations for the remaining molecules, given the structural similarities, and therefore reduce measurement times and simplify assignments. An important example of this is para-fluoroiodobenzene (4) with its small predicted dipole moment of 0.06 D. This simple approach did not consider the experimental uncertainties.

Fig. 3 Ratios of the experimental results for the rotational constants A, B, and C and theoretical results. The theoretical results were obtained from geometry optimisation at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. The ratios A_exp/A_theo,e, B_exp/B_theo,e and C_exp/C_theo,e are plotted for the molecules 0, 2, 26, 234, 2356, and 23456, which form the training set, highlighted in Fig. 1. The black dashed line indicates a ratio of exactly 1. Scaling factors p_A, p_B and p_C were calculated as the average of the ratios for a given rotational constant. Experimental uncertainties were not considered. Results for 0 were taken from ref. 21.

The rotational constants B and C computed at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level follow a clear trend, with the experimental results consistently surpassing the theoretical ones. For A, however, theory overestimates the experimental rotational constants. The ratios A_exp/A_theo,e show larger fluctuations, especially for 0, with a ratio under 0.991, and slowly approach the dashed black line, indicating perfect agreement between theory and experiment. It should be noted that the rotational constants from theory (A_theo,e) do not include any zero-point vibrational effects, since they are only equilibrium structures. This simple scaling approach indicates that scaling with a constant factor may be sufficient to describe these zero-point effects in the case of B and C. For A, however, zero-point effects are particularly significant for 0 and vary considerably between the molecules. In the inertial principal axis system of 0, the a-axis coincides with the C–I bond. Consequently, the A rotational constant is greatly influenced by zero-point vibrations of the four C–H bonds in ortho and meta positions relative to the iodine atom. As the number of fluorine atoms increases, the molecules become more rigid. In the diatomic case, the bond dissociation energies are 338.72(11) kJ mol⁻¹ for the C–H bond and 513.8(10.0) kJ mol⁻¹ for the C–F bond.¹⁰¹ This results in zero-point vibrations affecting A less, and the equilibrium rotational constants better resembling the experimental results. In general, zero-point effects are expected to reduce the rotational constants as the bonds are elongated, resulting in larger moments of inertia. Therefore, a ratio between experimental and theoretical equilibrium value smaller than one is expected, a result that is only observed for A. The fact that the ratios for A are close to one for the heavy systems, and that constant ratios over one are obtained for B and C highlights that the calculations are close to the experimental results for the wrong reasons as the inclusion of zero-point effects would push these ratios to even larger values.

3.2 Benchmarking rotational constants

In order to evaluate the performance of different levels of theory, numerous VPT2 calculations were performed, as described in Section 2.3. The def2-TZVP basis set was used in all cases, alongside D3(BJ) for the DFT-based methods. The zero-point corrected rotational constants A₀, B₀ and C₀ obtained from the VPT2 calculations are compared to their experimental counterparts. The rotational constants obtained from the scaled B3LYP-D3(BJ,abc)/ma-def2-TZVP described in Section 3.1 were added to the comparison and are indicated by the “scaled” addition to the functional. Even though they are different quantities, the equilibrium rotational constants A_e, B_e and C_e are also compared to the experiment as Fig. 3 already indicated a possible match in some cases. Here, the unscaled B3LYP-D3(BJ,abc)/ma-def2-TZVP was added to the comparison. The relative deviation (RD) was chosen as a performance metric. The RD values were calculated according to eqn (6). Here, x_i represents an observable, with the experimental value being denoted as x_exp,i for the molecular system i. Due to their small magnitudes, experimental uncertainties were not considered.

To enable a general comparison between the different levels of theory, both the RDs for the equilibrium rotational constants as well as their zero-point corrected counterparts are calculated. This is illustrated in Fig. 4 in the form of boxplots.

Fig. 4 Boxplots showing the three relative deviations RD(A_e), RD(B_e), RD(C_e), RD(A₀), RD(B₀), and RD(C₀) of the results of quantum chemical calculations from the experimental results for the rotational constants. The RD of the individual rotational constants is defined according to eqn (6). The zero-point corrected rotational constants A₀, B₀ and C₀ (right) were obtained from VPT2 calculations employing the functionals on the abscissa and utilising the def2-TZVP basis set. The equilibrium rotational structures A_e, B_e and C_e (left) were calculated at the same levels of theory. Details of the functionals used are given in Section 2.3. One exception to this is the method labelled with “scaled B3LYP-D3(BJ,abc)”. In this case, an experimental training set was used to scale the calculated rotational constants obtained from a pure geometry optimisation at the B3LYP-D3(BJ,abc)/ma-def2-TZVP of theory without a VPT2 calculation. This procedure is described in Section 3.1. This level of theory has a self-referencing problem due to the inclusion of its training set into the test set, as discussed in the text. The orange line in a box represents the median of a group of values. The position of the averages is indicated by purple symbols. Boxes extend from the first quartile to the third quartile. Maximum whisker lengths extend to the first/third quartile minus/plus one and a half times the height of the box. Fliers are indicated with dots. Exact values can be found in the SI (see Section 3).

For most of the compared functionals, the zero-point corrected rotational constants more closely match the experimental results than the equilibrium counterparts do. However, this is not the case for B3LYP and B2PLYP and the RD(A_e) of MP2, where better agreement is observed without the VPT2 correction. For B3LYP in particular, this was suggested in Section 3.1, where the theoretical equilibrium results were already underestimating the experimental rotational constants. As discussed, applying zero-point effects pushes the theoretical results towards smaller values, thereby increasing the deviation from the experimental results. Despite the results for A_e from MP2, B2PLYP without zero-point corrections outperforms all other functionals, even the scaled approach, and is the best-performing method overall. For A_e, B3LYP without zero-point effects, on the other hand, offers comparable performance to PBE0, B2PLYP, DSD-PBEP86, and MP2 without zero point effects. This highlights that, even though equilibrium rotational constants should not offer better agreement with experimentally obtained rotational constants by definition, as they neglect zero-point effects, comparable or even better agreement can be obtained by neglecting these effects for selected functionals, which significantly limits computational costs. However, despite the good agreement observed, the non-zero-point-corrected B3LYP and B2PLP functionals should not be considered “good”, given that they are correct for the wrong reasons, as outlined above. By definition equilibrium rotational constants should not agree better than their zero-point corrected counterparts.

Although taking zero-point effects into account does not necessarily improve the agreement between experiment and theory, it does decrease the spread of the deviations, as can be seen from the narrower box sizes, particularly for PBE0, DSD-PBEP86 and MP2. This limited spread highlights the importance of anharmonic effects in the systems studied.

Some general trends in the functional families can be identified when going from negative to positive RDs in the zero-point corrected rotational constants. The RDs observed for B2PLYP, which is based on the exchange and correlation terms of B3LYP with the addition of MP2 correlation, fall between those of its two building blocks. More positive RDs are observed for the range-separated hybrid functional CAM-B3LYP than for the non-range-separated B3LYP. The same is true of the pair formed by PBE0 and its range-separated counterpart, LC-ωPBE.

The worst-performing functional for all three rotational constants is LC-ωPBE, for which the global maximum RD(A₀) of 2% is observed for the molecules 0 and 4. Similar performance is observed for PBE0, CAM-B3LYP and M06-2X, all of which overestimate the experimental values. Indicating an underestimation of the experimental rotational constants, B3LYP, B2PLYP and DSD-PBEP86 all exhibit negative RDs. The overall second-best performance with applied zero-point corrections is observed for DSD-PBEP86, surpassed only by the scaled B3LYP as the best-performing method in this group. Given that this variant uses actual experimental results, and that the molecules are expected to behave similarly due to their structure, DSD-PBEP86 offers good performance, as it is nearly able to match the RD(ABC) distribution of the scaled method. In this regard, it should be noted that the current benchmarking procedure is subject to self-referencing, given that the training set used for the scaled approach is also included in the test set. Without removing the training set from the test set, it is impossible to determine the method with the best absolute performance. However, comparisons between the ab initio approaches remain valid.

Nevertheless, this comparison highlights that for similar molecular structures, scaling using previously obtained experimental results is not only beneficial for the results, but can also outperform methods requiring greater computational costs. Therefore, not only computational time is saved, but spectral assignments are also simplified. For reference, the DSD-PBEP86 calculation for 2345 took about 20 times longer than the B3LYP VPT2 calculation. In comparison to the B3LYP-D3(BJ,abc)/ma-def2-TZVP calculation, it took ca. 900 times longer. Furthermore, this result indicates similarities in the zero-point effects on the rotational constants among the fluoroiodobenzenes, as a scalar multiplication is sufficient in order to compensate for the absence of zero-point effects in the geometry optimisations, eliminating the need for additional VPT2 calculations, as is the case for the scaled B3LYP-D3(BJ,abc)/ma-def2-TZVP. The exact RDs for each of the rotational constants are shown in the SI (see Table S16).

It should be noted that this simple scaling approach does not benefit from a good agreement between theory and experiment prior to scaling, but rather from a small spread in deviations. For this reason, PBE0-D3(BJ)/def2-TZVP is a better scaling alternative than B3LYP-D3(BJ, abc)/ma-def2-TZVP, as it exhibits less spread in A_e while retaining computational demand. This highlights the importance of benchmarking even for semi-empirical scaling.

3.2.1 Planar moments. The experimentally determined planar moments P_aa, P_bb, and P_cc are summarised in Table S8. Given the molecular geometries, small values of P_cc are obtained, with a maximum deviation from 0 amu Ų of about −0.082 amu Ų for 25. The P_aa and P_bb values for species with C_2v symmetry follow the expected trends. For example, P_bb increases by approximately 1.779 amu Ų when going from 35 to 26. This illustrates the repulsion between the two ortho-positioned fluorine atoms and the iodine atom in 26, as opposed to the attraction between the two hydrogen atoms and the iodine atom in 35. In the latter case, P_bb is reduced due to a possible decrease in the H–C–I angle compared to the F–C–I angle in 26.

If perfect planarity of the systems studied was assumed, P_aa and P_bb would be sufficient for benchmarking theory against experiment. However, this can only be exactly true for theoretical equilibrium structures, as it neglects the subtle deviations caused by zero-point effects that affect both experimental and zero-point-corrected results. Plots showing the relative deviations for both P_aa,e and P_bb,e to their theoretical equilibrium counterparts can be found in the SI (see Fig. S3). No comparison is drawn for P_cc,e, since deviations from 0 amu Ų may be caused by numerical noise in the calculations. Additionally, the same plot is shown for the relative deviations to the zero-point-corrected values (see Fig. S4). Here, P_cc,0 was also added. In general, the performance of the functionals is retained when the number of parameters under investigation is reduced. However, the RDs (P_cc,0) show that significant relative deviations are obtained, exceeding −600% for 23456 when calculated with the M06-2x functional. When comparing experimental and zero-point-corrected values, the ab initio approaches are able to correctly predict almost all signs. Disagreements only occur for 4, where MP2 is the only method that correctly predicts a positive sign. All other functionals fail in this regard, which suggests problems with the vibrational contributions. This positive sign has also been observed for the lighter homologue, chloro-4-fluorobenzene.^102,103

3.3 Benchmarking nuclear quadrupole coupling constants

In order to benchmark experiment against numerous levels of theory, the new data on the fluoroiodobenzenes and the existing data on iodobenzene was used to expand Bailey's original training set. Computational details are provided in Section 2.3. To analyse the influence of relativistic effects, calculations were performed using DKH2 and X2C, as well as without any relativistic treatment. It should be noted that a possible source of error is that zero-point effects were not considered in the single-point calculations of the electric field gradients. These subtle effects are described by VPT2 calculations but given the large number of systems and functionals, they were not analysed as relativistic effects are expected to be significantly more impactful.¹⁰⁴ Experimental results for Bailey's training set,⁷⁶ consisting of HI,¹⁰⁵ CH₃I,¹⁰⁶ CH₂I₂,¹⁰⁷ CH₃CH₂I,¹⁰⁸ IN₃,¹⁰⁹ IF,¹¹⁰ ICN,¹¹¹ HOI,¹¹² CH₃CHICH₃,¹¹³ CH₂=CHI,¹¹⁴ g-CH₃–O–CH₂I,¹¹⁵ I₂,¹¹⁶ HCCI,¹¹⁷ OIO,¹¹⁸ CH₂ICl (³⁵Cl species)¹¹⁹ and (CH₃)₃SiI¹²⁰ were taken from the literature. Analogous to eqn (6), the RDs are calculated for the individual diagonal tensor components in the principal inertial axis system, denoted with RD(χ_aa), RD(χ_bb) and RD(χ_cc).

In addition, Bailey's original Q_eff = −739.8(1.2) mb,⁷⁶ resulting from the linear regression of the EFG calculations on the experimental structures, was used for comparison. The original version is labelled “Bailey”.

As this training set is included in the test set, self-referencing occurs. The rest of the test set comprises fluoroiodobenzenes and iodobenzene. Due to their structural similarities, the test set is weighted. Therefore, we have self-referencing, which benefits Bailey's parametrisation, whereas a weighted test set that puts Bailey's parametrisation at a disadvantage, as these structures are not part of the training set. These effects only impact the performance of Bailey's approach, not the ab initio approaches. Therefore, a direct comparison between them is not possible.

Results for all RDs are shown in Fig. 5, again in the form of boxplots. The PBE/def2-SVP results are omitted because the RDs show large, constant deviations of around ±100%, which reduces the level of detail visible in the remaining boxes. One possible cause of this is the small basis set with its effective core potential, which is unable to adequately describe the iodine atom. The largest RDs and therefore the worst results (generally over 20% deviation) were observed for the methods without relativistic corrections. To allow for a better comparison, the left column of Fig. 5 only shows the parametrised method, as well as those with enabled relativistic corrections. However, even though relativistic corrections are also absent for Bailey's method, the results are better since they are semi-experimental. To investigate whether the large RDs for the methods without relativistic corrections are caused by the basis set, the EFGs for B3LYP-D3(BJ,abc) with and without relativistic corrections were recalculated using the decontracted basis set. This was done to minimise the uncertainties introduced by the basis set, which was specifically contracted for the use with relativistic corrections. The results are shown in Fig. S5 of the SI. Similar results are obtained for both relativistic corrections. Without an applied correction, the results are improved by decontracting the basis set, but still perform significantly worse than the corrected ones. This highlights the importance of relativistic corrections for iodine.

Fig. 5 Boxplots of the relative deviations RD, defined according to eqn (6), of the three diagonal components of the nuclear quadrupole coupling tensor in the inertial principal axis system, χ_aa, χ_bb and χ_cc. The plots in the column on the right use parametrised results, obtained from linear regressions using all three diagonal components of all molecules. Numerous different functionals have been employed for the theoretical calculations of the electric field gradients, paired with the x2c-TZVPPall basis set, unless otherwise stated. “Bailey” denotes a special level of theory, explained in Section 2.3. This level of theory has a self-referencing problem due to the inclusion of its training set into the test set, as discussed in the text. Exact values are summarised in Table S19 of the SI.

In many cases, a relatively large spread of the fliers can be observed, particularly for B97-D3(BJ,abc). It should be noted that this benchmarking procedure is not solely a measure of the EFG calculation. In contrast to Bailey's original approach, in which EFG calculations were performed on experimentally determined structures, only single-point calculations on top of the optimised structures are performed. If a structure shows a large deviation from the experiment, this error is propagated to the EFG calculation. Following this argument, the RDs for “Bailey” highlights cases where the experimental structure and the result of the geometry optimisation at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory may disagree. Overall, the two relativistic corrections show a very similar performance. In the case of B97-D3(BJ,abc), B97M-D3(BJ,abc) and BP86-D3(BJ,abc), X2C is the better choice. For B3LYP and B2PLYP, DKH2 shows better agreement. The best overall performance is obtained for B3LYP-D3(BJ,abc)-DKH2.

To reveal possible projection errors, all NQCC tensors were rotated to the nuclear axis system, i.e. the NQCC tensor was diagonalised. Experimental rotation angles α_az between the inertial principal axis a and the nuclear principal axis z can be found in the SI (see Table S9). This approach reduces the influence of errors introduced by optimised structures that deviate from their experimental counterparts. The results are shown in Fig. S6 of the SI. Generally, the agreement between experiment and theory is better in this comparison as the spread of flyers is significantly reduced. Interestingly, B97M is the only functional that can consistently provide an adequate description of all the different molecules in the test set, including the open-shell OIO (2S + 1 = 2). For the latter, all other levels of theory show significant discrepancies between theory and experiment, especially for χ_xx for which an otherwise minimal RD of −28.43% is observed.

In the final step, all the different methods were parametrised using the experimental data from all 36 molecules (Bailey's training set + iodobenzene and iodofluorobenzenes), in order to find the method best suited to a new parametrisation specific to molecular structures optimised at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. This approach does not measure how well a method performs, but rather how consistently obtained results differ from the experimental ones. All results are visualised in the right column of Fig. 5. Consistent with previous findings, poor performance is obtained for OIO and CH₂I₂ due to error propagation following the geometry optimisation. Relativistic corrections again improve results for every functional. Decontracted results are not considered in this regression analysis. In contrast to the large RDs without the parametrisation, PBE/def2-SVP now shows competitive results. Although the deviations were large, they were also similar across the different molecules, which is beneficial for the parametrisation. Overall, the results obtained with B3LYP-D3(BJ,abc)-DKH2 are the best, surpassing Bailey's proposed level of theory.

The deviation between the effective nuclear quadrupole moment Q_eff, as determined by the various linear regressions, and the literature value for iodine's quadrupole moment, Q = −688.22 mb,⁷⁵ provides additional insight into the parametrisations. The quotients Q/Q_eff reveal by how much and in which direction the calculated values needed to be corrected in order to agree with the experiment. The results are shown in Fig. S7 of the SI. As expected from the poor initial PBE/def2-SVP results, a very small quotient of 0.0177(3) is obtained, indicating that the initial absolute values are consistently too small. The opposite effect is observed for the methods without relativistic corrections, albeit to a lesser extent. At Q/Q_eff = 0.917(12), Bailey's originally proposed level of theory stands out, as it requires the least correction of Q out of all the methods without relativistic corrections. Here, calculations were performed at the B1LYP/6-311G(df,p) level of theory, for which an f-polarisation function with a coefficient of 0.4 was added to the 6-311G(d,p) basis set for iodine. Of all the methods employing relativistic corrections, B3LYP-D3(BJ,abc)-DKH2 performs best and gives the closest value to one at Q/Q_eff = 1.0035(93).

3.4 Extended Townes–Dailey model

This subsection presents the results following the implementation of the extended Townes–Dailey model, starting with the ionic character I_TD. Following the definition of I_TD in Eq. 5, a decreased ionic character indicates more positive polarisation of the iodine atom. Given the geometries of the studied systems, the strength of the σ-hole and therefore the increased halogen bond donor ability should be directly proportional to a decreased valence p_z-orbital occupation n_z.

Fig. 6 shows the I_TD of all structures, sorted from smallest to largest. The exact values can be found in Table S11 of the SI. Three distinct sections can be identified, separated by two significant jumps in I_TD. From left to right, the first group consists of the molecules with two fluorine atoms in the ortho position, followed by the largest group with one ortho fluorine atom, and at last the group without fluorine atoms in the ortho position, which has the most positive I_TD. A general trend can be observed among all groups. Not only is the number of fluorine atoms important, but more so is their position. In all cases, the ortho fluorines have by far the largest contribution on the total electron-withdrawing effect the iodine atom experiences. The effect of the meta and para fluorines is significantly weaker, with meta having the stronger effect. Therefore, proximity to the iodine atom is more important than the mesomeric effects, in which the para fluorine's effect would be stronger. Two important edge cases can be identified, involving two pairs of molecules with the exact same number of ortho-, meta- and para-positioned fluorine atoms. The first such pair is formed by 23 and 25, and the second one by 234 and 245. In both cases, the variant with the fluorine atoms distributed more evenly/symmetrically around the aromatic system, namely 25 and 245, has a more negative I_TD. In these cases, the interaction between the fluorine atoms is minimised, as their distance is maximised, which may strengthen their effect on the iodine atom. In the case of 23, the fluorine atoms have an ortho interaction, whereas in 25, this interaction is switched to a para interaction. The same can be observed for 234 and 245.

Fig. 6 The results for the ionic character I_TD of the fluoroiodobenzenes and iodobenzene, obtained from the implementation of the extended Townes–Dailey model. The nomenclature for the molecules on the abscissa follows the nomenclature presented in Fig. 1. Error bars are always plotted but are not visible in all cases due to their small relative magnitude. Exact values are summarised in Table S11 of the SI.

Further insights into the valence orbital interaction of the iodine atom with its bonding partner can be obtained by examining the occupations n_x, n_y and n_z of the p-orbitals along their respective axes. For a definition of the axis system, see Fig. 2. The model's implementation calculates the uncertainties by maximising the difference between the total occupation, n_tot, and the best result with respect to the lowest objective value to any other solution of the model. However, given the correlation between n_x, n_y and n_z, this approach is not sensible as the uncertainties would exceed the results obtained, since the model is still underdetermined and has an infinite number of sets of solutions. Therefore, the uncertainties of the p-orbital occupations will be neglected in this systematic implementation. Fig. 7 shows the results as a bar plot, and the exact values are again summarised in Table S11 of the SI. The molecules are sorted in order of increasing I_TD, as shown in Fig. 6. Here, n_x represents the occupation of the valence p-orbital perpendicular to the aromatic ring, for which small changes of up to 0.042 electrons can be observed for different substitutions. Overall, the occupation of the p_x-orbital increases from left to right. As shown in Fig. S1 of the SI, a linear correlation is obtained between n_x and I_TD with an R² value of 0.939. This strongly suggests that the iodine atom interacts with the π-system of the aromatic ring and donates electron density via the interaction with its perpendicular p-orbital. Larger changes in occupation of up to 0.158 electrons are obtained for the p_z-orbital, which is oriented in the direction of the C–I bond. Here, a better linear correlation between n_z and I_TD is obtained, with an R² of 0.995. The iodine atom donates more and more electron density to the aromatic substituent via this σ-bond, as the withdrawing effects increase. No change is observed for the p_y-orbital, which lies in the plane defined by the ring. In all cases, this orbital is filled with two electrons. In summary, the iodine atom interacts with the π-system via its perpendicular p_x-orbital, whereby electron density is donated to the aromatic substituent. This interaction decreases from left to right. However, a larger change in interaction with the aromatic substituent occurs via its σ-bond to the carbon atom, where the iodine atom's acceptor interaction increases from left to right. Both of these changes lead to an increase in the ionic character.

Fig. 7 The results of the valence p-orbital occupations n_x, n_y and n_z of the fluoroiodobenzenes and iodobenzene, obtained from the implementation of the extended Townes–Dailey model, which is described in section 2.4.2. The nomenclature for the molecules on the abscissa follows the nomenclature presented in Fig. 1. Uncertainties are neglected. Exact values are summarised in Table S11 of the SI.

The implementation described above relies on two penalty terms, resulting in n_x and n_y having values close to two. Adding the same penalty term to n_z results in exactly the same I_TD and occupations, with only minor deviations in the uncertainties of I_TD. Therefore, the described approach does not influence the result by biasing certain orbitals. Enough information on the orbital occupations is contained within the measured NQCC tensor. However, the total removal of penalty terms does not produce a chemically reasonable result, as the optimisation aims to minimise the objective value o, resulting in n_z = 0 for the majority of the molecules. Since both the implementation with penalty terms for all, and the implementation with penalty terms only for two orbitals results in n_y = 2 for all structures, this observation can be used to explicitly determine the unique solution to the extended Townes–Dailey model. Setting n_y = 2 leaves only two parameters, resulting in the system no longer being underdetermined. The explicit expressions for n_x, n_y and n_z are given in eqn (7). Uncertainties are calculated via Gaussian error propagation. The exact results are summarised in Table S12 of the SI, showing identical results for I_TD, compared to the minimisation implementation. In some instances, the calculated occupations differ slightly between the two implementations, especially for n_x. A maximum deviation of 0.12 between the two is obtained in n_x for 236. These fluctuations may be due to rounding errors, since the original n_y values are not exactly 2. Overall, the minimisation approach combined with penalty functions is effective and can exactly reproduce the I_TD obtained from solving the model explicitly with added restraints. Therefore, it is expected that the results for other systems, for which no such restraints are warranted, will be accurate within the limitations of the extended Townes–Dailey model. This highlights that our implementation is more broadly applicable.

In order to compare the obtained results from the extended Townes–Dailey model, an intrinsic basis bonding analysis (IBBA) was conducted. As outlined in the computational details, the MINAO-PP basis set was selected as the default intrinsic minimal basis (IMB) for the iodine atom. The molecular Lewis structure emerges from this IMB. According to the MOLPRO manual, this basis set was constructed by stripping all the polarisation functions from the cc-pVTZ-PP basis set. Based on the electron configuration of the isolated iodine atom, the 5s-, 5p- and 4d-orbitals are expected to be the highest energy representatives of their respective orbital types. However, manual inspection of the IBBA outputs revealed that virtual space was applied up to the 7s-, 7p- and 6d-orbitals, with the 7s- and 7p-orbitals still showing meaningful occupations. Most likely the cc-pVTZ-DK basis set without f-orbital contributions is used for iodine instead of cc-pVTZ-PP. To obtain the valence orbital occupations within the framework of an isolated iodine atom, the number of electrons of a specific orbital type was summed up and then assigned to the orbitals in accordance with the Aufbau principle.

The resulting partial charge of the iodine atom can be converted to resemble the theoretical ionic character I_theo, to enable a comparison with the results from the extended Townes–Dailey model. This comparison is shown in Fig. 8, both with and without the treatment of relativistic effects employing DKH2. The results are organised into three different groups, depending on the amount of ortho-positioned fluorine atoms. In all cases, I_theo is positive, indicating a negative partial charge on the iodine, caused by an electron-donating effect of the aromatic substitute, regardless of the substitution pattern. Explicit treatment of relativistic effects leads to a maximum reduction in the ionic character of 0.93% and therefore a stronger electron-withdrawing effect. This consistently positive ionic character differs fundamentally from I_TD, where positive and negative results are observed. Apart from this difference in absolute references, linear trends are evident within each of the three groups for both levels of theory, as indicated by the dashed lines. Larger slopes are obtained with DKH2 for each group of ortho-positioned fluorine atoms. The group with two fluorine atoms in the ortho position has the largest slope, while the group without fluorine atoms in the ortho position has the smallest slope.

Fig. 8 The ionic characters I_theo of the measured fluoroiodobenzenes and iodobenzene, obtained from an intrinsic basis bonding analysis (IBBA) at the B3LYP-D3(BJ,abc)/cc-pVTZ-DK level of theory with (purple) and without (blue) the treatment of relativistic effects using DKH2, as described in Section 2.3, plotted against the ionic characters I_TD on the abscissa. These were obtained from the implementation of the extended Townes–Dailey model, as described in Section 2.4.2. The structures were optimised at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. Following the nomenclature in Fig. 1, the data for the molecules 0, 2, 3, 4, 23, 24, 25, 26, 34, 35, 234, 235, 236, 245, 246, 345, 2345, 2346, 2356, and 23456 are shown. Full circles represent the presence of two ortho-positioned fluorine atoms, half-filled circles represent one ortho-positioned fluorine atom, and empty circles represent the absence of fluorine atoms. Dashed lines indicate linear trends within these groups of molecules and are discussed in the text. Exact values are summarised in Table S13 of the SI. Uncertainties of I_TD are not shown.

In the case of iodine, the extended Townes–Dailey model's ionic character is obtained by deriving the interactions of the valence 5p-orbitals from the measured NQCCs. Therefore, all chemical information encoded in the tensor is projected onto the valence p-orbitals only. To analyse the feasibility of this restriction, valence occupations of iodine's 5p-orbitals, n_x,theo, n_y,theo and n_z,theo, were extracted from the IBBA analysis and compared with their experimental counterparts, n_x, n_y and n_z, following the extended Townes–Dailey model. This is illustrated in Fig. 9 for each combination. The use of DKH2 results in smaller occupations among all 5p-orbitals. This agrees with the findings of less positive ionic characters. However, the differences in withdrawn 5p-orbital electrons are consistent with and without relativistic corrections, while the differences in ionic characters change and become especially small with a minimum difference of 0.074% for 0 in the group without ortho-positioned fluorine atoms. One possible explanation for these results is that relativistic effects lead to the contraction of s-orbitals, hindering their interaction with the aromatic substituent in terms of electron donation to the iodine atom. Consequently, in the case of DKH2, the total electron-donating effect of the aromatic substituent is mediated via the p-orbitals, which may limit the amount of electron density that can be donated due to a bottlenecking effect. Without the relativistic treatment, the non-contracted s-orbitals are also available, in addition to the p-orbitals. This leads to more similar ionic characters, as the s-orbital interaction balances out the smaller p-orbital interaction. To further investigate this explanation, Fig. 10 shows the deviation in electron occupation of all s, p- and d-orbitals compared to the isolated iodine atom with 10 s-electrons, 23 p-electrons and 20 d-electrons. This resembles the strength of interaction for the different orbital types. Differences in deviations between both the treatment of relativistic effects and no treatment, are calculated. In both cases, the s-orbitals withdraw electron density from the aromatic substituent. In contrast, the p- and d-orbitals donate electron density, leading to a push–pull scenario. When DKH2 is employed, over 0.12 less electrons are withdrawn via the s-orbitals and interaction via the p-orbitals is also reduced by over 0.11 electrons. The decrease in electron density donated via the p-orbitals balances out the decrease in electrons withdrawn via the s-orbitals. This results in a similar total ionic character for DKH2 and no treatment. For the d-orbitals, only a slight decrease in donated electron density of approximately 0.01 is observed. Within this theoretical framework, the treatment of the relativistic effect for iodine is crucial, as the results on the orbital occupation deviations differ significantly with and without such treatment. It should be noted that the chosen basis set was specifically contracted for DKH2, which may affect the results when no such relativistic treatment is employed.

Fig. 9 The 5p-orbital occupations of iodine, n_x,theo, n_y,theo and n_z,theo, of the measured fluoroiodobenzenes and iodobenzene, obtained from an intrinsic basis bonding analysis (IBBA) at the B3LYP-D3(BJ,abc)/cc-pVTZ-DK level of theory with (purple) and without (blue) the treatment of relativistic effects using DKH2, as described in Section 2.3, plotted against the respective valence p-orbital occupation n_x, n_y and n_z on the abscissa, obtained from the implementation of the extended Townes–Dailey model, as described in Section 2.4.2. Structures have been optimised at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. Following the nomenclature from Fig. 1, the data for the molecules 0, 2, 3, 4, 23, 24, 25, 26, 34, 35, 234, 235, 236, 245, 246, 345, 2345, 2346, 2356, and 23 456 are shown. Full circles represent the presence of two ortho-positioned fluorine atoms, half-filled circles represent one ortho-positioned fluorine atom, and empty circles represent the absence of ortho-positioned fluorine atoms. Dashed lines indicate linear trends and are discussed in the text. Exact values are summarised in Tables S11 and S14 of the SI.

Fig. 10 The differences in deviations Δ(Δn_theo) of the iodine atom's orbital occupations, Δn_theo, to the isolated iodine atom (10 electrons in s orbitals, 23 electrons in p orbitals, 20 electrons in d orbitals), for the fluoroiodobenzenes and iodobenzene from an intrinsic basis bonding analysis (IBBA), performed at the B3LYP-D3(BJ,abc)/cc-pVTZ-DK level of theory with versus without the treatment of relativistic effects using DKH2. Results without relativistic treatments were subtracted from those with DKH2 enabled. In all cases, absolute deviations were larger for the treatment of relativistic effects. A positive sign indicates a gain in electron density. Structures have been optimised at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. Following the nomenclature from Fig. 1, the data for all different substitution patterns is shown. Orbital types are indicated by the colours shown in the legend. The molecules are sorted in ascending order of p deviation (identical in both cases). Exact values are summarised in Table S15 of the SI.

Apart from these deviations in the absolute 5p-orbital occupations with and without DKH2, trends in the comparison to the extended Townes–Dailey model are identical in both cases. For n_x, which is the p-orbital perpendicular to the aromatic ring, a linear correlation can be observed, with the previously defined groups of molecules being separated. For the p-orbital in the y-direction, lying in the ring's plane, n_y has a constant value of two, indicating no interaction for this specific orbital. However, the IBBA analysis disagrees with these findings, as similar interactions to n_x,theo are observed for those of n_y,theo. The majority of the electron density is transferred via the 5p-orbital along the C–I bond for both the calculation and the model. Once again, groups are separated by the number of ortho fluorine atoms, with a general linear correlation being observable.

The disagreement between the extended Townes–Dailey model and the IBBA regarding the occupation of the valence 5p-orbitals n_y and n_y,theo indicates that the interactions of the iodine atom are not adequately described by the valence p-orbitals alone. However, for the total ionic character, these deviations balance out, leading to a correlation between the calculation and the model. As shown in Fig. 10, the s- and d-orbitals of the iodine atom contribute substantially in the chemical bonding by exchanging electron density with the aromatic substituents. The total effect is mapped onto the valence p-orbitals, so the extended Townes–Dailey model can only describe these interactions implicitly. Therefore, even though the total interaction expressed via the ionic characters is restored, the limitation to just the valence p-orbitals may no longer be viable for heavy atoms such as iodine. It can also not be ruled out that other bonding analysis models may be in agreement with the extended Townes–Dailey model.

In order to investigate the claimed correlation between reduced valence p_z-orbital occupation and enhanced halogen bond donor ability, resulting in a stronger halogen bond, the planar complexes with the acceptor pyridine were calculated for all 19 fluoroiodobenzenes and iodobenzene. These calculations were performed at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. Although no imaginary frequencies were observed for any of these planar structures, depending on the substitution pattern, the halogen-bonded complexes are not always the global minimum structures. For a simple comparison, Fig. 11 plots experimentally obtained occupations n_z against theoretical zero-point-corrected dissociation energies D₀. While stabilising secondary effects, such as attractive fluorine–hydrogen interactions, cannot be ruled out and are possibly preventing a global correlation, linear trends are obtained within each group consisting of donors with the same number of fluorine atoms in the ortho positions. In this simple comparison, a decreased n_z does lead to a stronger halogen bond as indicated by the larger dissociation energy. Therefore, the valence orbital occupations following our implementation of the extended Townes–Dailey model may be used as a metric to gauge a molecule's halogen bond donor ability. The total ionic character, on the other hand, also includes π interactions, which should not affect the formation of the σ-hole. As seen above, these contributions are ca. 3.6 times weaker than the σ interaction. For this reason, a similar correlation can be observed between D₀ and I_TD, albeit to a lesser extent. This is shown in Fig. S2 of the SI. Therefore, the ionic character may also be used in order to gauge halogen bond donor ability, provided π interactions are negligible.

Fig. 11 Calculated zero-point-corrected dissociation energies D₀ of the planar halogen-bonded complexes of pyridine with the fluoroiodobenzenes, including iodobenzene, plotted against the valence p_z-orbital occupations n_z following the implementation of the extended Townes–Dailey model. Calculations were performed at the B3LYP-D3(BJ,abc)/ma-def2-TZVP level of theory. Full circles represent two ortho-positioned fluorine atoms, half-filled circles represent one ortho-positioned fluorine atom, and empty circles represent an absence of ortho-positioned fluorine atoms. Dashed lines indicate linear trends and are discussed in the text. Exact values are summarised in Table S11 and S18 of the SI.

4 Summary and conclusions

19 iodine-containing monomers were characterised both theoretically and experimentally. The main goal was to systematically study properties beneficial to halogen bond donors and how these properties can be influenced. Fluoroiodobenzenes were chosen for this study because they allow the investigation of both inductive and mesomeric effects mediated by the number and position of the fluorine atoms. The iodine atom was used as a probe in order to gauge these effects.

To derive ionic characters and p-orbital occupations from experimentally measured nuclear quadrupole coupling constant tensors, the extended Townes–Dailey model proposed by Novick¹⁶ was systematically implemented, ensuring the reliability of the results within the model's limitations. It was revealed that, in terms of the decrease in ionic character, the position of the fluorine atoms was more important than their number for the electron-withdrawing effect of the aromatic substituent on the iodine atom.

The experimental results were compared to theory in the form of an intrinsic basis bonding analysis. This revealed good agreement between the ionic characters, but a deviation in the p-orbital occupations. Contrary to the extended Townes–Dailey model, both p-orbitals perpendicular to the C–I bond interacted with the aromatic ring in the same manner. However, the extended Townes–Dailey model only predicted an interaction for the p-orbital perpendicular to the aromatic ring. This indicated that the extended Townes–Dailey model, limited to the 5p orbitals for iodine, may not be adequate for describing these interactions. Nevertheless, the total trends in ionic characters between the two approaches matched. The importance of relativistic corrections was revealed, as the IBBA results differed when these corrections were employed. One possible explanation for this was the bottlenecking effect of the now contracted s-orbitals, which led to less electron density being withdrawn from the substituent in the relativistic framework.

Additionally, theory was benchmarked against the wide range of experimental results obtained. Rotational constants were benchmarked within the group of fluoroiodobenzenes and iodobenzene via VPT2 calculations. This comparison highlighted the use cases of scaling within such a group of similar geometries. A training set was defined that included a representative of each degree of fluorine substitution. When comparing experiment and theory for all measured fluoroiodobenzenes, scaled B3LYP-D3(BJ,abc)/ma-def2-TZVP calculations outperformed all zero-point corrected ab initio approaches, including the best-performing DSD-PBEP86-D3(BJ)/def2-TZVP. It is not possible to give a definitive answer as to how strongly this result was influenced by the self-referencing problem of our test set, which only applies to the scaled approach. However, in the comparison of experimental rotational constants with equilibrium results calculated at an equivalent level of theory, the scaled approach was shown to be surpassed by B2PLYP-D3(BJ), an approach that demonstrated better agreement with the experiment when zero-point effects were excluded. A similar result was observed for B3LYP-D3(BJ), which exhibited an enhanced performance, albeit not reaching the level of DSD-PBEP86-D3(BJ)/def2-TZVP with zero-point correction. Hence, when comparing the equilibrium values to the experimental ones, B3LYP and B2PLYP can be considered right for the wrong reasons whereas DSD-PBEP86 is wrong for the right reason. If the zero point corrected values are compared DSD-PBEP86 is right for the right reasons while B3LYP and B2PLYP are wrong for the wrong reasons.

In summary, the available data on iodine-containing systems studied via rotational spectroscopy has substantially increased. At the time of writing, a total of 95 molecules and conformers are reported in ref. 76. In this study, 19 new molecules were characterised, reflecting an increase of approximately 20%. Theory was rigorously benchmarked against experiment, highlighting use cases for parametrisations with experimental results. The halogen bond donor properties of fluoroiodobenzenes were systematically investigated and probed via NQCC, employing an implementation of the extended Townes–Dailey model.

Future work could involve studying groups of fluorinated bromo- and chlorobenzenes, as well as the complexes formed by these donors, in a similar manner to that presented here. This would allow us to evaluate experimental results using the extended Townes–Dailey model in order to gauge donor/acceptor interactions, as well as investigate the importance of relativistic effects for lighter homologues through rigorous benchmarking. As the present test focused on rather rigid systems, the effect of projection errors cannot be conclusively determined. The systematic study of fluoroiodobenzenes and iodobenzene presented here paves the way for a similar investigation and benchmarking of a series of halogen-bonded complexes to that of Dohmen et al.,²⁷ incorporating relativistic and non-relativistic corrections.

Author contributions

Conceptualisation: JNH, BH, and DAO; data curation: JNH and BH; formal analysis: JNH, BH, and DAO; funding acquisition: DAO; investigation: JNH and BH; supervision: BH and DAO; visualisation: JNH and BH; writing – original draft: JNH; synthesis: JFK; writing – review & editing: JNH, BH, JFK, and DAO.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data have been made available in separate data repositories. These repositories have their own DOIs. They include the raw microwave spectra with electronic laboratory notebooks (https://doi.org/10.25625/KSE9WG), all computational outputs (https://doi.org/10.25625/7NTRX3) as well as the outputs of the experimental fits (https://doi.org/10.25625/BVGMII).^121–123 All python scripts used for data processing are available upon request.

Supplementary information (SI): the exact values corresponding to the figures presented in this paper. Additional experimental details for the synthesis of both 2,3,4,5-tetrafluoroiodobenzene and 2,3,4,6-tetrafluoroiodobenzene, as well as for the microwave experiments. We also provide a summary of the different levels of theory used throughout this work, alongside exemplary input files. An overview of all experimental fits is also provided. See DOI: https://doi.org/10.1039/d6cp00032k.

Acknowledgements

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 389479699/GRK2455 and 405832858 for computational resources. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – OB 535/1-1. This work used the Scientific Computer Cluster at the GWDG, the joint data center of the Max Planck Society for the Advancement of Science (MPG) and the University of Göttingen. We thank Prof. Dr. Martin A. Suhm for interesting discussions and valuable input and Prof. Dr. Sonia Melandri for helpful discussions regarding the measurement of the pentafluoroiodobenzene monomer. We also thank Prof. Dr. Ricardo A. Mata for valuable discussions about the implementation of IBBA in Molpro. Finally, we thank Prof. Dr. Manuel Alcarazo for the support with the syntheses of 2,3,4,5-tetrafluoroiodobenzene and 2,3,4,6-tetrafluoroiodobenzene. During the preparation of this work the authors used AI editing tools for proof reading and grammar in early drafts of the text. Additionally it was used to create/improve Python scripts used for data evaluation and figure creation. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of this publication.

Notes and references

1. G. Cavallo, P. Metrangolo, R. Milani, T. Pilati, A. Priimagi, G. Resnati and G. Terraneo, Chem. Rev., 2016, 116, 2478–2601.
2. G. R. Desiraju, P. S. Ho, L. Kloo, A. C. Legon, R. Marquardt, P. Metrangolo, P. Politzer, G. Resnati and K. Rissanen, Pure Appl. Chem., 2013, 85, 1711–1713.
3. M. M. Colin and H. Gaultier de Claubry, Ann. Chim., 1814, 90, 87–100.
4. M. Colin, Ann. Chim., 1814, 91, 252–272.
5. L. Fotović and V. Stilinović, CrystEngComm, 2020, 22, 4039–4046.
6. J. B. Davey, A. C. Legon and E. R. Waclawik, Phys. Chem. Chem. Phys., 2000, 2, 1659–1665.
7. E. R. Waclawik and A. C. Legon, Phys. Chem. Chem. Phys., 1999, 1, 4695–4700.
8. S. L. Stephens, N. R. Walker and A. C. Legon, Phys. Chem. Chem. Phys., 2011, 13, 21093–21101.
9. S. L. Stephens, N. R. Walker and A. C. Legon, Phys. Chem. Chem. Phys., 2011, 13, 20736–20744.
10. P. Metrangolo and G. Resnati, Chem. – Eur. J., 2001, 7, 2511–2519.
11. R. Bertani, P. Metrangolo, A. Moiana, E. Perez, T. Pilati, G. Resnati, I. Rico-Lattes and A. Sassi, Adv. Mater., 2002, 14, 1197–1201.
12. H. L. Nguyen, P. N. Horton, M. B. Hursthouse, A. C. Legon and D. W. Bruce, J. Am. Chem. Soc., 2004, 126, 16–17.
13. Z. Xu, Z. Yang, Y. Liu, Y. Lu, K. Chen and W. Zhu, J. Chem. Inf. Model., 2014, 54, 69–78.
14. F. Kniep, S. H. Jungbauer, Q. Zhang, S. M. Walter, S. Schindler, I. Schnapperelle, E. Herdtweck and S. M. Huber, Angew. Chem., Int. Ed., 2013, 52, 7028–7032.
15. C. H. Townes and B. P. Dailey, J. Chem. Phys., 1949, 17, 782–796.
16. S. E. Novick, J. Mol. Spectrosc., 2011, 267, 13–18.
17. D. Lv, A. Maris, L. Evangelisti, A. Maggio, W. Song, A. A. Elliott, S. A. Peebles, J. L. Neill, M. T. Muckle, B. H. Pate, R. A. Peebles and S. Melandri, Phys. Chem. Chem. Phys., 2021, 23, 18093–18101.
18. K. Johansson, H. Oldenberg and R. Selen, Ark. Fys., 1965, 29, 531–532.
19. A. Mirri and W. Caminati, Chem. Phys. Lett., 1971, 8, 409–412.
20. K. Etchison, C. Dewberry, K. Kerr, D. Shoup and S. Cooke, J. Mol. Spectrosc., 2007, 242, 39–45.
21. O. Dorosh, E. Białkowska-Jaworska, Z. Kisiel and L. Pszczółkowski, J. Mol. Spectrosc., 2007, 246, 228–232.
22. G. J. P. Perry, J. M. Quibell, A. Panigrahi and I. Larrosa, J. Am. Chem. Soc., 2017, 139, 11527–11536.
23. G. Deacon and D. Tunaley, Aust. J. Chem., 1979, 32, 737–753.
24. D. Bulfield and S. M. Huber, J. Org. Chem., 2017, 82, 13188–13203.
25. T. J. Balle and W. H. Flygare, Rev. Sci. Instrum., 1981, 52, 33–45.
26. U. Andresen, H. Dreizler, U. Kretschmer, W. Stahl and C. Thomsen, Fresenius’ J. Anal. Chem., 1994, 349, 272–276.
27. R. Dohmen, S. Arnold, J. Garrett, B. Kempken, B. Hartwig, B. Schröder, P. Pinacho, M. Schnell, G. G. Brown and D. A. Obenchain, J. Phys. Chem. A, 2025, 129, 860–873.
28. M. Schnell, D. Banser and J.-U. Grabow, Rev. Sci. Instrum., 2004, 75, 2111–2115.
29. H. M. Pickett, J. Mol. Spectrosc., 1991, 148, 371–377.
30. H. M. Pickett, J. Mol. Spectrosc., 2004, 228, 659–663.
31. H. Pickett, J. Pearson and C. Miller, J. Mol. Spectrosc., 2005, 233, 174–179.
32. Z. Kisiel, L. Pszczółkowski, I. R. Medvedev, M. Winnewisser, F. C. De Lucia and E. Herbst, J. Mol. Spectrosc., 2005, 233, 231–243.
33. Z. Kisiel, L. Pszczółkowski, B. J. Drouin, C. S. Brauer, S. Yu, J. C. Pearson, I. R. Medvedev, S. Fortman and C. Neese, J. Mol. Spectrosc., 2012, 280, 134–144.
34. Z. Kisiel, PROSPE – Programs for ROtational SPEctroscopy – info.ifpan.edu.pl, https://info.ifpan.edu.pl/kisiel/prospe.htm, 2001.
35. G. Winnewisser, J. Chem. Phys., 1972, 56, 2944–2954.
36. B. van Eijck, J. Mol. Spectrosc., 1974, 53, 246–249.
37. J. K. G. Watson, Vibrational spectra and structure, Elsevier Scientific Publishing Company, Amsterdam, 1977, vol. 6, pp. 1–89.
38. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
39. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108.
40. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
41. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
42. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, I. Miller, F. Thomas, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
43. H.-J. Werner, P. J. Knowles, P. Celani, W. Györffy, A. Hesselmann, D. Kats, G. Knizia, A. Köhn, T. Korona, D. Kreplin, R. Lindh, Q. Ma, F. R. Manby, A. Mitrushenkov, G. Rauhut, M. Schütz, K. R. Shamasundar, T. B. Adler, R. D. Amos, S. J. Bennie, A. Bernhardsson, A. Berning, J. A. Black, P. J. Bygrave, R. Cimiraglia, D. L. Cooper, D. Coughtrie, M. J. O. Deegan, A. J. Dobbyn, K. Doll, M. Dornbach, F. Eckert, S. Erfort, E. Goll, C. Hampel, G. Hetzer, J. G. Hill, M. Hodges, T. Hrenar, G. Jansen, C. Köppl, C. Kollmar, S. J. R. Lee, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, B. Mussard, S. J. McNicholas, W. Meyer, T. F. Miller III, M. E. Mura, A. Nicklass, D. P. O'Neill, P. Palmieri, D. Peng, K. A. Peterson, K. Pflüger, R. Pitzer, I. Polyak, M. Reiher, J. O. Richardson, J. B. Robinson, B. Schröder, M. Schwilk, T. Shiozaki, M. Sibaev, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, J. Toulouse, M. Wang, M. Welborn and B. Ziegler, MOLPRO, version 2024.1, a package of ab initio programs, 2024, see https://www.molpro.net.
44. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian∼16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
45. H. H. Nielsen, Rev. Mod. Phys., 1951, 23, 90–136.
46. V. Barone, J. Chem. Phys., 2004, 122, 014108.
47. A. Kanzow, B. Hartwig, P. Buschmann, K. G. Lengsfeld, C. M. Höhne, J. S. Hoke, F. Knüppe, F. Köster, J. K. van Spronsen, J.-U. Grabow, D. McNaughton and D. A. Obenchain, Phys. Chem. Chem. Phys., 2025, 27, 9491–9503.
48. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
49. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
50. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
51. B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett., 1989, 157, 200–206.
52. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
53. M. Ernzerhof and G. E. Scuseria, J. Chem. Phys., 1999, 110, 5029–5036.
54. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
55. H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540–3544.
56. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
57. O. A. Vydrov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 234109.
58. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2007, 120, 215–241.
59. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
60. T. Schwabe and S. Grimme, Phys. Chem. Chem. Phys., 2007, 9, 3397–3406.
61. S. Kozuch and J. M. L. Martin, J. Comput. Chem., 2013, 34, 2327–2344.
62. S. Kozuch and J. M. L. Martin, Phys. Chem. Chem. Phys., 2011, 13, 20104–20107.
63. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
64. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
65. S. Sæbø and J. Almlöf, Chem. Phys. Lett., 1989, 154, 83–89.
66. M. J. Frisch, M. Head-Gordon and J. A. Pople, Chem. Phys. Lett., 1990, 166, 275–280.
67. M. J. Frisch, M. Head-Gordon and J. A. Pople, Chem. Phys. Lett., 1990, 166, 281–289.
68. M. Head-Gordon and T. Head-Gordon, Chem. Phys. Lett., 1994, 220, 122–128.
69. S. Grimme, J. Chem. Phys., 2003, 118, 9095–9102.
70. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
71. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
72. J. Zheng, X. Xu and D. Truhlar, Theor. Chem. Acc., 2011, 128, 295–305.
73. B. M. Axilrod and E. Teller, J. Chem. Phys., 1943, 11, 299–300.
74. Y. Muto, Proc. Phys.-Math. Soc. Jpn., 1943, 17, 629–631.
75. P. Pyykkö, Mol. Phys., 2018, 116, 1328–1338.
76. W. C. Bailey, Iodine – nqcc.wcbailey.net, https://nqcc.wcbailey.net/Iodine.html, Accessed 05-08-2025.
77. M. N. Glukhovtsev, A. Pross, M. P. McGrath and L. Radom, J. Chem. Phys., 1995, 103, 1878–1885.
78. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
79. K. A. Peterson, D. Figgen, E. Goll, H. Stoll and M. Dolg, J. Chem. Phys., 2003, 119, 11113–11123.
80. P. Pollak and F. Weigend, J. Chem. Theory Comput., 2017, 13, 3696–3705.
81. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822–8824.
82. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 7406.
83. A. D. Becke, J. Chem. Phys., 1997, 107, 8554–8560.
84. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2015, 142, 074111.
85. B. A. Hess, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 33, 3742–3748.
86. M. Douglas and N. M. Kroll, Ann. Phys., 1974, 82, 89–155.
87. G. Jansen and B. A. Hess, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 6016–6017.
88. L. Visscher and K. Dyall, At. Data Nucl. Data Tables, 1997, 67, 207–224.
89. F. Neese, A. Wolf, T. Fleig, M. Reiher and B. A. Hess, J. Chem. Phys., 2005, 122, 204107.
90. D. Peng, N. Middendorf, F. Weigend and M. Reiher, J. Chem. Phys., 2013, 138, 184105.
91. Y. J. Franzke, F. Mack and F. Weigend, J. Chem. Theory Comput., 2021, 17, 3974–3994.
92. R. Dohmen, D. Fedosov and D. A. Obenchain, Phys. Chem. Chem. Phys., 2023, 25, 2420–2429.
93. W. Sun, D. Kargin, R. Pietschnig and M. Schnell, Mol. Phys., 2025, 0, e2585903.
94. G. Knizia, J. Chem. Theory Comput., 2013, 9, 4834–4843.
95. J. Dunning and H. Thom, J. Chem. Phys., 1989, 90, 1007–1023.
96. W. A. de Jong, R. J. Harrison and D. A. Dixon, J. Chem. Phys., 2001, 114, 48–53.
97. B. P. Dailey and C. H. Townes, J. Chem. Phys., 1955, 23, 118–123.
98. V. Jaccarino, J. G. King, R. A. Satten and H. H. Stroke, Phys. Rev., 1954, 94, 1798–1799.
99. W. Gordy and R. L. Cook, Microwave molecular spectra, Wiley, New York, 3rd edn, 1984.
100. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.
101. CRC Handbook of Chemistry and Physics, ed. W. M. Haynes, CRC Press, Philadelphia, PA, 97th edn, 2016.
102. S. A. Peebles and R. A. Peebles, J. Mol. Struct., 2002, 607, 19–29.
103. R. K. Bohn, J. A. Montgomery, H. H. Michels and J. A. Fournier, J. Mol. Spectrosc., 2016, 325, 42–49.
104. S. Stopkowicz, L. Cheng, M. E. Harding, C. Puzzarini and J. Gauss, Mol. Phys., 2013, 111, 1382–1389.
105. K. Chance, T. Varberg, K. Park and L. Zink, J. Mol. Spectrosc., 1993, 162, 120–126.
106. B. Osipov and M. N. Grabois, J. Mol. Spectrosc., 1985, 111, 344–351.
107. Z. Kisiel, L. Pszczólkowski, W. Caminati and P. G. Favero, J. Chem. Phys., 1996, 105, 1778–1785.
108. T. Inagusa, M. Fujitake and M. Hayashi, J. Mol. Spectrosc., 1988, 128, 456–468.
109. H.-O. Munz, H.-K. Bodenseh and M. Ferner, J. Mol. Struct., 2004, 695–696, 189–202.
110. H. S. P. Müller and M. C. L. Gerry, J. Chem. Phys., 1995, 103, 577–583.
111. J. Gripp, H. Dreizler, J. Gadhi, G. Wlodarczak, J. Legrand, J. Burie and J. Demaison, J. Mol. Spectrosc., 1988, 129, 381–387.
112. H. Ozeki and S. Saito, J. Chem. Phys., 2004, 120, 5110–5116.
113. J. Gripp and H. Dreizler, Z. Naturforsch., A, 1990, 45, 715–723.
114. M. Hayashi, C. Ikeda and T. Inagusa, J. Mol. Spectrosc., 1990, 139, 299–312.
115. M. Fujitake and M. Hayashi, J. Mol. Spectrosc., 1988, 128, 444–455.
116. A. Yokozeki and J. S. Muenter, J. Chem. Phys., 1980, 72, 3796–3804.
117. U. Andresen, N. Heineking and H. Dreizler, J. Mol. Spectrosc., 1989, 137, 296–299.
118. C. E. Miller and E. A. Cohen, J. Chem. Phys., 2003, 118, 6309–6317.
119. S. Bailleux, H. Ozeki, S. Sakai, T. Okabayashi, P. Kania and D. Duflot, J. Mol. Spectrosc., 2011, 270, 51–55.
120. I. Merke, A. Lüchow and W. Stahl, J. Mol. Struct., 2006, 780–781, 295–299.
121. J. N. Hasselhorn, B. Hartwig, D. A. Obenchain and J. F. Köster, Microwave Spectra for the Fluorinated Iodobenzenes, GRO.data, 2026, V1 DOI:10.25625/KSE9WG.
122. J. N. Hasselhorn, B. Hartwig, D. A. Obenchain and J. F. Köster, Computational Outputs of Geometry Optimisations, Frequency Calculations, Intrinsic Basis Bonding Analysis and Calculations of the Electric Field Gradients, GRO.data, 2026, V1 DOI:10.25625/7NTRX3.
123. J. N. Hasselhorn, B. Hartwig, D. A. Obenchain and J. F. Köster, Experimental Pickett Fits of the Fluoroiodobenzenes, GRO.data, 2026, V1 DOI:10.25625/BVGMII.

---