A Quantitative Molecular Orbital Perspective of the Chalcogen Bond

Abstract We have quantum chemically analyzed the structure and stability of archetypal chalcogen‐bonded model complexes D2Ch⋅⋅⋅A− (Ch = O, S, Se, Te; D, A = F, Cl, Br) using relativistic density functional theory at ZORA‐M06/QZ4P. Our purpose is twofold: (i) to compute accurate trends in chalcogen‐bond strength based on a set of consistent data; and (ii) to rationalize these trends in terms of detailed analyses of the bonding mechanism based on quantitative Kohn‐Sham molecular orbital (KS‐MO) theory in combination with a canonical energy decomposition analysis (EDA). At odds with the commonly accepted view of chalcogen bonding as a predominantly electrostatic phenomenon, we find that chalcogen bonds, just as hydrogen and halogen bonds, have a significant covalent character stemming from strong HOMO−LUMO interactions. Besides providing significantly to the bond strength, these orbital interactions are also manifested by the structural distortions they induce as well as the associated charge transfer from A− to D2Ch.


Introduction
The chalcogen-bond (ChB) is the net-attractive intermolecular interaction, often referred to as noncovalent interaction, between a Lewis-basic chalcogen-bond acceptor A and a Lewisacidic chalcogen-bond donor D 2 Ch featuring a chalcogen (group 16) atom Ch to which A binds. [1] Nearly 40 years ago, the first systematic study appeared of the chalcogen bond in which S···Y (e. g., Y = S, O, F, Cl, or Br) nonbonded atomic contacts were investigated. [2] Early studies generally characterized chalcogen bonds as being predominantly electrostatic in nature. [3] Later on, the significance of charge transfer from the occupied orbital of a Lewis base into an empty σ*-type orbital of a chalcogen molecule controlling the chalcogen bond strength was recognized. [4] Chalcogen-bonding has since found applications in various fields of chemistry, [1] including, supramolecular, [5] biochemistry, [6] spectroscopy [7] and catalysis. [8] In this study, we have computationally analyzed a range of chalcogen-bonded D 2 Ch···A À complexes (Ch = O, S, Se, Te; D, A = F, Cl, Br; see Scheme 1), using relativistic density functional theory (DFT) at ZORA-M06/QZ4P. One purpose of our work is to provide a set of consistent structural and energy data from which reliable trends can be inferred for a wide range of model systems. The primary objective is to achieve a detailed understanding of the nature of chalcogen bonds by studying the associated electronic structure and bonding mechanism and compare them with the better-known halogen bonds and hydrogen bonds. [9] To this end, we first explore how the geometries and energies of our model complexes D 2 Ch···A À vary as the chalcogen atom (Ch), or the chalcogen bond accepting Lewis base (A À ) are varied. To understand the origin of the computed trends, activation strain analyses [10] are performed on the formation of the chalcogen-bond complexes. As part of these analyses, the interaction energy and the underlying bonding mechanism are furthermore examined in the context of quantitative Kohn-Sham molecular orbital (MO) theory in combination with an energy decomposition analysis (EDA). [11,12] Our systematic and detailed analyses along the entire reaction profile for each of the chalcogen-bond complexation reactions provide in-depth insights. In particular, they demonstrate that chalcogen bonds are not at all purely electrostatic phenomena but are, to a substantial extent, covalent in nature.

Computational Details
All calculations were carried out using the Amsterdam Density Functional (ADF) 2017.103 program. [13] The equilibrium geometries and energies of chalcogen-bonded complexes were computed at DFT level using the meta-hybrid functional M06. [14] In addition, a large uncontracted relativistically optimized QZ4P Slater type orbitals (STOs) basis set containing diffuse functions was used. The QZ4P all-electron basis set, [15] no frozen-core approximation, is of quadruple-ζ quality for all atoms and has been augmented with the following sets of polarization and diffuse functions: two 3 d and two 4 f on oxygen and fluorine, three 3 d and two 4 f on sulfur and chlorine, two 4 d and three 4 f on selenium and bromine, one 5 d and three 4 f on tellurium and iodine. The molecular density was fitted by the systematically improvable Zlm fitting scheme. The scalar relativistic effects were accounted for by using the zeroth-order regular approximation (ZORA) Hamiltonian. [16] It has been shown that these computational settings give accurate bond lengths and energies. [17]

