What Is the Nature of Supramolecular Bonding? Comprehensive NBO/NRT Picture of Halogen and Pnicogen Bonding in RPH2···IF/FI Complexes (R = CH3, OH, CF3, CN, NO2).

We employ a variety of natural bond orbital (NBO) and natural resonance theory (NRT) tools to comprehensively investigate the nature of halogen and pnicogen bonding interactions in RPH2···IF/FI binary complexes (R = CH3, OH, CF3, CN, and NO2) and the tuning effects of R-substituents. Though such interactions are commonly attributed to "sigma-hole"-type electrostatic effects, we show that they exhibit profound similarities and analogies to the resonance-type 3-center, 4-electron (3c/4e) donor-acceptor interactions of hydrogen bonding, where classical-type "electrostatics" are known to play only a secondary modulating role. The general 3c/4e resonance perspective corresponds to a continuous range of interatomic A···B bond orders (bAB), spanning both the stronger "covalent" interactions of the molecular domain (say, bAB ≥ ½) and the weaker interactions (bAB ˂ ½, often misleadingly termed "noncovalent") that underlie supramolecular complexation phenomena. We show how a unified NBO/NRT-based description of hydrogen, halogen, pnicogen, and related bonding yields an improved predictive utility and intuitive understanding of empirical trends in binding energies, structural geometry, and other measurable properties that are expected to be manifested in all such supramolecular interaction phenomena.

1. Introduction

Recent computational [1,2,3,4] and experimental [5,6,7,8] studies have called attention to the multiplicity of halogen [9,10,11,12,13], chalcogen [14,15], pnicogen [16], and tetrel bonds [17,18,19] (corresponding to Groups 17, 16, 15, and 14, respectively) that appear necessary to supplement the familiar hydrogen bonds of intermolecular interaction theory [20]. Such proliferating “bond” forms naturally reignite age-old controversies concerning the relative importance of “covalent” versus “electrostatic” contributions in the general theory of intermolecular forces [21,22,23,24]. Of course, all protagonists concur that the only relevant potential terms in the non-relativistic Hamiltonian are of Coulombic form (hence, tautologically “electrostatic” in nature). However, the deeper issue is whether supramolecular complexation (like molecule formation) is primarily governed by exchange-type (exponential) interactions of the short-range quantal regime or by the exchange-free (power–law) behavior of the long-range classical limit. In short, we ask whether orbital-based (quantal) or “dipole–dipole”-type (classical) conceptions provide the more useful starting point for describing intermolecular forces. For present purposes, the former is exemplified by a natural bond orbital [25] (NBO)-based description, and the latter by is exemplified by symmetry adapted perturbation theory (SAPT) [26,27], electrostatic potential models (e.g., of CHELP [28,29] or “sigma hole” [30,31,32] types), or related classical force field [33,34] models.

In the case of hydrogen bonding, the primacy of short-range covalency forces is now widely recognized [20,35,36,37,38]. In the NBO framework [39], a general B:···H—A hydrogen bond is identified with the strength of a nB→σ*AH donor–acceptor interaction between the lone pair (nB) of the Lewis base B and the valence antibond (σ*AH) of the Lewis acid A—H. In natural resonance theory (NRT) [40,41,42,43] terms, such a nB→σ*AH interaction represents resonance mixing between the parent natural Lewis structure (NLS) I and the secondary “charge–transfer” resonance structure II (with respective NRT weightings w and w): B:···H—A ↔ B(+)—H···:A(−) in which the “covalent” (—) and “hydrogen” bond (···) symbols are interchanged [44]. As recognized by Coulson [45], the resonance hybrid (1a) is also the resonance-theoretic equivalent of the Pimentel–Rundle 3-center, 4-electron (3c/4e) molecular obital (MO) formulation of hypervalency [46,47]. In the simple two-resonance limit, the formal bAH, bBH bond orders are necessarily related by the characteristic “bond conservation” relationship [48], which merely reflects the resonance-averaging condition (w + w = 1) for the composite resonance hybrid. Even in the absence of NRT bond order descriptors, the importance of NBO nB→σ*AH delocalization can be assessed by 2nd-order perturbative energy (ΔE(2)nσ*), intermolecular charge-transfer (Qnσ*), or the deletion of this NBO interaction (ΔE$DELnσ*), as well as the variational recalculation of structure, energetics, and vibrational frequencies [49]. In contrast to the superficial “on-off” picture of H-bonding that is sometimes suggested by graphical software or inflexible analysis methods, the NBO/NRT descriptors all reflect the continuously variable character of resonance-type bonding, as quantified by fractional bond orders ranging over integer and sub-integer values. All such characteristic NBO/NRT signatures are found in all known H-bonded species, including the paradoxical “anti-electrostatic” H-bonds between like-charged ions [50] that most clearly demonstrate the relatively secondary influence of “electrostatic” modulation on the authentic H-bonding phenomenon.

In the present work, we take up the analogous covalency versus electrostatics questions for title complexes RPH2···IF/FI (R = CH3, OH, CF3, CN, NO2) that exhibit a wide range of binding strengths. A smaller subset of such species was previously examined [51] with an NBO/NRT analysis of MP2/aug-cc-pVTZ level computations. The present study extends the choice of species and the variety of NBO/NRT analysis tools to specifically address the broader questions of covalency versus electrostatics in halogen and pnicogen bonding, as well as the general relationships to H-bonding and other types of “X-ogen bonding” that signal a grand unified picture of all such donor-acceptor phenomena. In each case, we exhibit the mutually consistent correlations between NBO/NRT-based descriptors and experimentally measurable properties that illustrate the interpretative and predictive advantages of the unified NBO/NRT conceptual framework compared to alternative analysis methodologies [25].

