We begin with a look at the shapes of electron density—a simple but crude instrument. We focus on density around a Lewis acidic site. Figure 1 shows the isodensity plots for a number of molecules discussed in this work. Going from the left to right, the density contours show the deep depletion in BeO, still visible depletion in AuF, a slight depletion in ClF, especially in a short region of Cl, an axial flattening of density in I-F around the iodine, convex contours in SCO around sulfur with some axial flattening, and finally the convex contours in N₂ with no flattening. The axial flattening, also called polar flattening in literature (Sedlak et al., 2015), does not occur in N₂; rather, its density protrudes axially in the lone pair region. In fact, the quadrupole moment of N₂ is negative, thus indicating that the negative charge extends further along the molecular axis.

Now, how is this picture reflected in the electrostatic potential mapped on the threshold density?

Three limited cases are presented in Figure 2: a dipolar and two quadrupolar cases. It shows the contours of the electrostatic potential that emanate from a dipole and two quadrupoles of opposite signs. One may conclude that the positive “σ hole” in the left two and the negative “σ bump” offer no unique insights except for pointing to the sign of the leading multipole moment of the molecule. Arguably, when the density surface shows an axial indentation, as in AuF in Figure 2, the magnitude of electrostatic potential at this point is enhanced. Nevertheless, the insights from this picture are qualitative at best.

Recent years have seen an explosive growth (Web of Science Core Collection, 2022) of literature on the features of electrostatic potential painted on the density (Cavallo et al., 2016; Kolář and Hobza, 2016). The types of holes that emerged from these studies range from σ to π to δ (Angarov and Kozuch, 2018). Simultaneously, the names of bonds associated with these “holes” exploded up to 8 by the last count (Zierkiewicz et al., 2021). The “hole” properties such as the spatial extent of the positive region (Kolář et al., 2015), the value of the maximal electrostatic potential, V_smax, at a certain threshold density (Politzer et al., 2013), and the extent of the hole’s “polar flattening” (Sedlak et al., 2015) have been correlated with the bonding strength of a complex.Simultaneously, there were claims that only the density and its associated electrostatics represent “phenomena” and all the other computable and physically interpretable effects represent “unicorns” (Politzer et al., 2015). This literature avalanche gives an impression that all of the donor–acceptor interactions are electrostatically controlled, which is incomplete if not misleading. There are some dissenting voices from this chorus stressing the importance of the exchange effect (Adhikari and Scheiner, 2012; Stone, 2013), the ground-state charge-transfer interactions (Wang et al., 2014; Rezac and de la Lande, 2017; Thirman et al., 2018), orbital interactions (Pascoe et al., 2017), or covalency contributions (Bora et al., 2017).

The electrostatic effect and “σ holes”, which are expected to explain the directionality of these interactions, are purely classical concepts. To fully deconstruct the properties of the “hole” and to detect the existence of density depletions in the direction of the “hole”, we need a more sophisticated tool—a quantum effect of exchange. This requires the participation of at least a two-electron interacting partner.

Figure 3 shows three complexes, He-BeO, He-AuF, and He-IF. In the first complex, the density hole is evident and so deep that it can accommodate a He atom but not Ne. In the second one, the density still has some depletion but to a much lesser degree. In the third one, only a flattening is seen (Figure 1). We will show that regardless of the depth of the density shape, the same pattern persists if the energy components are considered. The systems are also chosen by steeply decreasing dipole moments from 6.6D (BeO) to 4.5D (AuF) to 2.7D (IF) (as per CCSD computation). Despite these widely varying electrical properties, the nature of bending potentials remains the same. It is the exchange repulsion term that displays a deep “hole” for the linear configuration and thus dictates the shape of the overall interaction energy (shown as SAPT2+ in Figure 3). The electrostatic term favors nonlinear geometries and disfavors the linear ones. It is of no importance since it is of the overlap character. The induction term’s importance reflects the size of the dipole moment: the effect is maximized in He-BeO and minimized in He-IF. The dispersion energy is an important attractive contributor but its angular variations are minimal.

The same behavior could be expected for the radial coordinate: Namely, the depletion of the exchange repulsion allows the subsystems to approach more closely and benefit from enhanced dispersion and induction effects.

