Weak Interactions between Trivalent Pnictogen Centers: Computational Analysis of Bonding in Dimers X 3 E∙∙∙EX 3 (E = pnictogen, X = halogen)

The nature of weak interactions in dimers X 3 E∙∙∙EX 3 (E = N  Bi, X = F  I) was investigated by wave function and density functional-based methods. Out of the twenty systems studied, ten are found to be bound at the CP-MP2 and LMP2 levels of theory. Detailed partition of the interaction energy into different components revealed that dispersion is the primary force holding the dimers together but there also exist an important ionic component whose contribution increases with increasing halogen size. As expected, standard density functionals fail to describe bonding in the studied systems. However, the performance of DFT methods can be easily improved via empirical dispersion correction though full agreement with high level ab initio results was not obtained. Total binding energies calculated at the SCS-MP2 and LCCSD(T) levels of theory yield an energy scale of 10  15 kJ mol  1 which compares to a weak hydrogen bond and demonstrates that E∙∙∙E interactions, and P∙∙∙P interactions in particular, can be considered relevant for determining molecular structure in the solid state. In addition to high-level energy estimates, results from detailed bonding analysis showed that group 13 dimetallenes are structural analogues of the studied dimers, and as such contain a slipped  -interaction which is anti-bonding in nature.


Introduction
Intermolecular forces, van der Waals forces in particular, play an important role in determining the physical properties of matter. 1 For example, the melting point, boiling point, vapor pressure, viscosity, surface tension and solubility of a substance are all related to the strength of attractive forces between molecules. In addition, van der Waals forces impart chemical properties by influencing the reactivity and the chemical stability of both molecules and their complexes. Furthermore, intermolecular forces are largely responsible for crystallographic packing in solids, and they often fix the structures of molecules into specific conformations also in other phases. In case of macromolecules like enzymes, proteins and DNA, this is crucial for their biological activity. Thus, van der Waals forces have an unparalleled influence on life as we know it.
The standard definition of van der Waals forces includes the dipole-dipole, dipole-induced dipole and dispersion (London) forces though the term is occasionally used synonymously with dispersion. 2 The dispersion force arises from the transient multipoles that all molecules possess: as a result of the nonstatic nature of any electronic distribution, there is a high chance that electron density will be unevenly distributed in space which further gives rise to attractive forces between molecules even if they are nonpolar. 1 Though its strength increases with increasing molecular size, the dispersion interaction is generally the weakest of all intermolecular forces and it is effective only over short distances as the energy term follows r 6 dependence. Illustrative examples of non-covalent interactions with a dominant dispersion contribution are   and CH 4 interactions. 4 Theoretical calculations have shown that dispersion also gives an important contribution to the energy of a hydrogen bond between water molecules. 5 A particularly strong manifestation of the dispersion force is the aurophilic interaction: the propensity of Au(I) centers to form dimers, oligomers, chains and two dimensional sheets via gold-gold contacts. 6 These interactions are strengthened by relativistic effects and an energy range of 3050 kJ mol 1 has been established for a variety of Au(I)•••Au(I) bonded systems by temperature-dependent NMR. Other pairs of metal ions having a closed d 10 (or d 10 s 2 ) configuration show similar, albeit weakened bonding. 7 As a whole, the metallophilic interactions constitute a relatively new design element in supramolecular chemistry that has received intense interest over the last decade, resulting in its utilization in the synthesis of functional multidimensional systems. 8 Lately, we reported the crystal structure of the aminotetra(phosphine) ligand p-C6H4[N(PCl2)2]2 (1) in which short intermolecular contacts exist only between phosphorus atoms at the terminal NPCl2 moieties (P•••P distance 3.54 Å). 9 The absence of other intermolecular contacts and the slightly undulating two-dimensional sheet structure observed for 1 (see Figure 1) indicate that the weak dispersion interaction has a strong influence in determining the packing of molecules in the solid state. Interestingly, the fluorine analogue of 1 displays a totally different packing wherein the P•••P interactions are replaced by intermolecular F•••H, F•••P and F•••C contacts which are clearly of electrostatic origin. To investigate this phenomenon deeper, preliminary theoretical calculations were performed for simplified model dimers X3P•••PX3 (X = F, Cl). In agreement with experimental observations, a binding energy of approximately 5 510 kJ mol 1 was calculated for the chlorine derivative, whereas no P•••P interaction was found in case of its fluorine analogue.