Note that we ignore the weaker effects of the London dispersion interaction [52], a pure correlation effect that (unlike the n-σ* interaction) is strictly absent in uncorrelated Hartree–Fock theory. Though an empirical dispersion correction [53] is often found to improve the performance of density functional theory (DFT) calculations, such correction is only an additive term in the total energy expression, with no direct effect on Kohn–Sham orbitals, electron density, or other properties. Effects of this correction are therefore “invisible” to NBO/NRT analysis at any specific geometry, but, in principle, they can lead to shifts on the potential energy surface at which an analysis is performed. We evaluated the B3LYP-D3 properties for CH3PH2···IF (see Supplementary Materials for full details) and found that the changes are relatively negligible in magnitude (smaller, e.g., than variations among alternative basis levels discussed below) as well as opposite in sign to that expected for the authentic London dispersion attraction (acting instead to reduce ΔEbind). Such “corrections” are not considered further in the present work.

In the following, we first describe the species, computational levels, and NBO/NRT methods to be employed. The latter include NBO ΔE(2)nσ*, charge–transfer Qnσ*, and $DEL-deletion ΔE$DELnσ* measures of donor–acceptor attraction [49], the corresponding pairwise steric-exchange opposition (ΔE(pw)nσ) of filled orbitals [54] that completes the 3c/4e description of overall bonding energetics, and the various NRT-based regularities and correlations with measurable properties that are considered diagnostic of H-bonding [22]. We then document the systematic application of these NBO/NRT methods to the title species to exhibit their mutual coherence and consistency with well-established features of H-bonding. The concluding summary emphasizes how a balanced NBO/NRT description of leading donor–acceptor interactions can lead to a unified conceptual picture of supramolecular bonding that encompasses the entire range of chemically significant (“non-innocent”) complexation and ligation phenomena.

2. Computational Methods and Results

Our studies began with a larger data set of dihalogen complexes RPH2···XY that included all possible (X, Y = F, Cl, Br, I) dihalogen species of the first four periodic rows. However, from a qualitative conceptual perspective, it soon became apparent that the properties of such complexes vary in a smooth and chemically reasonable fashion with dihalogen polarity, as maximized in the iodine–fluorine (IF) species. Accordingly, we focus here on the “polar extremes” of RPH2···IF versus RPH2···FI complexation, corresponding to opposed signs of any envisioned “dipole–dipole” contributions to intermolecular binding. A direct comparison of these two extremes thereby allows one to recognize the important modulating influence of classical dipole–dipole forces (i.e., in strengthening or weakening the underlying quantal interactions), while also verifying that resonance-type covalency forces yield robust supramolecular bonding even when the presumed classical electrostatic prerequisites for such complexation are profoundly violated. While the computational results presented here focus only on the extremal RPH2···IF versus RPH2···FI limits of this broader picture, the Supplementary Information (SI) includes results for all remaining dihalogen species to illustrate the generality of the resonance-type 3c/4e interaction picture to be sketched below.

Though the conceptual picture we describe is insensitive to many details of the chosen computational method and basis set [55], the inclusion of iodine (for which common Pople-type 6-311++G** basis sets are unavailable) mandates the adoption of a relativistically corrected basis set (e.g., of the pseudopotential LANL2DZ type). For present qualitative purposes, we therefore adopt an unconventional but simple “mixed-PP” basis set of 6-311++G** for all RPH2 monomers and LANL2DZ for dihalogen monomers of each complex. We employ simple B3LYP density functional methodology throughout this work, but alternative MP2 evaluations at different basis levels were also performed for comparison (see SI). All structures were optimized, and frequencies were computed to confirm equilibrium geometries corresponding to minimum energy points. The binding energy, ΔEbind = E(AB) − E(A) − E(B), was calculated as the direct energy difference between optimized dimer and relaxed monomers. All electronic structure calculations were performed with the Gaussian 16 program package [56]. NBO/NRT analyses were performed with NBO 7.0 software [57], and both NBO/NRT analyses and orbital imagery were obtained with the NBOPro@Jmol program [58].

Table 1 summarizes the key structural and binding energy parameters of RPH2···IF and RPH2···FI species as obtained at the adopted “B3LYP/mixed-PP” level, showing the wide range of tuning by various R-substituents. Despite the fact that all complexes are appreciably bound (by circa 3–15 kcal/mol, in defiance of superficial “dipole–dipole” or Coulombic point–charge expectations), interesting structural contrasts are immediately evident between these two classes of species. Focusing first on the orientation angles (Θ) in Table 1, one can see that the RPH2···IF species all maintain nearly linear P···IF alignment (ΘPIF ≈ 180°) and a circa 120° bend angle between RP and IF axes (ΘRPI ≈ 120°). In contrast, the RPH2···FI species are noticeably tilted away from P···FI linearity (ΘPFI ≈ 160°), instead adopting increasingly near-linear alignment with respect to RP and FI axes (ΘRPF ≈ 170°). Such structural tendencies suggest that the principal underlying “linearizing” (H-bond-like) interaction differs in the two cases and that other RPH2···XY dihalogen species of intermediate polarity will likewise exhibit intermediate orientational preferences between these two competing tendencies. Indeed, this is the broader picture obtained from the full set (X, Y = F, Cl, Br, I) of optimized dihalogen structures, with the homopolar RPH2···X2 species forming the approximate “dividing line” between RPH2···IF-like versus RPH2···FI-like structural propensity. These structural features of RPH2···IF-like (“P···IF”) versus RPH2···FI-like (“RP···F”) complexes are also consistently exhibited at other theory levels, as shown in SI Tables S1–S3. All such comparisons indicate the adequacy of the adopted B3LYP/mixed-PP model and the extremal RPH2···IF versus RPH2···FI species for describing the full range of RPH2···XY complexation (X, Y = F, Cl, Br, I; R = CH3, OH, CF3, CN, NO2), including general consistency with previous MP2/aug-cc-pVTZ-level calculations.