The cavity in the interaction potential of He-BeO was first demonstrated in the work of Hapka et al. (2013) in the form of zero-interaction energy boundary showing a deep indentation for the linear geometry on the Be side. They estimated the binding energy of nearly 1,900 cm⁻¹ or 5.4 kcal/mol. Ne is similarly bound to BeO but penetrates less into the cavity by 0.267 Å according to Nunzi et al. (2017). The Ar-BeO study estimated the well depth of 3,990 cm⁻¹ or 11.4 kcal/mol (Tebai et al., 2014), quite unusual for the rare gas binding to a molecule. The SAPT conclusion is that the minimum in He-BeO is largely due to the exchange depletion reinforced by the induction effect. The work of (Nunzi et al., 2017) goes even further by postulating that the unique strength of He-BeO as compared to Ne-BeO, for example, is due to the ability of BeO to back-donate to the vacant 2p orbitals of He—and the inability to do so to Ne because its 2p orbitals are already occupied. Our observation is that (1) Ne is simply too large to penetrate as deep into the density depletion in BeO as He could and (2) the induction effect is actually stronger in Ne than in He.

We examine this intriguing hypothesis of back-transfer below (Figure 4). First, one expects that if there is an overlap between the occupied π MOs of BeO with vacant 2p orbitals of He, it should favor the linear configurations. Second, it should manifest itself in the Molecule → Atom induction contribution. These contributions (along with a residual δ_HF) are shown in Figure 4 for both He and Ne complexes with BeO taken at the same distance R = 3.0 Å, which is the equilibrium for Ne-BeO.

It appears that BeO → He/Ne contributions are strongly attractive, favoring the linear structure. Their respective exchange counterparts cancel very small fractions of these attractions.

The magnitudes of the induction terms of the BeO → He/Ne type are consistent with ca. 2:1 ratio of static dipole polarizabilities of Ne to He. In both cases, the BeO → He/Ne induction attraction strongly favors the linear structure. What if the donation to virtual 2p orbitals shows up in the residual term δ_HF? Here, in neither case is δ_HF the most attractive for the linear geometries. Thus, in the SAPT language, the only indication to support the BeO-to-He donation is in very strong BeO → He induction. One observation, however, is worth emphasizing: Namely, in the Ne-BeO complex, we can see very strong, common sense eluding Ne → BeO contribution of ca. 2 kcal/mol at skewed geometries. Thankfully, the exchange counterpart of this term varies as almost its mirror image and nearly completely wipes out the common-sense-defying effect. It is useful to interdict to explain why we use the term “common-sense-defying.” This is because the induction effect allows electrons of one monomer to occupy the already occupied orbitals of the partner, thus violating the Pauli principle, if unconstrained by the exchange. The lesson from this is that employing the induction effect alone is error prone. The induction energy must be accompanied by its exchange counterpart, if unphysical effects are to be avoided.

To determine the effect of the nonsphericity of halogen on the shape of the interaction potential, we examine a series of electron donors (Lewis bases) with one Lewis acid ClF, the latter showing the prominent “σ hole,” to establish in what circumstances is electrostatic effect the controlling one. This question was previously explored in the literature (Stone, 2013), but we intend to place it in a wider context.