Figure 1. Intermolecular P•••P interactions in 1.
There exists one more crystal structure in the Cambridge Crystallographic Database in which a halogenated tri-coordinate phosphorous moiety (PX2) forms an intermolecular P•••P interaction shorter than the sum of van der Waals radii. 10 Analogous close contacts have also been observed for other group 15 atoms: NF2 (3.07 Å), 11 AsI2 (3.37 Å), 12 SbF3 (3.92 Å), 13 SbCl3 (3.78 Å), 14 SbCl2 (3.673.83 Å) 15 and BiCl2 (3.893.98 Å). 15a Although the presence of abnormally short pnictogen•••pnictogen distances in these structures has been noted, neither their origin nor possible strength has been speculated in any of the published reports. In addition to the available X-ray data, there exist spectroscopic evidence that, in the liquid phase, the pnictogen trihalides EX3 (E = P, As, Bi; X = Cl, Br) form ethane like dimers at certain temperature ranges. 16 These dimers are structurally equivalent to the local C2h symmetric bonding arrangement that 6 the various EX2 moieties adopt in the solid state. 17 Short E•••E interactions have also been observed between nonhalogenated pnictogen centers. For example, a diphenyldiphosphino methane analogue with pentafluorophenyl substituents on phosphorus yields a dimeric structure with a P•••P distance of 3.32 Å. 18 In addition, distibenes and dibismuthines are known to exist as associated dimers in the solid state. 19 To our knowledge, there are no systematic studies on the nature of van der Waals interactions in X3E•••EX3 dimers (E = pnictogen, X = halogen) though E•••E attraction in related dimers (H2E-EH2)2 has been analyzed in detail. 20 In contrast, the ammonia dimer has received considerable interest. 21 In spite of the somewhat contradicting results concerning its global minimum, it has been unequivocally confirmed -both by theory and experiment-that the potential energy surface is extremely flat and that molecular interactions are dominated by hydrogen bonding. A few theoretical studies discussing the phosphine dimer have also been reported. 21a,22 Its minimum energy structure at the MP2 level of theory is of bifurcated type (C2h) with P•••P distance just little shy of the sum of van der Waals radii and binding energy of only 2 kJ mol 1 . Hence, unlike in ammonia, the interactions in phosphine are solely due to dispersion as expected taking into account the small electronegativity difference between phosphorus and hydrogen.
Taken as a whole, the available experimental and theoretical evidence indicates that the P•••P dispersion interaction between two -PCl2 moieties selectively determines the packing of molecules in the crystal structure of 1. 9 This raises an important question whether the observed P•••P interactions -or E•••E interactions in general-are strong enough to be more widely utilizable in crystal engineering or in building supramolecular 7 assemblies. In order to obtain a conclusive answer, the first theoretical analysis of bonding in doubly bifurcated C2h symmetric X3E•••EX3 dimers (26) is presented herein.
Our main objective is to find out which of the studied systems are bound by dispersion and how strong the interaction is in each particular case. Next, the dependence of intermolecular interaction energy on the E•••E distance is mapped using high-level calculations and the results obtained for the model systems are extrapolated to explain the structural features observed experimentally. As a conclusion, a link from X3E•••EX3 dimers to the electronic structures of group 13 and 14 dimetallenes will be made.

Computational Details
The geometries of X3E•••EX3 dimers were fully optimized in C2h symmetry using RHF, MP2, counterpoise corrected 23 CP-MP2 and local correlation LMP2. 24 For comparison, optimizations were also performed at the DFT level. The standard PBE1PBE 25 and TPSSTPSS 26 exchange correlation functionals as well as the PBEPBE 25a-c functional 8 augmented with empirical dispersion correction, 27 PBEPBE-D, were used. Counterpoise correction was not applied in DFT calculations due to the fact that as long as basis sets of at least polarized triple- quality are used, basis set superposition effects remain negligible within the DFT formalism of electronic structure theory.
Frequency analyses were performed for all stationary points found expect for those calculated at the CP-MP2 level due to the high computational cost of the method (numerical evaluation of energy derivatives). Single point energy calculations employing the optimized structures derived from the LMP2 calculations were done at LCCSD(T) and SCS-LMP2 28 levels of theory to obtain accurate interaction energies. For selected dimers, a relaxed potential energy scan in which the monomer separation was systematically varied over a range of values was performed employing the SCS-LMP2 Hamiltonian.
The cc-pVTZ and aug-cc-pVTZ basis sets were used for nuclei up to row 5; 29 small core ECP basis set of similar valence quality, namely (aug-)cc-pVTZ-PP, were used for the heavy nuclei Sb, I and Bi. 30 Density fitting approximation with auxiliary basis sets of triple- quality was used to speed up all local correlation calculations. 31 The local correlation methods employed the Pipek-Mezey localization scheme 32