Table 1

Optimized dihalogen RPH2···IF and RPH2···FI complexes (B3LYP/mixed-PP level), showing key structural and energetic descriptors: Bond lengths (RP···X, RIF, RRP; Å) and changes from monomer values (ΔRIF, ΔRRP; Å), orientation angles (ΘRPX, ΘPXY; degrees), dihalogen infrared frequency shift (ΔνIF; cm−1), and binding energy (ΔEbind; kcal mol−1).

RPH2···IF	R P···I	R IF	ΔRIF	R RP	ΔRRP	ΘRPI	ΘPIF	ΔνIF	ΔEbind
CH3PH2···IF	2.859	2.071	0.064	1.845	−0.027	115.4	178.2	−84.4	14.60
OHPH2···IF	2.854	2.066	0.059	1.652	−0.036	117.1	179.2	−79.6	12.36
CF3PH2···IF	2.992	2.038	0.031	1.894	−0.002	122.1	177.9	−50.1	7.18
CNPH2···IF	3.031	2.034	0.027	1.783	−0.014	121.2	177.6	−45.0	5.93
NO2PH2···IF	3.019	2.031	0.024	1.880	0.000	123.1	175.3	−42.8	5.10
RPH2···FI	RP···F	RIF	ΔRIF	RRP	ΔRRP	ΘRPF	ΘPFI	ΔνIF	ΔEbind
CH3PH2···FI	2.323	2.114	0.107	1.873	0.001	166.1	164.4	−211.9	3.25
OHPH2···FI	2.075	2.197	0.190	1.684	−0.004	164.8	157.4	−320.6	5.85
CF3PH2···FI	2.624	2.025	0.018	1.900	0.004	172.7	165.1	−52.6	4.43
CNPH2···FI	2.653	2.020	0.013	1.805	0.008	173.1	164.6	−37.3	4.85
NO2PH2···FI	2.339	2.055	0.048	1.902	0.022	177.1	160.0	−116.2	6.79

Figure 1 depicts the numerical structural data of Table 1 in graphical form, showing the distinct structural motifs that characterize each extremal “type” of association. Though the graphical imagery primarily allows for the visualization of the angular features of each complex, the progressive shortening of phosphorus–halogen bond length RP···F or RP···I (as well as concomitant elongation of dihalogen bond length RIF) with increasing binding energy ΔEbind is also discernable in the details of the figure panels. From both Table 1 and Figure 1, one can recognize that the structural effects of different R-substituents are modest compared to that of reversing dihalogen polarity.

Figure 1

Optimized structures of RPH2···IF and RPH2···FI complexes (cf. Table 1).

Figure 2 presents a graphical display of the trends in tuning binding energy ΔEbind with various R-substituents for the strongly bound complexes of the RPH2···IF type. The graph includes comparison values for the alternative theory levels given in Tables S1–S3 of SI. From Figure 2, one can observe that significant quantitative differences in binding energies are calculated at the various theory levels, but the general trend with respect to the R-substituent variation is reasonably consistent in all cases. A slight exception is seen in the reduced cyano substituent effect in the “unmixed” B3LYP/LANL2DZ, MP2/LANLDZ theory levels, which may reflect the need for diffuse functions (present in the mixed-PP model but absent in uniform LANL2DZ) to better describe the strong electron-withdrawing effects of this substituent.

Figure 2

Calculated trend of ΔEbind with substituent R in the RPH2···IF complexes for four theoretical levels. Note the reversed ordering of the cyano substituent effect on binding energy when diffuse functions are absent (“pure” LANL2DZ; see text).

Figure 3 shows a corresponding plot of dihalogen frequency shift ΔνIF with respect to the R-substituent (including comparison results from alternative methods). As shown in the figure, ΔνIF is strongly red-shifted for the most strongly bound complexes, similar to the well-known spectroscopic signature of H-bonding. Figure 4 exhibits the correlation plot of ΔνIF (Figure 3) versus ΔEbind (Figure 2) in the RPH2···IF species, showing the high correlation (Pearson |χ|2 = 0.99) between these experimentally observable quantities that confirms the similarity to H-bonding.

Figure 3

The trend of frequency shift, ΔνIF, with the R-substituent in the RPH2···IF complexes for four theoretical levels (cf. Figure 2).

Figure 4

Correlation plot of frequency shift ΔνIF versus binding energy ΔEbind for the RPH2···IF species (B3LYP/mixed-PP level), showing the excellent least-squares regression fit (dashed line) and Pearson correlation coefficient (|χ|2 = 0.99) for these experimentally measurable properties.

Returning to the structural data of Table 1 and Figure 1 for the strongly bound RPH2···IF complexes, one may note that the intermolecular P···I distances vary in the range 2.86–3.03Å. Table 2 displays empirical van der Waals (vdW) radii (Pauling values [59]) for the P, I, and F atoms of present interest, showing that the optimized RP···I distances in these complexes are all more than 1Å inside nominal van der Waals contact. This deep penetration into the exponentially repulsive sub-vdW region implies the presence of powerful intermolecular attractive forces (also of exponentially increasing strength) to achieve the net binding energies displayed in Table 1 and Figure 2. In addition, the circa 0.2Å variation in penetration distance with the R-substituent testifies to the impressive strength of substituent effects that further modulate ΔEbind binding affinity and RP···I separation in the deep sub-vdW region.