In Figure 5, four Lewis bases are shown interacting with ClF, from very weak to ionic, as functions of the bending angle of ClF around the Cl atom (denoted φ) from linear to perpendicular. It is noteworthy that in neither case is the induction effect controlling the shape of the bending potential and neither is the dispersion energy nor the electrostatic. The Ne-ClF complex’s (Figure 5A) total SAPT2+ interaction energy is the most binding for the linear geometry owing to the exchange repulsion’s minimum of ca. 0.5 kcal/mol compared to the perpendicular arrangement. Indeed, this is the true manifestation of the hole. The electrostatic effect is overlap-dependent and varies in the opposite way. Two σ donors described in what follows, N₂ and NH₃, are considered at the same intersystem distance N---Cl of 2.5 Å, which equals the equilibrium separation in the latter. In the BF₄ ⁻-ClF complex, the same intersystem distance 2.5 Å is employed because it is close to the complex’s equilibrium separation. In a non-polar N₂ interacting with ClF, the exchange hole amounts to 15 kcal/mol, dominating the total SAPT2+ energy. In fact, one curve appears to be the vertical translation of the other. The only other term that favors the linear geometry is δ_HF (see below for more discussion). In the polar NH₃ interacting with ClF, considered at the same intersystem distance, the exchange hole amounts now to 18 kcal/mol. The exchange effect again controls the shape of the bending potential as well as the well depth of this complex. Surprisingly, given the participation of the two polar subsystems, the electrostatic energy is almost constant with angular variation between −19.2 kcal/mol at 0° and −19.0 kcal/mol at 90°. Similarly, the induction effect (more specifically its perturbation theory part, E_Ind-PT) is also nearly constant (to 0.6 kcal/mol). To emphasize: the electrostatic and induction effects are isotropic, or nearly so, around the Cl end of the molecule. This appears to contradict the electrostatic predominance due to “σ hole” and, in particular, attributing any predictive capabilities to the value of V_smax. The second most anisotropic term is the residual δ_HF contributions to which we will return below.

Finally, in the ion-dipole complex BF₄ ⁻-ClF, one would expect that the electrostatic energy should drive the interaction. However, it is still the exchange hole of over 9 kcal/mol that shapes the potential curve. The electrostatic energy varies only by ca. 3 kcal/mol between 0° and 90°. All the other terms remain constant or nearly so (cf. dispersion). Naturally, at longer separations, one would expect the electrostatics to dominate; this is because the electrostatic effect varies more slowly with the intersystem distance than the exchange effect.

To summarize, with the variety of electron donors, ranging from a rare gas to quadrupolar to dipolar to ionic interacting with ClF, the bending potential around the equilibrium separation is controlled by the exchange hole. The effects of the electrostatic “σ hole” of ClF manifesting themselves in the electrostatic and induction interactions are only secondary.

The complexes considered in the section When is the Electrostatic Effect Important? above were halogen-bonded, which are believed to be stabilized also by the charge-transfer contribution, as predicted by the Mulliken’s interpretation of the complexes of this type (Mulliken and Person, 1969) (before they became known as halogen-bonded); although the subject remains controversial. There have been several schools of thought on inclusion of charge-transfer energies in the interaction energies beginning with Morokuma partitioning of the HF interaction energy in the early 70s (Morokuma, 1971; Kitaura and Morokuma, 1976), the Stone–Misquitta monomer-to-dimer basis set extension of the induction effect in SAPT (Stone and Misquitta, 2009) [see also recommendation in Psi4 (Parrish et al., 2017)], Misquitta’s attempt to “regularize” one-electron V operator terms (Misquitta, 2013), constraining intermonomer DFT and then relaxing it (Wu and Van Voorhis, 2005; Lao and Herbert, 2016), and the energy decomposition scheme ALMO-EDA in its latest variant (Thirman and Head-Gordon, 2015), just to name a few. Finally, there are also views that the charge transfer is fully included in the induction effect (Politzer et al., 2013; Sedlak et al., 2015).

In the halogen bond context, the charge-transfer contribution is understood as energy lowering resulting from the orbital interaction involving the lone pair of Lewis base (n) with the LUMO of halogen molecule σ ∗ of the type (Thirman et al., 2018) E C T = 〈 n | H ^ | σ ∗ 〉 E n − E σ ∗ (2)

This interaction favors the direction where n and σ ∗ overlap best, i.e., the linear geometry. It also benefits from the denominator’s energy rising of n and lowering σ ∗ resulting from the interaction. Since perturbation theory is based on the product of unperturbed monomers’ wave functions, the latter effect will be absent. If one performs the full dimer HF energy calculations, the orbitals that are still located largely on the respective monomers will be relaxed. Therefore, the term δ H F , which is obtained as δ H F =   Δ E H F − E e s ( 10 ) − E e x c h ( 10 ) − E i n d , r ( 20 ) − E e x − i n d , r ( 20 ) (3) where Δ E H F stands for the supermolecular interaction energy at the HF level, should capture such effects. A conventional view of δ H F is that it sums up induction and exchange-induction terms of orders higher than the second. However, it also includes other residual effects including a possible divergence of SAPT due to the weak antisymmetry forcing.