Results and discussion
To obtain a proper description of the dispersion interaction, it is important to use a method capable of treating dynamic electron correlation. In most cases, MP2 yields a semi-quantitative picture at best, whereas quantitative accuracy necessitates the use of either coupled cluster or configuration interaction based methods. 37 In addition, the employed basis set must be at least valence triple- in quality and augmented with both polarization and diffuse functions. The latter set is especially important in this regard as addition of functions with small exponents allows a better description of correlation effects in outer regions i.e. parts of the electron density that are distant from the nuclei. In addition, significant error can also be introduced to the results if the basis set superposition error (BSSE) is not treated appropriately. 38 There are many more primitives available to describe an interacting molecular assembly (dimer) than its individual components (monomers) which leads to overestimation of binding energies and interaction distances that are too short. In case of weakly bound systems, the most common approach to treat BSSE a posteriori is the counterpoise correction (CP) method which requires additional energy calculations for the individual components in the orbital basis of the entire molecular assembly. 23 Another approach for obtaining structures and energies that are virtually BSSE free is to use local correlation methods that avoid the pitfall of their canonical equivalents simply by construction. 24 That is, in local correlation methods the orbitals are localized and excitations are restricted only to certain domains which are subspaces of the projected atomic orbitals, spatially close to the orbitals from which the electrons are being excited.
The ability of DFT to describe different van der Waals interactions has been a topic of significant amount of discussion and therefore deserves a comment as well. The general conclusion is that the various standard density functionals perform the better the stronger the interaction. 39 For example, the most popular gradient corrected and hybrid functionals are able to describe hydrogen bonding with accuracy comparable to the MPX (X = 24) methods. On the other hand, since the current exchange correlation functionals are local, they are by nature incapable of modeling dispersion which is a purely nonlocal phenomenon. 1 Both rigorous theoretical solutions as well as ad hoc corrections to resolve the "DFT dispersion problem" have been suggested. Two of the more recent methods showing satisfactory and uniform performance are the DFT-D 27 ansatz involving empirical energy correction terms and the range-separated hybrid 40 approach in which the correlation is divided into long and short range parts which are treated by MP2 and local/semi-local density functionals, respectively. In the current contribution, the former of the two aforementioned methods was chosen.

Optimized geometries
Since Hartree-Fock is unable to describe dynamic electron correlation, it comes as no surprise that it predicts all doubly bifurcated C2h dimers to be unbound with the exception of 2a. For this system, the (unphysical) stabilization via basis set superposition error suffices to balance the repulsive N•••N interactions, leading to optimized distance of 3.66 Å and interaction energy of 0.5 kJ mol 1 . When using counterpoise corrected optimization procedure, 2a becomes unbound at the RHF level.
The inclusion of electron correlation changes the picture dramatically as MP2 predicts 11 of the studied dimers to be minima on the potential energy hypersurface.
However, three of these structures are unbound when BSSE is corrected in geometry optimizations either via the counterpoise procedure or using the LMP2 method (see Table   1). In addition, even though 2c and 2d are transition states with respect to NX3 rotation at the MP2 level of theory, they appear as local minima when the LMP2 method is used. This is due to the fact that the pure MP2 method overbinds the dimers, bringing the halogen centers to close proximity and increasing the repulsive X•••X interactions. As expected, the MP2 method predicts abnormally short E•••E distances also for all other systems, the deviations from CP-MP2 and LMP2 results increasing together with the basis set size and interaction strength. Conversely, the CP-MP2 and LMP2 results are in reasonable agreement with each other, the latter values generally trailing the former from above by around 0.05 Å. Since the basis sets and convergence criteria used in CP-MP2 and LMP2 calculations were of similar quality, the observed discrepancy is attributed to the fact that the ionic component of the correlation energy is underestimated in the local picture due to restrictions in the allowed excitations .24 In light of the MP2 results discussed above, the performance of the two standard density functionals seems extremely poor. The hybrid PBE1PBE method predicts only three of the studied dimers to be minima on the potential energy hypersurface (3c, 3d and 5d), whereas only two minima are found using the TPSSTPSS GGA functional (2c and 3c). In addition, the DFT optimized structures have significantly elongated E•••E distances which in the case of TPSSTPSS functional even exceed the sum of van der Waals radii for the corresponding atoms. These findings reflect the well-known fact that standard density functional methods show varying performance for systems dominated by the dispersion interaction. 27,39,40 While reasonable results can in some cases be obtained, they do not result from proper description of the underlying physics.