Analysis of the Bonding Mechanism
Insight into the bonding mechanism is obtained through activation strain analyses of the various chalcogen bond formation reactions. These complexation reactions are computationally modeled by increasing the distance between A À and the Ch atom of the D 2 Ch fragment, allowing the system to geometrically relax at each point. The D 2 Ch···A À distance is increased from the actual bond length value in the chalcogenbonded complex (r Ch···A À ) to a value of 4.300 Å. Thus, each analysis starts from an optimized D 2 Ch···A À complex, which is then transformed to the D 2 Ch molecule and a halide at a relatively large distance.
These complexation reactions are analyzed using the activation strain model. The activation strain model of chemical reactivity [10] is a fragment-based approach to understand the energy profile of a chemical process in terms of the original reactants. Thus, the potential energy surface ΔE(ζ) is decomposed along the reaction coordinate ζ (or just at one point along ζ) into the strain energy ΔE strain (ζ), which is associated with the geometrical deformation of the individual reactants as the process takes place, plus the actual interaction energy ΔE int (ζ) between the deformed reactants [Eq. (1)].
In the equilibrium geometry, that is, for ζ = ζ eq , this yields an expression for the bond energy ΔE(ζ eq ) = ΔE strain (ζ eq ) + ΔE int (ζ eq ). The PyFrag program was used to facilitate the analyses along the reaction coordinate ζ of the bond formation processes. [18] The interaction energy ΔE int (ζ) between the deformed reactants is further analyzed in the conceptual framework provided by the quantitative Kohn-Sham MO model. [11] To this end, it is decomposed in three physically meaningful terms [Eq. (2)] using a quantitative energy decomposition scheme developed by Ziegler and Rauk. [12] DE int ðzÞ ¼ DV elstat ðzÞ þ DE Pauli ðzÞ þ DE oi ðzÞ (2) The usually attractive term ΔV elstat corresponds to the classical Coulomb interaction between the unperturbed charge distributions of the deformed reactants and has four components [Eq. (3)]: i) the electrostatic repulsion between the electron densities of fragments 1 and 2, ΔV elstat,1112 ; ii) the electrostatic attraction between the nucleus of fragment 1 and the electron density of fragment 2, ΔV elstat,n112 ; iii) the electrostatic attraction between the electron density of fragment 1 and the nucleus of fragment 2, ΔV elstat,11n2 ; and iv) the electrostatic repulsion between the nuclei of fragments 1 and 2, ΔV elstat,n1n2 . The Pauli repulsion energy (ΔE Pauli ) comprises the destabilizing interactions between occupied orbitals of the reactants and is responsible for steric repulsion. The orbitalinteraction energy (ΔE oi ) accounts for charge transfer, that is, the interaction between occupied orbitals of one fragment with unoccupied orbitals of the other fragment, including the interactions of the highest occupied and lowest unoccupied MOs (HOMOÀ LUMO), and polarization, that is, empty-occupied orbital mixing on one fragment, due to the presence of another fragment.
The electron density distribution is analyzed using the Voronoi deformation density (VDD) method for computing atomic charges. [19] The VDD atomic charge on atom X in a molecule (Q X VDD ) is computed as the (numerical) integral of the deformation density in the volume of the Voronoi cell of atom X [Eq. (4)]. The Voronoi cell of atom X is defined as the compartment of space bounded by the bond midplanes on and perpendicular to all bond axes between nucleus X and its neighboring nuclei.
Here, the deformation density is the difference between 1(r), i. e., the electron density of the overall molecule or complex, and 1 promolecule (r) = Σ Y 1 Y (r), i. e., the superposition of spherical average-of-configuration atomic densities 1 Y (r) of each atom Y in the fictitious promolecule without chemical interactions, in which all atoms are considered neutral. The interpretation of the VDD charge Q Ch VDD is rather straightforward and transparent: instead of measuring the amount of charge associated with a particular atom Ch, Q Ch VDD directly monitors how much charge flows out of (Q Ch VDD > 0) or into (Q Ch VDD < 0) the Voronoi cell of atom Ch due to chemical interactions.
The VDD scheme can also be used to directly compute how much charge flows into or out of an atomic Voronoi cell X in an overall complex (e. g., [D 2 Ch···A] À ) relative to two (poly)atomic molecular fragments (e. g., D 2 Ch and A À ), instead of spherical atoms, as shown in [Eq. (5)].
ΔQ X VDD is a measure of how the atomic charge of atom X changes due to the bonding between the fragment. In this work, [Eq. (5)] is used to compute the flow of electrons from the halide A À to the chalcogen-bond donating molecule D 2 Ch (see ΔQ D2Ch VDD in Table 1).