Table 2

Van der Waals radii of P, I and F atoms (Pauling, Ref. [59]).

Atom	P	I	F	P + I	P + F	I + F
VDW radii/Å	1.9	2.15	1.35	4.0	3.2	3.5

The computed structural and binding characteristics noted above identify the propensities and patterns to be explained and provide helpful clues to the nature of the underlying attractive forces. Even the more weakly bound complexes of the P···FI type exhibit optimized P···F distances (2.1–2.6Å) that lie significantly (>0.5Å) inside vdW contact (circa 3.2Å, Table 2). Indeed, in both strong RPH2···IF and weak RPH2···FI complexes, the structural features hint at striking resemblances to familiar H-bonding phenomena, including deep sub-vdW penetration (with concomitant covalent bond RIF elongation) and pronounced “linearization” tendency (albeit along different 3-center axes: P···IF in the RPH2···IF case versus RP···F in the RPH2···FI case). We take up NBO/NRT investigation of these apparent similarities and bonding relationships in the following sections.

3. NBO/NRT Descriptors

As suggested in (1a,b) above, an NBO/NRT analysis of B:···HA hydrogen-bonding phenomena reveals a host of characteristic correlations [22,23,60] that implicate resonance-type nB→σ*AH “charge transfer” interactions as the fundamental electronic origin of H-bonding, consistent with the earliest results of NBO analysis [61,62,63]. This conclusion also supports a considerable variety of other experimental and theoretical evidence [64,65,66,67,68] for the charge–transfer nature of H-bonding. The question here is whether similar NBO/NRT diagnostics may provide analogous associations with n→σ* (lone pair to valence antibond) donor–acceptor interactions of the present species.

The analogies of the preceding section suggest appropriate NBO/NRT descriptors for the corresponding ···X— (“X-ogen”) species that can be similarly visualized in simple orbital overlap terms. For RPH2···IF complexes, the key feature of envisioned P···I—F “halogen bonding” is the donor–acceptor delocalization of electronic charge from an occupied donor NBO of P (lone pair nP) to the favorably polarized acceptor NBO (antibond σ*IF) of the P···I—F triad. For RPH2···FI complexes, where the nP-σ*IF overlap is unfavorable, the analogous feature of envisioned F···P—R “pnicogen bonding” is the donation from a fluorine lone pair (nF) to the proximal antibond (σ*PR) of the phosphine monomer. Each such n→σ* delocalization can be alternatively quantified in terms of energetic (perturbative ΔE(2)nσ* or variational deletion ΔE($DEL)nσ*), charge transfer (Qnσ*), or fractional NRT bond order (bP···I or bP···F) descriptors.

3.1. Energy Descriptors

Table 3 summarizes estimates of the ΔEnσ* stabilization energy associated with each type of delocalization in the strongly bound RPH2···IF species. The first two columns display 2nd-order perturbative ΔE(2)nσ* estimates for nP→σ*IF (“halogen bond”) and nI→σ*PR (“pnicogen bond”) contributions. The final two columns of Table 3 display corresponding ΔE($DEL)nσ* variational deletion estimates, obtained as the variational energy raising when the specific n→σ* NBO interaction is deleted from the total energy evaluation. The two estimates are seen to be in reasonable qualitative agreement for each interaction type. The ΔE(2)n(P)→σ*(IF) values range from 19 to 37 kcal/mol (in the same order as ΔEbind), whereas ΔE(2)n(I)→σ*(FI) is relatively negligible in each case, and a similar pattern is seen in the ΔE($DEL)nσ* values. We therefore expect that the stabilization energy of the nP→σ*IF interaction is the dominant attractive contribution to the structure and binding of RPH2···IF complexes, similar to the dominance of the nB→σ*AH interaction in H-bonding.

Table 3

ΔEnσ* stabilization energy estimates(kcal/mol) for halogen-type [n(P)→σ*(IF)] and leading pnicogen-type [n(I)→σ*(PR)] interactions of RPH2···IF complexes, showing perturbative ΔE(2)nσ* values in the first two columns and variational deletion ΔE($DEL)nσ* values in the final two columns (The connected R-atom of each σ*PR acceptor NBO is shown underlined in the species listing.).

Species	ΔE(2) n(P)→σ*(IF)	ΔE(2) n(I)→σ*(PR)	ΔE($DEL) n(P)→σ*(IF)	ΔE($DEL) n(I)→σ*(PR)
CH3PH2···IF	36.96	1.32	51.87	1.47
OHPH2···IF	35.99	2.12	51.88	2.60
CF3PH2···IF	22.42	1.62	28.56	1.69
CNPH2···IF	20.73	1.40	24.80	1.73
NO2PH2···IF	19.12	2.69	23.73	2.64

For the weakly bound RPH2···FI species, the situation is more complex. Due to the unfavorable polarization of σ*FI away from the nP lone pair, the single nP→σ*IF interaction that dominates bonding in RPH2···IF species becomes only one of several such competing contributions to net binding. Specifically, for the opposed F···P—R (pnicogen bonding) motif, three possible nF lone pairs (one of sigma type and two of pi type) and three σ*PY acceptors (one σ*P and two σ*PH) are within proximal interaction range as potential contributors to a complex resonance mixture. A particular interaction such as nF(σ)→σ*P therefore represents only one of the nine related n→σ* interactions that may exert leverage on angular structure and binding energy. In such a case, only a more nuanced resonance-type description (see the discussion of NRT bond orders below) can describe the complete bonding picture that provides useful correlations with experimentally measurable properties. We therefore defer further discussion of RPH2···FI complexes to a later subsection.