13
In contrast to the poor performance of the standard DFT approaches, the PBEPBE functional augmented with empirical dispersion correction factors (PBEPBE-D) performs much better (see Table 1). For dimers containing nitrogen and phosphorus, the DFT-D approach gives optimized geometrical parameters similar to those at the LMP2 level of theory. Interestingly, DFT-D predicts 2d to be a transition state with respect to formation of cyclic, dipole bound, conformation. In addition, the results for both arsenic and antimony series deviate significantly from the MP2 values: calculated E•••E distance for 4c is overestimated by 0.30 Å and not any of the structures represents a true minimum on the potential energy surface. 41 In total, only six local minima were found at the PBEPBE-D level of theory cf. 10 for LMP2. The differences between PBEPBE-D and MP2 for systems containing heavier nuclei can be, at least in part, attributed to inefficiencies in the atomic parameterization employed. To the best of our knowledge, there are no published reports detailing the use of DFT-D for systems with main group elements from rows four and five. Another possibility is that LMP2 leads to false minima in case of dimers 4 and 5 due to overestimation of the binding energy (dispersion). However, we note that single point calculations at SCS-MP2 and LCCSD(T) levels of theory show that the error in LMP2 results is rather uniform in all systems (see below). Geometry optimizations and subsequent frequency calculations at couple cluster level would be needed to fully resolve this issue.
The LMP2 and CP-MP2 results in Table 1  In order to unambiguously distinguish a local minimum from the global, it would be necessary to perform a full conformational search. We note here that, though not explicitly the topic of the current study, we tried to locate other minima for systems 3bd using both LMP2 and PBEPBE-D methods in combination with aug-cc-pVTZ-(PP) basis sets. Different conformations were used as starting points for geometry optimizations that were performed without symmetry constraints. Only the cyclic, E•••X dipole bound, conformation was found to be a true minimum on the potential energy hypersurface.
Nevertheless, the E•••E bonded doubly bifurcated geometry is lower in energy for all dimers 3bd at LMP2, PBEPBE-D as well as SCS-LMP2 levels of theory. Hence, at least for the phosphorous dimers, which are of primary interest in the current study, the doubly bifurcated geometry represents the global minimum.
Assuming that the interaction energy is directly proportional to E•••E distance, the I3P•••PI3 dimer is expected to be the strongest bound system; at 2.91 Å, the P•••P separation in 3d is as much as 0.70 Å shorter than the sum of van der Waals radii. As already discussed in the introduction, short intermolecular interactions between -EX2 moieties have been observed in the solid state for all group 15 atoms, and the published X-ray data is in fine agreement with the numbers given in Table 1 for E = N, P and As. 1012 The data in Table 1 can also be compared with those reported for dimers (H2E-EH2)2 (E = As, Sb) 20 and it is easily seen that the calculated values for X3E•••EX3 dimers fall in the same range as the previous data. However, several crystal structures displaying short E•••E contacts between -SbCl2 or -BiCl2 units have been reported 1315 despite the fact the respective model dimers are predicted to be unbound at CP-MP2 and LMP2 levels. As will be shown below, the short E•••E distances reported for these systems arise from the combined effect of molecular environment and the shape of the potential energy surface which together facilitate close packing of molecules despite the fact the interaction is repulsive in the gas phase. Molecular structures, as obtained from the X-ray data, display the combined effect of all intermolecular interactions. Hence, crystal structures do not always reflect the geometry and bonding characteristics present in the isolated system and the importance of theoretical approaches in validating the origin of short intermolecular contacts in each particular case cannot be underestimated. Table 2 lists the interaction energies of the dimers at different levels of theory using the aug-cc-pVTZ(-PP) basis sets. The CP-MP2, LMP2 and PBEPBE-D calculations utilized molecular geometries optimized with the respective method. However, interaction energies calculated with the SCS-LMP2 and LCCSD(T) methods were obtained from single point calculations employing the LMP2 minimum geometries.