Chalcogen Bond Strength and Structure
The results of our ZORA-M06/QZ4P calculations are shown in Table 1 for a representative selection of oxygen-, sulfur-, and tellurium-bonded model complexes D 2 Ch···A À , covering D, A = F, Cl, and Br (the complete dataset for all model systems is provided in Tables S1-S2). In the first place, we note that all model reactions are associated with single-well potential energy surfaces (PES), that is, there is no energy barrier separating the reactants from their resulting complex. In the cases where D ¼ 6 A, C s symmetric complexes with D 1 À Ch bond lengths longer than D 2 À Ch and with bond angles Θ 1 ¼ 6 Θ 2 are formed. For the cases where D = A, C 2v symmetric complexes with equal bond distances r ChÀ D 1 = r Ch···A and equal bond angles Θ 1 = Θ 2 are formed (see Table 1). In general, chalcogen bonds become stronger on descending group 16 in the periodic table, in agreement with previous ab initio results. [4b,17] The heavier D 2 Ch···A À chalcogen bonds (i. e., Ch = S, Se, and Te) become weaker and longer as the accepting halide (A À ) varies from F À to Br À . In the case of the telluriumbonded complexes D 2 Te···A À , for example, the ΔE weakens from around À 73 kcal mol À 1 for A À = F À to around À 38 kcal mol À 1 for Table 1. Activation strain analyses (in kcal mol À 1 ) of a representative set of D 2 Ch···A À at the equilibrium geometries (in Å, deg.) [a]  A À = Br À (see Table 1). However, the oxygen chalcogen bonds D 2 O···A À display a more complex dependency of ΔE upon variation of the accepting halide A À . From A À = F À to Cl À , the oxygen-bond strength still weakens, similar to the situation for the heavier chalcogen bonds. However, thereafter, from A À = Cl À to Br À , the oxygen-bond strength does not weaken but instead becomes stronger. This is most clearly seen in the series constituted by the complex F 2 O···A À between an oxygen molecule and a halide ion. Here, ΔE for the oxygen bond strength varies along A À = F À , Cl À , and Br À with values of À 21.9, À 9.9, and À 11.5 kcal mol À 1 , respectively (see Table 1). When the substituent D is varied from D = F to D = Br, the heavier chalcogen-bond strength (i. e. Ch = S, Se, and Te) changes only slightly (see Table 1 and S1 in the Supporting Information). For example, along the series from F 2 Te···F À to Br 2 Te···F À , the tellurium bond strength varies only from a ΔE value of À 72.4 to À 72.0 kcal mol À 1 , the tellurium-bond distance r Ch···A decreases in value from 2.054 to 2.040 Å, and the stretch Δr D 1 À Ch increases in value from 0.162 to 0.227 Å. An exception to this is again the oxygen bond D 2 O···A À , which becomes weaker and longer as D is varied from F to Br (see Table 1). For example, along the series from F 2 O···F À to Br 2 O···F À , the oxygenbond strength weakens from a ΔE value of À 21.9 to only À 12.9 kcal mol À 1 , the oxygen-bond distance r Ch···A increases in value from 1.784 to 2.162 Å, and the stretch Δr D 1 À Ch decreases in value from 0.408 to 0.153 Å.