In Figure 5, we purposely showed perturbation-theory-derived induction effects, E_ind-PT, as separate from δ H F to test the hypothesis of δ H F including the charge-transfer effects.

In Ne-ClF, δ H F (Figure 5A) is near zero for the linear geometry and more attractive for the bent geometries. Upon closer inspection, we may notice that in the bent geometries, E_ind-PT becomes positive—an unphysical behavior. Thus, δ H F in this case corrects the SAPT’s convergence. In BF₄ ⁻-ClF (Figure 5D), δ H F is constant to a fraction of kcal/mol while E_ind-PT varies by 1 kcal/mol in favor of the linear geometry. Therefore, δ H F is expected to be of induction origin. We should contrast these two cases with what happens in NH₃-ClF. In this instance, when we sum up all the perturbational SAPT terms, the well depth is too shallow by half and δ H F provides the missing half. δ H F is strongly attractive for the linear structure and vanishes for the perpendicular one. This is the typical dependence that one would expect from the charge-transfer effect from n to σ*. In N₂-ClF, considered at the same distance, δ H F binding exceeds that from E_ind-PT and again is the most attractive at the linear geometry. However, it does not vanish at the perpendicular geometry indicating the presence of some residual effects in strongly repulsive part of the interaction SAPT2+ potential.

To conclude, δ H F may contain, depending on the type of a complex, a variety of effects: residual induction and exchange in a weak one with a rare gas, Ne-ClF, or strong but with similarly nearly spherical donor BF₄ ⁻-ClF. In the paradigm “charge-transfer”, halogen-bonded complex, NH₃-ClF, one-half of the binding energy in the minimum resides in δ H F , i.e., the term whose angular dependence is consistent with the orbital interaction between n and σ*. Thus, we posit that when charge-transfer energy contributions are warranted (e.g., by the right symmetry and energetical balance of donor–acceptor properties), they would manifest themselves in SAPT in its residual term δ H F .

Finally, we should mention that the use δ H F is only compatible with SAPT based on HF subsystem zero-order wave functions because only in this case, as proved by Moszynski et al. (1996), do the induction terms summed up to the infinite order in V converge to the HF interaction energy (sans the first-order energy). As observed before (Modrzejewski et al., 2012), if SAPT is based on the DFT monomer description, the inclusion of δ H F , which is common practice in the literature to improve results, is problematic.

A new and unique bond has been attributed to the 16-group containing atoms O, S, Se, etc., the chalcogen bond (Pascoe et al., 2017; Wang et al., 2009). In this section, we examine an example chalcogen-bond donor SF₂. The electrostatic potential of SF₂ as shown in Angarov and Kozuch (2018) features two positive regions on both sides of the C₂ axis roughly perpendicular to the S–F bonds. These positive regions are supposed to attract basic regions of a Lewis base. To determine what is really behind this attraction, we show in Figure 6 the SAPT results for the complex SF₂-He. The Jacobi coordinate R = 3.5Å corresponds to the sum of van der Waals radii of He and S and also to the radial minimum in the SAPT2+ potential surface. The angular potential (here derived at the SAPT0 level of theory) is shallow with a weak minimum at theta = 50°. The minimum closely follows the angular variation of the exchange repulsion. Both the electrostatic and induction effects contribute minimally, especially in the minimum region. The only attractive term of some importance is the dispersion energy albeit the least attractive around the theta = 50° angle. Our two-electron probe thus indicates that the controlling effect is the exchange depletion—an exchange hole. The electrostatics and induction are unimportant. Would the same principles apply in the case of nonspherical Lewis base? The work of Angarov and Kozuch (2018) shows that when the Lewis base is a σ-donor HCN molecule, the angle of approach is also tilted ca. 53° away from the C₂ axis just as in the SF₂-He discussed above. Their results show only the total interaction. In a quick look at the actual balance of the energy components, we show in Figure 7 the interaction with another σ-donor N₂: N₂-SF₂. The SAPT2+ curve again follows closely the exchange depletion at the tilted geometry with theta = 50°, just as in SF₂-He. The electrostatic energy favors the axial geometry (theta = 0°) and only coincidentally matches the value of total interaction energy around the minimum. The dispersion energy shows a similar pattern to SF₂-He: the most attractive at theta = 0 and the least attractive in the tilted geometry.