Total interaction energy
The counterpoise corrected MP2 method predicts the dimers to be more strongly bound than its local variant. This was anticipated based on the trends in optimized geometrical parameters and the fact that only 90 % of the total MP2 correlation energy is within the reach of the LMP2 method (see above). Indeed the difference of CP-MP2 and LMP2 values is roughly 10 % of the counterpoise corrected MP2 energy throughout the series. Nevertheless, quantitative differences between the two sets of numbers are small, around 25 kJ mol 1 , and both methods display the same qualitative trend. That is, for all pnictogen atoms the interaction energy decreases with increasing halogen size. As expected, the most strongly bound system is 3d with 2d and 4d not too far behind. It is known that both the canonical and local variants of MP2 theory overestimate the strength of the dispersion interaction. 37 The benzene dimer is a standard textbook example where MP2 fails dramatically: irrespective of the relative orientation of the monomers, the binding energy is overestimated by more than 70% as compared to the results at coupled cluster level of theory. 34d This error can be corrected with spin component scaling (SCS) which is a simple empirical modification that scales the correlation energy contributions from parallel and antiparallel electron pairs by separate factors, correcting the overestimation at no extra computational cost. 28,42 However, the SCS formalism is not strictly ab initio. A more rigorous alternative to it is to use a coupled cluster based approach e.g. (L)CCSD(T) The performance of these two methodologies was assessed in the present study and the results are given in Table 2.
The SCS-LMP2 interaction energies are in excellent agreement with the computationally much more demanding LCCSD(T) values, which illustrates that the empirical scaling correction works perfectly also for molecules containing atoms from the lower part of the periodic table. Like in the case of DFT-D, we were unable to find references demonstrating the performance of SCS-based approaches in modeling dispersion interactions involving heavier main group elements. As the numbers in Table   2 show, the canonical and local MP2 methods predict interaction strengths which are almost twice that calculated at the SCS-LMP2 and LCCSD(T) levels. The error is naturally largest for the most strongly bound systems I3E•••EI3. Though the qualitative trend is reproduced well at all levels of theory employed, it is evident that neither of the two classical variants of MP2 is able to uncover the subtle features in the electronic structures of the systems in question.
The interaction energies calculated at the PBEPBE-D level show an interesting trend. Numeric values for the nitrogen and phosphorus dimers are in good agreement with data at SCS-LMP2 and LCCSD(T) levels, whereas results for the arsenic and antimony analogues agree better with the CP-MP2 and LMP2 methods. Particularly noteworthy is the result for 4c whose interaction energy is lower at the PBEPBE-D level despite the fact that the optimized geometry is a transition state with a significantly elongated As•••As interaction as compared to CP-MP2 and LMP2 data. Thus, the data in Table 2 clearly indicates that, akin to MP2, the empirical DFT-D approach overestimates dispersion in heavy element systems and therefore overbinds the dimers.
To ensure that the SCS-LMP2 and LCCSD(T) results are reasonably well converged with respect to basis set size, the interaction energies for the smallest dimers were recalculated using quadruple- level basis sets for all atoms. As evident from Table   2, the numeric values for the two sets of calculations differ by less than a few kJ mol 1 .
The final conclusion we can draw from Table 2 is that the binding energy for the majority of dimers 25 in a doubly bifurcated geometry is between 1015 kJ mol 1 . This is considerably more than the binding energy in a phosphine dimer (around 2 kJ mol 1 ) 21a,22 and comparable in strength to E•••E attraction in related dimers (H2E-EH2)2. 20 For comparison, the binding energy in both T-shaped (CH) and parallel displaced () benzene dimers is around 11 kJ mol 1 . 37d However, it needs to be recognized that in reallife systems E•••E interactions will be observed between two -EX2 moieties (e.g. -NPCl2 in 1). Hence, the values reported in Table 2 should be considered lower estimates of interaction energies in experimentally characterizable systems. In part, this explains why the vast majority of known molecular systems with terminal -PCl2 functionalities do not form P•••P interactions in the solid state. The primary reason is naturally the shape of the molecule which is crucial for achieving a lattice arrangement in which E•••E interactions are possible.