Bond Analyses with Variation of Ch
The strengthening of chalcogen bonds D 2 Ch···A À , as Ch varies along O, S, Se, and Te, with no change in the donating atom (D) and the accepting halide (A À ), is related to the increasing electronegativity difference across the DÀ Ch bonds as Ch descends in the periodic table, which is translated into two main effects. Firstly, this causes the Ch atom to become increasingly positive along O, S, Se, and Te (see VDD atomic charges in Table 2), resulting in a greater electrostatic attraction. Secondly, this causes, among other effects that will be explained later, the σ* DÀ Ch antibonding 4a' acceptor orbital to have higher amplitude on Ch (see Figure 1), resulting in stronger HOMOÀ LUMO orbital interactions.
The trend in bond strength ΔE is mainly determined by the interaction energy ΔE int . For example, from F 2 O···F À to F 2 Te···F À , ΔE is strengthened from À 21.9 to À 72.4 kcal mol À 1 while ΔE int is strengthened from À 50.2 to À 80.3 kcal mol À 1 (see Table 1). The trend in ΔE is further enhanced by the strain energy (ΔE strain ), which becomes less destabilizing from Ch = O to Te. However, the differences are smaller than the differences in ΔE int . For example, from F 2 O···F À to F 2 Te···F À , ΔE strain is weakened by 20.4 kcal mol À 1 (from 28.3 to 7.9 kcal mol À 1 ; see Table 1), while ΔE int becomes 30.1 kcal mol À 1 more stable.
To understand the origin of these trends, we have carried out activation strain analyses along the entire reaction coordinate ζ, projected onto the stretch in D 1 À Ch bond, Δr ChÀ D 1, that occurs as the chalcogen-bond accepting A À atom approaches the D 2 Ch molecule (see Theoretical Methods section). The resulting activation strain diagrams (ASD) including EDA terms of the interaction are shown for a representative example series, namely, F 2 O···F À to F 2 Te···F À , in Figure 2 (for the complete dataset, see Table S2). Again, the trend in bond energy ΔE is mainly determined by ΔE int (ζ), which strengthens when going from Ch = O to Te (Figure 2, left). On the other hand, the ΔE strain (ζ) curves almost coincide. However, the strain curves reach a final point at ζ eq , that is, the equilibrium geometry of the complex; and here the strain energy ΔE strain (ζ eq ) becomes more destabilizing from Ch = Te to O. Note that the trend in strain energies at the equilibrium geometries along the series of F 2 Ch···F À complexes (see Table 1) arises mainly from changes in the steepness of the interaction curves, not from the relatively minor variation in the strain curves (see Figure 2). Thus, as the F 2 Ch···F À interaction gets weaker along Ch = Te, Se, S and O, the interaction curve becomes shallower and the balance between strain and interaction curve, i. e., the stationary point of the complex, occurs at longer and longer FÀ Ch distances and, consequently, more destabilizing ΔE strain (ζ eq ) (see Table 1).
To understand the trends in ΔE int (ζ), we further decomposed the ΔE int into the individual energy components (Fig-Table 2. Bond lengths (in Å), bond angle (in deg.), VDD charge (in a.u.), orbital energies (in eV) and the homolytic bond dissociation energy without ZPE (in kcal mol À 1 ) of isolated D 2 Ch fragments. [a] D 2 Ch ure 2, right). The electrostatic energy ΔV elstat (ζ) is the least stabilizing for Ch = O and then strengthens along S, Se, and Te. This can be understood by the increasing differences in electronegativity across the DÀ Ch bonds when going from O to Te, resulting in a larger positive charge on Ch. For example, the VDD atomic charge on Ch in F 2 O, F 2 S, F 2 Se, and F 2 Te amounts to + 0.09, + 0.18, + 0.28, and + 0.31 a.u., respectively, and becomes even more positive as the D 1 À Ch bond elongates (see Figure 3a). Nevertheless, our analyses reveal that the chalcogen bonding mechanism is absolutely not purely electrostatic but instead has a relatively large covalent component (ΔE oi ), stemming mainly from the HOMOÀ LUMO interaction between the occupied halide np y atomic orbital (AO) and the σ* DÀ Ch antibonding 4a' acceptor orbital (see Figure 1). The associated charge transfer from A À to D 2 Ch is reflected by the ΔQ D2Ch VDD , which is negative, i. e., D 2 Ch gains charge from A À upon complexation, for all D 2 Ch···A À complexes (see Table 1). For example, ΔQ D2Ch VDD is À 0.37 a.u. for F 2 O···F À and À 0.32 a.u. for F 2 Te···F À . The HOMOÀ LUMO charge transfer nature of the chalcogen bond is also clearly reflected by the associated deformation density. This is illustrated by the 3D plots of the deformation densities associated with chalcogen-bond formation in F 2 S···F À and F 2 Te···F À (see Figure 4). As can be seen, there is charge depletion on the Lewis base F À (and in between the Ch···F À bond due to the Pauli repulsion [10a] ) and charge accumulation on D 2 Ch. Note the 3D shape of the regions of charge depletion and accumulation: they reflect the shape of the 2p-type lone pair from which the F À Lewis base donates and the σ* DÀ Ch antibonding 4a' acceptor orbital on D 2 Ch into which this charge is donated, respectively, in the HOMOÀ LUMO interaction. For the chalcogen bonded complexes, the orbital interaction term ranges from 37 % for F 2 Te···F À to as much as 76 % for Br 2 O···F À of the total bonding interactions (ΔE oi + ΔV elstat ; see Table S2 in the Supporting Information). As can be seen in our energy decomposition diagram, the orbital interaction curves ΔE oi (ζ) become more stabilizing from Ch = O Figure 1. Schematic molecular orbital diagram for a) isolated D 2 Ch fragments at C 2v symmetry (blue: a 1 ; green: a 2 ; red: b 1 ; black: b 2 ) and b) D 2 Ch···A À complexes. The first column in (b) refers to the isolated D 2 Ch fragment and the second column refers to the D 2 Ch fragment deformed to its C s symmetric geometry in the complex (blue: a'; red: a''), in which one DÀ Ch bond has been elongated. See Figure S1 in the Supporting Information for the 3D isosurfaces of the orbitals. to Te (Figure 2, right). The stronger orbital interaction for the heavier chalcogens is the result of the larger LUMOÀ HOMO overlap (i. e. h4a' j np y i; see Figure 1 for the MO diagram which shows the np y orbital of A À pointing towards the D 1 À Ch bond of the D 2 Ch fragment) as Ch becomes more electropositive. For example, in the Cl 2 Ch···Cl À series, h4a' j np y i increases from 0.12 to 0.20 to 0.22 to 0.24 along Ch = O, S, Se, and Te in the equilibrium geometry (see Table S2 in the Supporting Information). The larger percent contribution of the covalent component on oxygen bonds is simply because the electrostatic attraction is relatively weak, caused by the smaller positive charge on O (see Table 2). Whereas ΔE oi (ζ) becomes more stabilizing from Ch = O to Te, it becomes comparable in magnitude for all chalcogens in the equilibrium geometry ΔE oi (ζ eq ). This is a consequence of the fact that the F 1 À Ch bond expansion becomes more pronounced when going from Ch = Te to O. The increasing F 1 À Ch bond expansion causes the σ* DÀ Ch antibonding 4a' acceptor orbital (see Figure 1a) to drop further in energy for lighter chalcogens, resulting in a smaller HOMOÀ LUMO gap and hence more stabilizing donor-acceptor interactions. This effect can be observed in Figure 3a, which shows the energies of the σ* FÀ Ch antibonding 4a' acceptor orbitals along the reaction coordinate. For Ch = S, Se, and Te, the energy of the σ* FÀ Ch antibonding 4a' acceptor orbital converges to an energy value of À 3.8 eV as the chalcogen bond is formed. For Ch = O, on the other hand, the σ* FÀ O antibonding 4a' acceptor orbital energy quickly drops to a value of À 6.4 eV, because the overlap between the F and O AOs is more sensitive to the DÀ Ch distance than for the more diffuse AOs of heavier Ch (see Figure S2 in the supporting information). However, it is counteracted by the orbital overlap between the σ* DÀ O antibonding 4a' acceptor orbital and the np y donor orbital, which is significantly worse for Ch = O than for other chalcogen systems (see Table S2 in the Supporting Information).