What these two cases have in common is that the Lewis base approaches the C₂ axis of SF₂ at a tilted angle and this happens regardless of the polarity state of the base. This may create an impression that the base is electrostatically attracted to the “σ hole” on the side of the sulfur atom (or, more accurately, above the S–F bond). In reality, this position is preferable due to the depletion in the exchange repulsion. Could this be viewed as a unique chalcogen bond? To answer this question let us recall a 30-year-old perturbation theory result of another chalcogen: H₂O-Ar (Chałasiński et al., 1991). In this complex, Ar is also tilted with respect to the C₂ axis of H₂O, the fact that was confirmed experimentally by the high-resolution vibration–rotation–tunneling spectroscopy (Cohen and Saykally, 1993). The perturbation-theory paper (Chalasinski et al., 1991) states that the angle is determined by the depletion in the exchange repulsion. Incidentally, the electrostatic potential of H₂O does not feature a sigma hole in this region. The σ holes, if we can call them such, reside on H atoms!

The final question concerns the cause of the exchange depletion in the systems studied here and in many other cases discussed in the vast literature of novel non-covalent bonds. The electron density depletion shown in Figure 1 may be obvious in some molecules, such as BeO, but even in the halogen bond donors, the density depletions are not obvious from the density contours. The first-order exchange repulsion detects regions that are, to a varying degree, Pauli-forbidden and thus serve as a very sensitive probe of the doubly occupied space. The regions that are not occupied invariably involve nodes and nodal planes. It is, for example, axiomatic and taught in undergraduate courses that the nodal planes help distinguish bonding from antibonding orbitals or that the T_2g symmetry positions in the octahedral symmetry are more ligand-favorable than E_g because they lie on the intersections on the d_π orbitals of the transition metal. The same principle applies to the exchange holes of the systems studied here (and elsewhere). A good example is the case of SF₂ shown in Figure 8. On the left, there is the electrostatic potential mapped on the density, and on the right, the highest occupied molecular orbital (HOMO) of the molecule. The depletions seen above in the exchange repulsion (Figure 7) occur on the intersection of the nodal plane of the π orbital in plane of the molecule and another node that runs perpendicular to the S–F bond.

All Lewis acids examined in this work have the same features: the regions of depleted density occur on the intersections of nodal planes (the illustrations are included in the Supplementary Material). One easily explainable exception to this rule among the present systems involves AuF (Supplementary Figure S2). In AuF, as seen in Figure 3, the exchange depletion occurs in the axial position on the Au side whereas the HOMO orbital is d_σ along the molecular axis. However, due to the relativistic contraction of the 6s orbital in Au, the axially extending part of s + d_z2 hybrid is minimized. Close-lying d_π and d_δ occupied orbitals all feature one and two nodal planes, respectively, on the Au side of the molecule.

Molecular orbitals are, after all, one-electron functions and, to an extent, noumena. However, their symmetries have a deeper meaning as the irreducible representations of the nuclear point group. Our point here is that the depletion of the density on nodal planes has consequences to bonding whether covalent or noncovalent, indeed.

The circumstances of an actual electron vacancy as in a Cl radical, often mentioned as a protagonist of the halogen bond (Politzer et al., 2010; Kolář and Hobza, 2016), are entirely different. A ²P state of Cl is triply degenerate and it takes an interaction with another system to remove this degeneracy. The Pauli repulsion along the singly occupied orbital with a closed-shell system will be smaller than along the doubly occupied one giving rise to multiple, coupled potential energy surfaces (Klos et al., 2004). If that system is a positive point charge, as in the case of the electrostatic potential, the interaction in the direction of singly occupied orbital will also be more repulsive as shown in (Politzer et al., 2010) for electrostatic reasons. The open-shell case brings us back to the original meaning of the electronic hole (Koch et al., 1987), which, by and large, is not governed by the nodal surfaces but by an actual electron vacancy. A transplantation of this concept to the closed-shell case of halogen bonding has been a gross simplification.