3.2. Graphical Orbital Overlap Imagery

A powerful aspect of NBO energetic descriptors is their intimate relationship to the intuitive visual imagery of the orbital overlap, the NBO counterpart of the well-known “Mulliken approximation” [69,70,71,72] that underlies semi-empirical MO theory. For such visualization purposes, “pre-orthogonal” NBOs (PNBOs) [49] are employed that maintain free-space symmetries (cartoon-like “textbook” shapes) of non-interacting NBOs prior to the actual distortions of the molecular environment. Table 4 displays such PNBO visualizations for the various donor–acceptor n→σ* interactions of halogen (nP→σ*IF; first column) or pnicogen (nI(π)→σ*P, nI(σ)→σ*P; final two columns) types, with corresponding numerical ΔE(2)nσ* values (Table 3) inset in the panels.

Table 4

“Pre-orthogonal” natural bond orbital (PNBO) orbital overlap diagrams (with perturbative ΔE(2)nσ* stabilization energy (kcal/mol) as inset) for RPH2···IF complexes (cf. Table 3), showing the usefulness of Mulliken-type orbital overlap visualizations in “guesstimating” donor–acceptor interaction strengths.

Species	nP→σ*IF	nI(π)→σ*PR	nI(σ)→σ*PR
CH3PH2···IF
OHPH2···IF
CF3PH2···IF
CNPH2···IF
NO2PH2···IF

For example, in the three panels for CH3PH2···IF (first row of Table 4) one can see that the “end-on” (sigma-type) overlap of σ*IF with the nitrogen lone pair nP (left panel) is much stronger than the “edge-on” (pi-type) overlap of σ*P with the in-plane iodine nI(π) (middle panel) or nI(σ) (right panel) lone pairs, in accordance with the ΔE(2)nσ* values shown in each panel. Even a casual glance at the overall patterns of the orbital overlap in Table 4 can therefore make clear why the slight variations of orbital shape with R-substitution do not materially alter the strong propensity for halogen-type (nP→σ*IF) rather than pnicogen-type (nI(π)→σ*P, nI(σ)→σ*P) bonding throughout the RPH2···IF series, as numerically displayed in Table 3.

Though the donor–acceptor ΔEnσ* stabilization energies of Table 3 display a qualitative ordering pattern similar to that of the ΔEbind binding energies of Table 1, it is evident that these estimates significantly “overshoot” the final net binding energy. This is to be expected, because each stabilizing donor–acceptor interaction is opposed by steric repulsions of the corresponding ΔEnσ donor–donor interaction. Table 5 displays orbital overlap diagrams and inset values of the pairwise estimates ΔE(pw) n-σ for repulsive “steric exchange energy” [54] between occupied NBOs of the RPH2···IF series. As shown in the panels, the two leading intermolecular repulsions refer to the nP-nI steric clash between lone pairs of proximal P and I atoms (left) and nP-σIF clash between a lone pair and bond of the halogen bond P···I—F triad. However, weaker nI-σPR, nI-σPH clashes of the alternative pnicogen bond triads (I···P—R, I···P—H) also contribute to sub-vdW steric opposition, thereby further offsetting the apparently “excessive” NBO estimate of donor–acceptor attraction needed to achieve the final equilibrium geometry and net binding energy.

Table 5

ΔE(pw)donor–donor repulsion energy estimates (kcal/mol) for leading nP-nI and nP-σIF steric clashes of RPH2···IF complexes, obtained from steric-keyword analysis.

Species	ΔE(pw)n(P)-n(I)	ΔE(pw) n(P)-σ(IF)
CH3PH2···IF
OHPH2···IF
CF3PH2···IF
CNPH2···IF
NO2PH2···IF

PNBO visualizations for occupied nP, nI, and σIF NBOs (not shown) similarly serve to rationalize the trends in repulsive donor–donor interactions (Table 5). Furthermore, it becomes evident from graphical comparisons that all such orbital visualizations for the presently considered species are highly analogous to those for the corresponding interactions of H-bonding [22], which is in accordance with the patterns and relationships inferred from numerical descriptors.

3.3. Structural, Energetic, and Spectroscopic Effects of Donor–Acceptor Deletion

Still another widget from the NBO toolbox can be used to exhibit the profound effects of nP→σ*IF delocalization on supramolecular bonding properties. Similar to the variational $DEL-deletion ΔE($DEL)n(P)→σ*(IF) evaluations of n→σ* donor–acceptor strength in the equilibrium species (Table 3), one can variationally reoptimize the structure for the associated E($DEL) n(P)→σ*(IF) potential energy surface in which n→σ* delocalization is absent, as though nature (in accordance with common textbook presentations) failed to include such interactions.

Table 6 displays results for the variational reoptimization of a simple 2D model of RPH2···IF dissociation, with the nP→σ*IF interaction deleted. Starting from the equilibrium values of the actual RPH2···IF species, the distance variables RP···I, RIF (for fixed idealized angular geometry) were optimized to find the new binding energy (ΔEbind) and shifts (ΔRP···I, ΔRIF, ΔΔEbind) that describe the hypothetical “resonance-free” RPH2···IF that lacks nP→σ*IF donor–acceptor attraction. From the resulting values, one sees that characteristic structural and energetic signatures of RPH2···IF bonding are lost when the nP→σ*IF interaction is “turned off;” the RPH2 and IF monomers retreat (by circa 1.3Å) to beyond-vdW separation, the characteristic RIF elongation (associated with vibrational red-shifting) from isolated monomer geometry essentially disappears, and the net binding energy drops precipitously (by >80%) to a remnant value that might plausibly be associated with residual “dipole–dipole” forces. Note that the monomer dipole moments or other aspects of electron density distribution are scarcely altered in the $DEL-reoptimized geometry, and the “lost” energy of nP→σ*IF stabilization is negligibly small (1.59 kcal/mol) at the reoptimized separation distance. From the results of such $DEL-deletion calculations, one recognizes nP→σ*IF stabilization to be the unique “smoking gun” that is both necessary and sufficient to bring the supramolecular complex to the actual short-range geometry and other signatures of halogen bonding, which is in direct correspondence to the analogous nB→σ*HA role in H-bonding.