Bond Analyses with Variation of A À
Our analyses show that the weakening of heavier chalcogen bonds D 2 Ch···A À (Ch = S, Se, Te), as the accepting group varies from A À = F À to Br À , is directly related to the concomitant reduction in electron-donating capacity of the np-type HOMO and thus the Lewis basicity of the A À halide. [20] We recall the chalcogen bonds display both an electrostatic component (ΔV elstat ) and a covalent component (ΔE oi ). The latter stems mainly from the HOMOÀ LUMO interaction between the occupied halide np atomic orbital (AO) and the σ* DÀ Ch antibonding 4a' acceptor orbital (see Figure 1). Both ΔV elstat and ΔE oi are weakened as the halide HOMO becomes more diffuse and effectively lower in energy from A À = F À to Br À (see Table S2). [20b] Consequently, the interaction energy (ΔE int ) and, thus, the net chalcogen-bond strength ΔE becomes less stabilizing along A À = F À to Br À (see Table 1 and Table S1 in the  Supporting Information). This is very similar to what was found for hydrogen bonds DH···A À and heavier halogen bonds DX···A À (X = Cl, Br, I). [9] The key to understanding why oxygen bonds D 2 O···A À show a more complex, partially opposite trend (i. e., the expected weakening from A À = F À to Cl À but thereafter a strengthening along A À = Cl À to Br À ) is contained in the counteracting effects evolving from DÀ O bond stretching induced in the triatomic D 2 O molecule as it interacts with the halide A À . Interestingly, activation strain analyses reveal, again, that interaction energies recover the original trend in total energies, that is, ΔE int (ζ) weakens from A À = F À to Br À . This can be seen in Figure 5 which shows the activation strain and energy decomposition diagrams along the reaction coordinate ζ projected onto the stretch Δr D 1 À Ch for two representative series. Each diagram in Figure 5 refers to one particular F 2 O or F 2 Te molecule forming a chalcogen bonding with A À = F À , Cl À , and Br À . The ΔE strain curves within each subgraph coincide because they refer to the same DÀ Ch bond in the same triatomic molecule being stretched as the complexation reaction progresses. Consequently, the trend A À = F À to Br À in the total F 2 O···A À and F 2 Te···A À energy profiles ΔE in each subgraph is directly determined by the trend in the corresponding ΔE int curves.
The reason why the oxygen bonds D 2 O···A À do not experience a weakening in ΔE int from A À = F À to Br À , as all other chalcogen bonds, is promoted by a combination of factors: i) a weak DÀ O bond that is easily stretched; ii) a strong interaction with an approaching halide A À ; and iii) a σ* DÀ Ch antibonding 4a' acceptor orbital that drops in energy, more quickly than for other D 1 À Ch bonds due to a more sensitive overlap between the D 1 and O AOs, as the D 1 À O bond elongates (see Figure 3 and Figure S2 in the supporting information). The latter generates a stronger driving force for D 1 À Ch stretching in D 2 Ch···A À because this deformation enhances the orbital interactions and thus ΔE int . Note that, for D 2 O···A À , ΔE oi is the strongest bonding component and that the ΔE oi (ζ) curves directly reflect the electron-donating capacity of the np-type HOMO of the A À halides, that is, the ΔE oi curves become more stabilizing from A À = Br À to F À (see Figure 5). Indeed, D 1 À Ch stretching is most pronounced if this bond in the neutral fragment is weaker, that is, for the weaker chalcogen bonds (e. g., ca. 38 kcal mol À 1 for FÀ O, ca. 35 kcal mol À 1 for ClÀ O and ca. 34 kcal mol À 1 for BrÀ O; see Table 2). In this case, it is able to affect the trend in overall bond strength ΔE. The D 1 À O stretching in oxygen-bonded complexes is most pronounced in the Cl 2 O···A À series, along which the Cl 1 À O stretch Δr D 1 À Ch varies Figure 5. Activation strain (left panel) and energy decomposition (right panel) analyses of a) F 2 O···A À and b) F 2 Te···A À (black, A À = F À ; blue, A À = Cl À ; red, A À = Br À ). Full Papers doi.org/10.1002/open.202000323 between 0.5 and 0.9 Å, but it is already relevant in the F 2 O···A À series in which the F 1 À O stretch Δr ChÀ D 1 varies between 0.4 and 0.6 Å from A À = F À to Br À (see Table 1).