Energy partitioning
Local correlation methods offer the possibility to decompose the correlation contribution to the interaction energy into different classes, either by considering from which localized orbitals the excitations arise, or into which virtual orbitals these are made. 24 Following this route, it is possible to obtain the relative importance of dispersion to the total correlation energy. Another possibility to decompose interaction energy into physically meaningful classes would be to use symmetry adapted perturbation theory (SAPT) 43 and DFT-SAPT 44 in particular.
The partitioning of the correlation contribution (at CCSD/aug-cc-pVTZ(-PP) level of theory) in dimers 25 is shown in Figure 2; contributions from exchange-dispersion type excitations are not shown since they are negligible for all systems in question. As anticipated, the dispersive component constitutes the dominant contribution in all dimers.
Interestingly, its magnitude in the strongest bound systems 3d and 4d is approximately the same as observed for aurophilic interactions. For example, the calculated dispersion energies in [HAuPH3]2 and [ClAuPH3]2 are 43 kJ mol 1 and 46 kJ mol 1 , respectively. 45 For comparison, the dispersion energy in a benzene dimer varies from 7 to 11 kJ mol 1 depending on the orientation of the monomers. 46 The analysis presented in Figure 2 shows that there is also a significant ionic contribution whose relative importance increases with increasing atomic size; for 3d and 20 4d the dispersive and ionic attractions are in fact comparable in magnitude. The data can be contrasted with values reported for [HAuPH3]2 and [ClAuPH3]2 for which the ionic contributions are 37 kJ mol 1 and 28 kJ mol 1 , respectively. 45 Consequently, the total correlation contribution to the interaction energy in dimers 3d and 4d is comparable to that typically observed for aurophilic interactions, yet the total binding energy is an order of magnitude smaller for the former. The difference is readily ascribed to repulsive interactions which are vastly different in the two systems. It should also be noted here that together with Table 2, the analysis presented in Figure 2 shows that there is a strong increase of the E•••E interaction energy with the softness of the halide X. This is entirely analogous with the corresponding increase of the metallophilic attraction in perpendicular [XAuPR3]2 dimers. 47

Potential energy surface scan
As mentioned above, E•••E distances less than the sum of van der Waals radii have been observed for the doubly bifurcated As, Sb and Bi dimers in the solid state even though the corresponding model systems are found to be unbound in the gas phase. 1215 These observations indicate that the repulsive forces between monomers remain small for a reasonably wide range of E•••E distances and that they are negligible compared to other interactions present in the solid state. To gain more insight into this matter, SCS-LMP2 relaxed potential energy scans were performed for selected dimers and the calculated interaction energies were plotted as a function of the E•••E distance ( Figure 3). Figure 3, the interaction energy for 2 and 3 increases very modestly within  0.5 Å from the equilibrium distance Re. For all four dimers, the increment is less than 10 kJ mol 1 which means that the interaction energy remains negative even at R = Re-0.5 Å. If the monomers are brought further closer to each other, the interaction energy starts to increase rapidly and at R = Re-0.8 the repulsive contributions overweight the attractive in all dimers. Nevertheless, the energetically favored area is quite long for all systems in question. For dimers that are unbound in the gas phase, that is 5b and 6b, the interaction energy naturally approaches zero from above as the monomer separation is increased. However, even for these systems, the monomers can be brought into close 22 proximity without significant increase in interaction energy. This readily explains the experimentally observed distances that are slightly less than the sum of van der Waals radii. For example, the crystallographically determined Sb•••Sb distance for 5b is 3.78 Å. 14 According to Figure 3, this corresponds to interaction energy of around 3 kJ mol 1 , a value which is easily offset by e.g. crystal packing forces.