Table 6

Relaxed geometrical parameters (RP···I, RIF; Å), binding energy (ΔEbind; kcal/mol), and dihalogen stretching frequency (νIF; cm−1) (with associated shifts from true equilibrium structure in parentheses) for $DEL-type variational reoptimizations, with the nP→σ*IF charge–transfer interaction deleted.

Species	RP···I (ΔRP···I)	RIF (ΔRIF)	ΔEbind (ΔΔEbind)
CH3PH2···IF	4.20 (1.34)	2.007 (−0.064)	2.32 (−12.28)
OHPH2···IF	4.29 (1.43)	2.006 (−0.060)	1.31 (−11.05)
CF3PH2···IF	4.34 (1.35)	2.003 (−0.035)	0.96 (−6.22)
CNPH2···IF	4.37 (1.34)	2.003 (−0.031)	0.63 (−5.30)
NO2PH2···IF	4.40 (1.39)	2.003 (−0.028)	0.38 (−4.72)

Still another use of $DEL-deletion techniques is to prepare relaxed-scan potential curves that illustrate how intermolecular n→σ* interactions lead to net binding against the essentially repulsive potential curve for hypothetical monomers that interact without benefit of resonance-type stabilization. Figure 5 illustrates such potential curves for CH3PH2···IF dissociation, showing the potential curve for the fully interacting monomers (squares, black line, “ΔEbind”; cf. Table 1) versus that for “resonance-free” monomers that lack intermolecular n→σ* interactions (triangles, blue line, repulsive “ΔEr”). The “missing” resonance-type (purely attractive) interactions are shown as the difference curve (diamonds, red line, attractive “ΔEa”). The blue curve exhibits a feeble attractive well beyond vdW contact at RP···I ≈ 4.0Å (presumably due to “dipole–dipole” attraction) but rises steeply as a repulsive steric wall in the sub-vdW region. Already in the region of the initial sub-vdW penetration, the red curve (n→σ* resonance attraction) has achieved sufficient stabilization to oppose thermal fluctuations (kT ≈ 0.6 kcal/mol at 300K) and yield significant sub-vdW binding in the black (full potential) curve down to RP···I ≈ 2.9 equilibrium separation. Inside this equilibrium distance, the intermolecular n-σ steric repulsion overcomes n-σ* attraction to give the steeply repulsive inner wall of the full potential, resonance-shifted inside circa 2.5Å. Where such nP-σIF versus nP-σ*IF “cross-over” occurs will evidently depend on the polarization of the dihalogen bond and is therefore expected to shift to a progressively larger RP···I separation (and weakened ΔEbind) for other members of the dihalogen series, thereby appearing as a secondary “electrostatic” modulating effect on the overall halogen bonding phenomenon.

Figure 5

Potential curves for CH3PH2···IF dissociation (B3LYP/mixed-PP level), showing (i) repulsive “resonance-free” potential ΔEr (blue, triangles) with all resonance-type intermolecular delocalizations deleted, (ii) the stabilizing potential for attractive resonance-type interactions ΔEa (red, diamonds) that were deleted in the resonance-free model, and (iii) the combined “full” potential ΔEbind (sum of (i), (ii)) (black, squares).

3.4. Charge and Polarity Descriptors

Table 7 displays various NBO charge/polarization descriptors for the RPH2···IF species (net intermolecular charge transfer QCT, atomic charges qP, qI, and qF of the P···I—F triad, IF bond ionicity iIF [73]), compared with those for isolated IF monomer. As shown for the dihalogen monomer in the final column of the table, I–F bond ionicity (−1 ≤ iIF ≤ +1, related to dipole moment direction from I to F) tends to increase steadily with strength of RPH2···IF binding or net QCT. This trend reflects the expected dominance of the nP→σ*IF interaction, which is evidently enhanced (cf. overlap diagrams of Table 5) if σIF bond ionicity increases toward F, thereby polarizing the σ*IF antibond toward I (and adjacent nP). The combined effect of repolarization (iIF) and net charge transfer (QCT) is to give a somewhat irregular pattern to individual atomic charges qI, qF, but all such IF charge descriptors are seen to properly “add up,” as they must according to the strict logic of a natural population analysis.

Table 7

Total RPH2 → IF charge transfer (QCT; e), atomic charges (qP, qI, qF), and σIF bond ionicity (iIF) in RPH2···IF complexes, compared with isolated IF monomer.

Species	Q CT	q P	q I	q F	i IF
CH3PH2···IF	0.2405	0.4084	0.3227	−0.5632	0.496
OHPH2···IF	0.2432	0.7345	0.3144	−0.5576	0.486
CF3PH2···IF	0.1413	0.2974	0.3857	−0.5270	0.485
CNPH2···IF	0.1255	0.3493	0.3924	−0.5179	0.480
NO2PH2···IF	0.1137	0.4496	0.4013	−0.5150	0.478
IF	-	-	0.4691	−0.4691	0.469