ChemistryOpen
We conclude that, in general, chalcogen bonds D 2 Ch···A À become weaker along A À = F À to Br À because the larger radii and lower np AO energies of the halides lead to weaker electrostatic attraction and weaker orbital interactions. The trend in D 2 O···A À oxygen bond strength is partially inverted, that is, ΔE becomes more stabilizing along A À = Cl À and Br À because of a subtler interplay of factors. Notably, a significant stretching of the relatively weak DÀ O bonds in the D 2 O···A À equilibrium structures lowers the σ* DÀ O antibonding 4a' acceptor orbital and thus amplifies the donor-acceptor orbital interactions.

Bond Analyses with Variation of D
The strength of the heavier chalcogen bonds D 2 Ch···A À varies little when going from D = F to Br because the ClÀ Ch and BrÀ Ch bonds are significantly weaker than the FÀ Ch. This allows the ClÀ Ch and BrÀ Ch bonds to stretch to a higher extent and, therefore, to have more stabilizing electrostatic attraction and orbital interactions. For the oxygen bonds D 2 O···A À , the bond energy is weakened along the same variation because the DÀ Ch bond strength are all comparable (see Table 2). In both cases, the trend in bond strength ΔE is determined by the interaction energy ΔE int . For example, from F 2 O···F À to Br 2 O···F À , ΔE int is weakened from À 50.2 to À 17.4 kcal mol À 1 , respectively, whereas from F 2 Te···F À to Br 2 Te···F À , the bond energy only changes from À 80.3 to À 81.9 kcal mol À 1 , respectively (see Table 1). The strain energy (ΔE strain ) is not negligible, but it does not offset the trend set by ΔE int . Our activation strain analyses explain the above differences between oxygen and heavier chalcogen bonds (see Figure 6).
Starting with some general observations, we find that for oxygen, as well as heavier chalcogen bonds, the ΔE strain curves are most unfavorable when D = F and gradually become less destabilizing as the donating atom is varied along D = F, Cl, and Br (see Figure 6). Furthermore, for all D 2 Ch···A À complexes, the ΔE int curves become less stabilizing along D = F, Cl, and Br. The resulting energy profiles of D 2 Ch···A À depend on the balance between both ΔE strain and ΔE int , but the interaction energy curves already show a very similar trend to ΔE.
The slope and shape of the ΔE strain curves is of course directly related to the D 1 À Ch bond strength of the neutral fragment, which in general becomes stronger as the polarity across the DÀ Ch bond increases [21] (see Table 2). From F 2 Ch to Br 2 Ch, where Ch is S, Se or Te, the halogen-chalcogen bond strength decreases significantly from a value of ca. 93 to 50 kcal mol À 1 ( Table 2). The corresponding halogenÀ oxygen bonds are all much weaker, and variations in the homolytic bond dissociation energy (BDE) are also much smaller. From F 2 O to Br 2 O, the bond strength decreases from 38.0 to 33.7 kcal mol À 1 . Thus, for the heavier chalcogen-bonded complexes, where Ch is S, Se, or Te, the ΔE strain curves show a pronounced reduction in slope from F 2 Ch to Br 2 Ch, which, in the corresponding chalcogen-bonded complexes F 2 Ch···A À to Br 2 Ch···A À , translates into an increasing stretch Δr D 1 À Ch of the neutral fragment. As the stretch Δr D 1 À Ch becomes larger from equilibrium structures F 2 Ch···A À to Br 2 Ch···A À , the ΔE int curves have been able to descend further, to lower, more stabilizing energies. This stabilization is, of course, related to the ΔV elstat and ΔE oi . Note that the electrostatic attraction and orbital interaction curves become less stabilizing along D = F, Cl, and Br, but turn out to have comparable strength in the equilibrium structures, because the D 1 À Ch bonds have been increasingly stretched in the latter, that is, in Cl 2 Ch···A À and Br 2 Ch···A À . The bonding components ΔV elstat and ΔE oi are the most stabilizing for D = F because of the larger difference in electronegativity across the DÀ Ch bonds (vide supra). However, the ΔE oi is able to further stabilize for D = Cl and Br because, in the equilibrium structure of the chalcogen-bonded complexes, the ClÀ Ch and BrÀ Ch bonds expand to a higher extent, resulting in a stronger stabilization of their σ* DÀ Ch antibonding 4a' acceptor orbitals (see Figure 3b). Furthermore, the VDD atomic charge on Ch becomes increasingly more positive as the D 1 À Ch bond expands, which translates into more stabilizing ΔV elstat for D = Cl and Br in the equilibrium geometry. The final result is, thus, a comparable stability among heavier chalcogen bonds D 2 Ch···A À complexes when the substituent D is varied from F to Br.