Group-13 structural analogues
The dimers investigated in the current contribution can be considered group 15 counterparts of group 13 (RE=ER) and group 14 (R2E=ER2) dimetallenes. The analogy between 26 and group 13 dimetallenes is particularly noteworthy: both set of compounds have C2h symmetric minima, weak E-E bonding interactions and monomeric structures in solution. 48 In fact, even the frontier MOs for these two types of compounds 23 are similar and show the expected and "slipped -type" orbitals i.e. two dative bonds ( Figure 4). The significance of the slipped -type orbital to bonding in RGaGaR and, in particular, [RGaGaR] 2 has been a topic of much debate. 48d,49 The orbital has a node between the two gallium centers which implies that it should be considered an antibonding combination of two lone-pair orbitals. However, when plotted with smaller contour values, the orbital appears to have some -bonding character as the lone pair on one gallium atom merges in with the contribution from the GaR bond on the other and vice versa. Thus, although digallenes RGaGaR have been described with a formal bond order of two, in reality it is closer to one. This is in good agreement with the GaGa bond length determined for ArGaGaAr (Ar = 2,6-Dipp2C6H3), 2.627 Å (cf. GaGa single bond range 2.33-2.54 Å). 48a It is evident that bonding in dimers 26 cannot be explained simply by using orbital overlap as electron correlation methods are mandatory to describe the interaction in the first place. However, if one plots the slipped -type orbital in Figure 4 with smaller contour values, it displays bonding contributions analogous to those presented for digallenes. Hence, following the reasoning presented for digallenes and digallynes, 49a,c,d the dimers 26 can be considered to formally contain a double bond! Clearly such interpretation makes no sense chemically. We therefore decided to have a closer look at bonding in digallenes to determine the true nature of the slipped -interaction.
It appears that the majority of calculations done for digallenes and digallynes have been done using either DFT of MP2. 49 However, if the slipped -type orbital displays significant bonding character, it should be present already at the HF level of theory. Our calculations show that when using the HF method, the GaGa bond length in RGaGaR is heavily dependent on the type of substituent used. The bond length in the parent digallene is 2.794 Å, whereas it is 2.918 and 3.313 Å in Me and Ph substituted species, respectively; when using 2,6-Ph2C6H3 as a substituent, the digallene is unbound at the RHF level. These results clearly demonstrate that the slipped -interaction is in general destabilizing and effectively cancels the -bond in all real-life systems with aryl substituents. Hence, were there no electron correlation, bonding in all but H and Me substituted digallenes would be negligible. However, unlike in the case of dimers 26, digallenes can be readily described with standard density functional theory, which indicates that something else than dispersion keeps the molecules from falling apart.
Preliminary calculations reveal that digallenes are bound by singlet diradical character 25 which, in wave function terminology, results in transfer of electrons from the antibonding -slipped HOMO to the -bonding LUMO. Since density functional theory describes singlet diradical character, and electron correlation in general, through exchange-correlation functional, 50 it is understandable that this aspect to bonding in digallenes has not been recognized before. 51

Conclusions
High-level quantum chemical calculations performed for model systems 26 provided insight into the nature of weak pnictogen•••pnictogen interactions observed in the solid state. The results conclusively show that the molecules adhere together primarily via dispersion but there is also a significant ionic contribution whose importance increases with increasing halogen size. The bonding interactions maximize under phosphorus and I3P•••PI3 is predicted to be the most stable dimer. Its interaction energy is approximately 15 kJ mol 1 which compares to the strength of CH and  interactions in benzene dimers. Particularly notable is the magnitude of the dispersive component in I3P•••PI3 which is on par with values typically observed in aurophilic interactions. Based on the computational as well as prior experimental evidence, we conclude that the E•••E interactions discussed herein are sufficiently strong to be considered relevant in building supramolecular assemblies and constitute an underexplored inorganic structural motif that can be used in crystal engineering. Experimental efforts to validate the presented computational results for heavier halogenides of 1 are currently ongoing.

26
Alongside with the calculated properties of the model dimers, the reported results have significant theoretical value. To the best of our knowledge, the present study is the first in which DFT-D and SCS-LMP2 are used to describe dispersion interaction involving heavy main group atoms. Hence, the current contribution is as an important benchmark. It was shown that the results given by DFT-D are in agreement with highlevel ab intio data only in case of dimers which contain the lighter pnictogens nitrogen and phosphorus. For heavier elements, DFT-D apparently overestimates the importance of dispersion, but it nevertheless predicts the doubly bifurcated dimers to be first-order transition states on the energy hypersurface. In contrast, the SCS-LMP2 method yields accurate interaction energies at low computational cost throughout the entire periodic table which implies that the empirical spin component scaling approach works perfectly even for heavy elements arsenic and antimony. As a final remark, the orbital analysis showed that digallenes are group 13 counterparts of dimers 26 and, as such, they contain a slipped -type interaction which is anti-bonding in nature. However, diradical character, not dispersion, holds group 13 dimetallanes together which underlines the fact that they are only structural analogues of 26.