In the RPH2 monomer, the phosphorus charge qP reflects the still more complex effects of sigma-type induction versus intra- and intermolecular lone pair delocalizations. One can see in the table the evident effect of bonding the phosphorus to more electronegative oxygen (qP ≈ 0.73) or nitrogen (qP ≈ 0.45) rather than carbon (qP ≈ 0.3–0.4). However, also important in the present context is the intra- versus intermolecular competition for charge donation from the “busy” phosphorus lone pair nP to various pi-type acceptor orbitals of the R-substituent (e.g., πCN, πNO) versus the intermolecular nP→σ*IF interaction of principal interest.

3.5. NRT Bond Order Descriptors

Perhaps the most useful and general descriptors of supramolecular bonding are NRT bond orders {bAB}, which balance the nuances of multiple resonance structure contributions to give a single composite measure of A···B “connectivity.” As described elsewhere [40,41,42,43], NRT bond orders differ from those of the Pauling type, particularly in their close association with the first-order reduced density matrix [74] (rather than assumed valence bond (VB) wavefunction) of the system. A vastly improved convex solver [75] for the NRT variational objective function is a leading feature of the NBO 7.0 program [57].

The first three columns of Table 8 show calculated NRT bond orders for the key P/I/F atoms in both the RPH2···IF and RPH2···FI series. The values vary rather continuously over the range 0 ≤ bAB ≤ 1, with no apparent “dividing line” between what should be described as the “covalent bond” (stick) versus the “noncovalent interaction” (dots). Most surprising in these species are the non-vanishing bond orders for the “long-bond” [76,77,78] interactions of P^F (in P···I—F triads) or P^I type (in P···F—I triads). Though negligibly small in the most strongly bound RPH2···IF species, certain P^I long bonds achieve bond orders that rival or exceed those of familiar bP···F or bIF linkages (e.g., in OHPH2···FI).

Table 8

Natural resonance theory (NRT) bond orders and bond order sums for the bP···I halogen bond, bIF covalent bond, and bP^F long bond of RPH2···IF complexes (P···IF triad; upper rows) or the bP···F pnicogen bond, bIF covalent bond, and bP^I long bond of RPH2···FI complexes (F···PR triad; lower rows). Note the near-constant “bond conservation” of the summed bond orders within each 3c/4e triad (final column).

Species	b P···I	b IF	b P^F	bP···I + bIF	bP···I + bIF + bP^F
CH3PH2···IF	0.47	0.60	0.00	1.07	1.07
OHPH2···IF	0.45	0.59	0.00	1.04	1.04
CF3PH2···IF	0.31	0.64	0.11	0.95	1.05
CNPH2···IF	0.28	0.67	0.10	0.94	1.04
NO2PH2···IF	0.32	0.67	0.10	0.99	1.09
Species	bP···F	bIF	bP^I	bP···F + bIF	bP···F + bIF + bP^I
CH3PH2···FI	0.20	0.59	0.27	0.80	1.06
OHPH2···FI	0.34	0.45	0.36	0.79	1.14
CF3PH2···FI	0.09	0.89	0.07	0.98	1.05
CNPH2···FI	0.08	0.92	0.05	1.00	1.05
NO2PH2···FI	0.16	0.80	0.12	0.96	1.07

The final two columns of Table 8 show two-term or three-term bond order sums that test the validity of resonance-type “bond order conservation” [48] rules (cf. (1a,b) for H-bonding). The two-term sum (column 4) exhibits significant deviations from the expected unit value in a simple two-resonance model for a mono-valent central atom of the 3c/4e triad. These deviations indicate the need for a three-term sum (three-resonance model) that includes long-bond contributions (column 5). Remaining deviations from unity in the three-term sums of column 5 indicate additional conjugative or hyperconjugative interactions with the R-substituent that require still higher-order resonance couplings. All these results suggest the close relationship of NRT bond orders to the familiar resonance-theoretic concepts of the organic chemist, as well as their usefulness in expected empirical correlations with experimentally measurable structural, energetic, and spectroscopic properties, as discussed in the following section.

4. Mutual Correlations of NBO/NRT Descriptors with Measurable Properties

The principal NBO descriptors of donor–acceptor interactions (such as ΔE(2)nσ*, ΔE($DEL)nσ*, Q nσ*,...) intrinsically focus on orbitals of a specific n→σ* “delocalization” from a specific Lewis structural bonding pattern. Such orbital-specific descriptors are valuable when a single “parent” Lewis structure and “child” n→σ* delocalization clearly dominate quantum mechanical descriptions of the structural, energetic, and spectroscopic properties of the chosen system. However, for many chemical systems of interest (including the RPH2···FI species considered here), multiple n→σ* interactions come into play, corresponding to contributions of alternative resonance–structural bonding patterns that may no longer be “child-like” compared to a reference parent pattern. Such cases demand the more nuanced NRT descriptors such as interatomic bond orders {bAB} that balance the many possible orbital interactions contributing to the overall resonance hybrid.

As could be anticipated from the empirical origins of bond order descriptors in the pre-quantum mechanical conceptions of “mesomerism” theory [43], the “averaged” (fractional) character of NRT bond orders {bAB} and their association with atoms (rather than orbitals) makes them better adapted to describe empirical correlations with measurable properties. Well-known examples include bond order–bond length (bAB-RAB) [79], bond order–bond energy [80], and bond order–bond frequency (Badger’s rule) [81] relationships. Here we wish to briefly examine the mutual correlations of NBO/NRT descriptors in the single strong n→σ* limit (where orbital-type descriptors may suffice) as well as the more general correlations of NRT bond orders with experimentally measurable properties, spanning the full set of RPH2···IF/FI species considered in this work.