Chalcogen Bonds Versus Halogen and Hydrogen Bonds
Our analyses highlight that chalcogen bonds, halogen bonds, and hydrogen bonds are all similar in nature. [9] Each of these bonds in our set of model systems has a significant covalent component in addition to electrostatic attraction, and can range in strength roughly between À 6 and À 70 kcal mol À 1 (see Figure 7). Chalcogen bonds and halogen bonds have a larger range in polarities in DÀ Ch and DÀ X than in DÀ H bonds and are in general stronger than hydrogen bonds because of more stabilizing orbital interactions (see Table S3 for bond energies ΔE of a representative series of XB and HB ΔE computed at ZORA-M06/QZ4P). However, chalcogen bonds and halogen bonds also have more destabilizing Pauli repulsion because the lone-pair HOMO of the Lewis base overlaps with more closed shells, in particular, with the σ DÀ Ch bonding 3a' and 2a' FMOs or σ DÀ X bonding FMO with a higher amplitude on Ch and X, respectively, than the amplitude of σ DÀ H bonding FMO has on H (see Figure 7; see also Ref. 9). Our analyses provide a unified picture for chalcogen bonds, halogen bonds, and hydrogen bonds based on quantitative Kohn-Sham molecular orbital theory, which proves that these intermolecular interactions cannot be described by a pure and simple electrostatic model.

Conclusions
Chalcogen bonds in D 2 Ch···A À range between 6 and 73 kcal mol À 1 in strength, becoming stronger as the chalcogen atom becomes more electropositive, along Ch = O, S, Se and Te, and also as the halide becomes a stronger Lewis base, along A À = Br À , Cl À and F À . The trend upon variation of the substituent along D = F, Cl, Br is less pronounced, as are all trends for the relatively weak oxygen bonds. This follows from our bonding analyses based on relativistic density functional theory (DFT) calculations at ZORA-M06/QZ4P.
Our activation strain and quantitative Kohn-Sham MO bonding analyses reveal that the chalcogen bonds in D 2 Ch···A À are similar in nature to halogen bonds in DX···A À and hydrogen bonds in DH···A À (Ch = O, S, Se, Te; D, X, A = F, Cl, Br). Chalcogen bonds are far from being solely electrostatic phenomena. Similar to halogen and hydrogen bonds, chalcogen bonds have a sizeable covalent component, ranging up to 80 % of the bonding components (ΔV elstat + ΔE oi ), stemming from HOMOÀ LUMO interactions between the np-type lone pair on the bond accepting fragment A À and the LUMO with strong DÀ Ch σ* anti-bonding character on the bond donating fragment D 2 Ch.
Chalcogen bonds become stronger for heavier Ch because of the greater difference in electronegativity across the DÀ Ch bonds, causing: i) the σ* DÀ Ch antibonding 4a' acceptor orbital to have higher amplitude on Ch, enhancing HOMOÀ LUMO orbital interactions; and ii) the Ch to become more positively charged, resulting in greater electrostatic attraction when descending in group 16 of the periodic table. The chalcogen bonds also become stronger for lighter A À because the electron-donating capacity of the np-type HOMO (i. e. Lewis basicity) of the halides increases ascending group 17 in the periodic table. The trends for oxygen bonds, as well as along various D, are less pronounced because of counteracting effects or small values in bond strength.