Starting with the RPH2···IF series, in which the nP→σ*IF interaction plays a clearly dominant role, we first consider the mutual correlations of ΔE(2)nσ*, ΔE($DEL)nσ*, Qnσ*, bP···I descriptors for this interaction. Taking perturbative ΔE(2)nσ* as the base descriptor for the group, we display (a) ΔE(2)nσ*-ΔE($DEL)nσ*, (b) ΔE(2)nσ*-Qnσ*, and (c) ΔE(2)nσ*-bP···I correlation diagrams in successive panels (a)–(c) of Figure 6. The displayed least-squares regression lines and Pearson |χ|2 correlation coefficients demonstrate the high quality of these mutual correlations: In the range 0.94–0.99 for NBO-specific descriptors, but slightly lower (0.94) for correlation with the NRT bond order bP···I (presumably due to the proper inclusion of secondary donor–acceptor interactions with R-substituents that are present only in the latter). The results show that descriptors of this group could be chosen rather interchangeably for correlations with measurable experimental properties when a single donor–acceptor interaction such as nP→σ*IF is clearly dominant.

Figure 6

Correlation diagrams for (a) ΔE(2)nσ*-ΔE($DEL)nσ*, (b) ΔE(2)nσ*-Qnσ*, and (c) ΔE(2)nσ*-bP···I descriptors of RPH2···IF species, showing least-squares regression line and Pearson correlation coefficient |χ|2 for each pairing.

For more general correlations with measurable experimental properties such as (a) intermolecular RP···I or RP···F, (b) intramolecular RIF, (c) binding energy ΔEbind, or (d) infrared stretching frequency νIF, we employed the NRT bond order bAB (of specified atoms A, B) as base descriptor in the correlation diagrams of Figure 7a–d, extending the correlations to both RPH2···IF (black; approximate two-resonance) and RPH2···FI (red: Multi-resonance) series. As examples of bond order­bond length correlations, Figure 7a displays the intermolecular bP···X-RP···X correlations (X = I, F) for each series, with reasonably high correlations (|χ|2 = 0.92) in each case. Figure 7b similarly displays the intramolecular bIF-RIF correlation for the dihalogen monomer, with still higher correlation coefficients (|χ|2 = 0.94 for RPH2···IF, |χ|2 = 0.98 for RPH2···FI species). As examples of bond order–bond energy correlation, Figure 7c displays the bP···X-ΔEbind correlations (X = I, F), with similar correlation coefficients (|χ|2 = 0.94) in each series. Finally, as bond order–bond frequency examples, Figure 7d displays the bIF-ΔνIF correlations, again with high correlation coefficients (|χ|2 = 0.95 for RPH2···IF, |χ|2 = 0.98 for RPH2···FI species). All these correlations suggest the high predictive utility of NRT bond orders for both the tuning effects of R-substituents and the polarity variations of the various dihalogen monomers that govern the broad range of binding energies in these supramolecular species.

Figure 7

Bond order–bond property correlation diagrams for (a) bP···X-RP···X, (b) bIF-RIF, (c) bP···X-ΔEbind, and (d) bIF-ΔνIF descriptors of RPH2···IF (black) and RPH2···FI (red) species, showing least-squares regression line and Pearson correlation coefficient |χ|2 for each pairing.

5. Summary and Conclusions

We have computationally investigated the nature of supramolecular bonding in a series of R-substituted phosphine···dihalogen complexes (RPH2···IF/FI, R = CH3, OH, CF3, CN, NO2), focusing on the electronic origin of what might be identified as a “halogen bond,” “pnicogen bond,” or some combination of both. For this purpose, we employed a broad variety of natural bond orbital (NBO) and natural resonance theory (NRT) descriptors, searching for relationships to hydrogen bonding or other known forms of intermolecular attraction (often loosely termed “noncovalent”).

Our results show that resonance-type “n→σ*” delocalization (leading to fractional intermolecular bond orders) is the dominant feature of halogen-type (e.g., nP→σ*IF) or pnicogen-type (e.g., nF→σ*PR) bonding motifs, just as for the dominant nB→σ*HA interaction of B:···H—A hydrogen bonding. We obtained evidence both for resonance-type “bond conservation” rules and the complete set of mutually consistent correlations of NBO/NRT descriptors with experimentally measurable properties of RPH2···IF/FI complexes that closely parallel those previously demonstrated for H-bonded complexes. We also demonstrate that the removal of n→σ* interactions obliterates the signature features of halogen/pnicogen bonding, whereas such features persist (albeit in muted form) even if envisioned “dipole–dipole” contributions are reversed (by reversing dihalogen monomer direction). Thus, our results establish that resonance-type n→σ* stabilization is both necessary and sufficient for the characteristic structural, energetic, and spectroscopic features of halogen or pnicogen bonding, as was previously demonstrated for H-bonding [22,23]. The present results are fully consistent with previous studies of orbital-interaction effects in halogen or pnicogen bonding [82,83]. Connections can also be seen to more nuanced interpretation of the “σ-hole” in terms of chemical “lone pair” interactions with the “extension of one of the covalent bonds,” “dative sharing (coordinate covalence),” and other quantal phenomena [84].

Based on these results, we anticipate that similar unique associations with resonance-type n→σ* “fractional bonding” will be found for chalcogen bonds, tetrel bonds, and all other such supramolecular bonding phenomena. Our conclusions thereby extend “resonance covalency” concepts to the entire supramolecular domain of sub-integer bond orders, challenging the common “electrostatic” assumptions that underlie current empirical force field modeling and textbook expositions of supramolecular chemistry.