Abstract

The present work starts with providing a description of the halogen bonding (XB) interaction between the halogen atom of MH3X (where M = C–Pb and X = I, At) and the N atom of HCN. This interaction leads to the formation of stable yet very weakly bound MH3X⋯NCH complexes for which the interaction energy () between MH3X and HCN is calculated using various symmetry-adapted perturbation theory (SAPT) methods combined with the def2-QZVPP basis set and midbond functions. This basis set assigns effective core potentials (ECPs) not only to the I or At atom directly participating in the XB interaction with HCN but also to the M atom when substituted with Sn or Pb. Twelve SAPT methods (or levels) are taken into consideration. According to the SAPT analysis of , the XB interaction in the complexes shows mixed electrostatic-dispersion nature. Next, the accuracy of SAPT is evaluated by comparing with CCSD(T) reference data. This comparison reveals that high-order SAPT2+ method and the much less computationally demanding SAPT(DFT) method perform very well in describing of the complexes. However, the accuracy of these methods decreases dramatically if they are combined with the so-called Hartree-Fock correction.

1. Introduction

Weak intermolecular interactions play an important role in establishing the structure and stability of a broad range of chemical systems, from simple molecular complexes to macromolecular assemblies [13]. Computational methods based on electronic structure calculations are useful in determining the magnitude of weak intermolecular interactions as well as in understanding their nature [46]. However, the ability to describe weak intermolecular interactions accurately and efficiently differs significantly between various families of computational electronic structure methods. Among these families, the symmetry-adapted perturbation theory (SAPT) [711] is capable of providing weak interaction energies accurate even up to 2–4% for benchmark systems of rare-gas dimers [9] but the computational cost of the corresponding calculations in principle scales like with the system size. This makes such calculations very time-consuming, especially if they are carried out for large chemical systems. The high computational cost of SAPT calculations is associated, to a great extent, with the inclusion of intramolecular electron correlation effects. It is particularly true for the traditional formulation of SAPT which uses many-body electron correlation wave functions for molecules forming a complex [7]. An alternative approach to handling the intramolecular electron correlation within SAPT is to utilize the density functional theory (DFT) for the molecular wave functions [1217]. The resulting SAPT(DFT) formulation scales better than the traditional formulation, and, therefore, the former can be applied to systems with hundreds of atoms [18]. In order to adjust the computational demands of SAPT calculations to the size of a given system, various methods (or levels) have been introduced for the traditional SAPT formulation [10, 19]. In general, these levels differ in their scaling behavior through the omission of some energy correction terms from the SAPT expansion of intermolecular interaction energy.

Weak intermolecular interactions sometimes occur between molecular sites, of which one contains an atom of element with high nuclear charge. Halogen bonding (XB) [20] involving iodine constitutes a good example of such interactions (historically speaking, the XB formation in H3N⋯I2 was probably the very first observation of XB interaction [21]). From the computational point of view, describing weak intermolecular interactions involving heavier elements is a real challenge for electronic structure methods for two main reasons. Firstly, the large number of electrons in atoms of heavier elements makes electronic structure calculations very costly. Secondly, the incorporation of relativistic effects in the calculations is essential to obtain accurate energies of weak interactions involving heavier elements. These difficulties can be largely overcome by the use of effective core potentials (ECPs) for atoms of heavier elements (usually with the nuclear charge ) [22]. Within the ECP approximation the explicit treatment of core electrons of heavier atoms is replaced by introducing potential operators fitted to the results of relativistic all-electron calculations for such atoms. Therefore, the ECP approximation is able to implicitly account for leading relativistic effects at low computational cost. The use of ECPs has also been implemented for SAPT calculations of interaction energies [23].

In this work a systematic examination of the accuracy of a combination of SAPT and ECP for the description of weak XB is provided relative to the level of SAPT. The accuracy of a variety of SAPT methods is established across a series of 10 heterodimers containing MH3X (where M = C–Pb and X = I, At) complexed with HCN. These complexes exemplify the X⋯N type of XB in which the halogen atom possesses a very high nuclear charge (53 and 85 for iodine and astatine, resp.). Thus, the ECP treatment of I- and At-atom core electrons is employed here for the SAPT description of the XB interaction in the MH3X⋯NCH complexes.

2. Methods

The geometries of the MH3X⋯NCH complexes have been optimized at the MP2/def2-QZVPP level of theory [2427]. These geometries are characterized by the absence of imaginary vibrational frequencies, and, therefore, they correspond to minima on the potential energy surfaces of the complexes. These geometries will also be kept fixed in all subsequent calculations of the interaction energy () between the MH3X and HCN parts of the complexes. The MP2 method has been combined with the resolution-of-the-identity approximation [28]. The def2-QZVPP basis set makes use of Stuttgart relativistic ECPs [26] to describe the core electrons of Sn, Pb, I, and At. These calculations have been carried out with the TURBOMOLE 6.6 program [29].

The values of between MH3X and HCN in the complexes have been calculated using such SAPT methods as SAPT0, SAPT2, SAPT2+, SAPT2+, SAPT2+3, and SAPT(DFT) [10, 19]. The first five methods are based on the traditional (or regular) wave function-based formulation of SAPT. They truncate the SAPT expansion of either at the second or at the third order in the perturbation operator of intermolecular interaction, and they differ in the degree of the completeness of included intramolecular electron correlation effects. Details of the truncation of SAPT expansion for individual SAPT methods are shown in the following definitions:The explanation of SAPT energy correction terms on the right sides of (1) can be found in a number of reviews covering the fundamentals of SAPT [711]. Suffice it to say here that SAPT0, SAPT2, and SAPT2+ treat the intermolecular interaction perturbation up to second order, while SAPT2+ and SAPT2+3 additionally include third-order terms in the perturbation operator of intermolecular interaction. Regarding the intramolecular electron correlation, it is completely neglected by SAPT0 and included up to second order for electrostatic, exchange, and induction contributions to in SAPT2 and additionally for the dispersion contribution by SAPT2+, SAPT2+, and SAPT2+3. The PBE0 density functional [30] with the Fermi-Amaldi asymptotic correction [31] and the Tozer-Handy splicing scheme [32] have been used in SAPT(DFT). To compute the asymptotic correction, either available experimental or calculated ionization potentials have been employed (see Supplementary Section  1 in Supplementary Material available online at https://doi.org/10.1155/2017/9031494). Both the wave function-based SAPT methods and SAPT(DFT) can be supplemented with the so-called Hartree-Fock (HF) correction ( or , depending on the SAPT method to which the correction is added). For SAPT0, SAPT2, SAPT2+, SAPT2+, and SAPT(DFT) such a correction is determined aswhereas a different form is employed for the SAPT2+3 methodThe term in (2) denotes the supermolecular Hartree-Fock interaction energy in the complex for which the SAPT interaction energy is calculated. The HF correction is a way to capture some higher-order terms not explicitly evaluated by SAPT. All SAPT results have been produced using the dimer-centered plus basis set (DC+BS) scheme, with the def2-QZVPP basis set extended with a set of midbond (mb) functions [33]. The resulting basis set will be denoted by def2-QZVPP+mb further in the text. The set of mb functions is of the 3s3p2d2f form, with the exponents 0.9, 0.3, and 0.1 for and p and 0.6 and 0.2 for and . These functions are centered midway between the X and N atoms of the complexes. All SAPT calculations have been done using the SAPT 2012.2 [34] program interfaced with DALTON 2.0 [35]. The SAPT(DFT) calculations have additionally been carried out with the SAPT(DFT) implementation available in the MOLPRO 2012.1 program [36], but these calculations do not end with successful reproduction of the results obtained with SAPT 2012.2.

Apart from the SAPT description of the XB interaction in MH3X⋯NCH, additional supermolecular estimations of are provided in this work for comparison purposes. A number of popular wave function theory methods, mainly those belonging to the family of Møller-Plesset methods [24], will be taken into account. MP2 [37], SCS-MP2 [38], SCS(MI)-MP2 [39], MP2C [40], MP2.5 [41], MP3 [37], MP4(SDQ) [37], MP4 [37], and CCSD [42] are selected. The resolution-of-the-identity approximation is used wherever possible. The aforementioned methods are combined with the def2-QZVPP+mb basis set in order to calculate the magnitude of the XB interaction in the MH3X⋯NCH complexes. The resulting values are corrected for the basis set superposition error by utilizing a counterpoise correction [43]. All the supermolecular calculations of have been performed using TURBOMOLE 6.6.

The accuracy of SAPT methods is evaluated by comparing their values calculated for MH3X⋯NCH with the corresponding reference data. The reference values of have been obtained from the CCSD(T) method [42] combined with the def2-QZVPP+mb basis set. The CCSD(T)/def2-QZVPP+mb values of are corrected for the basis set superposition error [43]. Additionally, the CCSD(T) values extrapolated to the complete basis set (CBS) limit [4447] will also serve as reference values. The CCSD(T)/CBS level of theory is commonly perceived as the “gold standard” for predicting reliable energies of weak intermolecular interactions [48]. This level of theory was previously used to provide benchmark values in a set of 40 complexes that were stabilized by a variety of weak intermolecular interactions in which halogens participated [49]. Details of the CCSD(T)/CBS calculations performed for MH3X⋯NCH are given in Supplementary Section  3.

3. Results and Discussion

3.1. Fundamental Properties of the Complexes

The geometry optimization of the complexes with the assumed X⋯N XB interaction leads to structures in which the HCN molecule is oriented precisely along the M–X bond axis (see Figure 1). The values of the distance between the X and N atoms () for the complexes in their optimized geometries are listed in Table 1. These values are rather large ( Å) but in the case of the complexes with X = I they are still shorter than the sum of the I- and N-atom van der Waals radii. This sum is equal to 3.70 Å if the van der Waals radii proposed by Alvarez [50] are taken. The fulfillment of the geometrical criterion formulated as  Å suggests the existence of XB in the MH3I⋯NCH complexes. Such a geometrical criterion cannot be applied to the complexes with X = At because the van der Waals radius of astatine is unknown [50].

The structure of the MH3X⋯NCH complexes results from the highly directional ability of the X atom, covalently bonded to the M atom, to attract the N-atom lone electron pair of HCN, acting as a nucleophilic species in this case. A widely accepted mechanism of XB formation postulates that XB originates from the interaction between two centers possessing opposite electrostatic potentials [51]. For the MH3X⋯NCH complexes the center exhibiting a positive electrostatic potential is located on the X atom in the region on the extension of the M–X bond (i.e., the -hole), whereas the negative electrostatic potential characterizes the region around the lone electron pair of the N atom. The maximal value of electrostatic potential on the molecular surface delineated by the electron density contour of 0.001 au () in the -hole region on the X atom is a convenient tool for quantifying the -hole [52]. The calculated values of are presented in Table 1. An increase in values is observed if the X atom changes from I to At, and thus it becomes heavier and more polarizable [52]. Moreover, values increase with the growing electron-attracting character of the MH3 group [53]. The electron-attracting character can be represented by the so-called group electronegativity [54] which is calculated for the MH3 groups (for details, see Supplementary Section 4). The electronegativity of the MH3 groups increases gradually while going up group 14 from M = Pb to M = C, and, therefore, the most positive values are achieved for CH3X in either of the sequences of the MH3X⋯NCH complexes with the X atom kept fixed.

The CCSD(T)/CBS values shown in Table 1 are negative for all MH3X⋯NCH complexes. This implies their energetic stabilization and this also constitutes an energetic criterion for XB occurrence in the complexes. The values of fall in a narrow range from −0.55 to −2.36 kcal/mol. The XB interaction in MH3X⋯NCH can be classified as weak in view of the fact that the strongest XB interactions reported in literature exhibit interaction energies exceeding −40 kcal/mol [55]. The CH3I molecule forms weaker XB with HCN than with NH3 [56] or CH2O [57]. The weakness of the XB interaction in the MH3X⋯NCH complexes, as indicated by their values, is connected with the large distances. If pairs of the complexes with the M atom fixed are considered, MH3At always forms stronger XB than MH3I. This is reflected in the geometrical relation for each such pair. The XB interaction becomes weaker and weaker when the M atom traverses down group 14 for two sequences of the complexes with the X atom kept fixed and varying M. The clear trend in is only roughly followed by a trend in for these sequences. These trends do not correlate with each other ideally but the monotonous weakening of the XB interaction is often accompanied by the gradual elongation of . It is noteworthy that correlates well with for all MH3X⋯NCH complexes: the greater is developed in the -hole region, the stronger the XB occurs.

In order to establish the nature of the XB interaction in MH3X⋯NCH, the energy correction terms from the SAPT(DFT) expansion of have been grouped into four principal components corresponding to electrostatics, induction, dispersion, and exchange. Details of the grouping scheme can be found in Supplementary Section5. The values of the electrostatic (), induction (), dispersion (), and exchange () components of obtained from SAPT(DFT) are listed in Table 1. The most important aspect of these results is the fact that and play the central role in the attraction between MH3X and HCN. These two components are responsible for at least 88% of the total attraction occurring in the complexes. Thus, the XB interaction in MH3X⋯NCH shows mixed electrostatic-dispersion nature. This finding is in agreement with the analogical conclusion drawn in an earlier work taking various types of XB into account [58]. It is interesting to examine trends in the attractive components for the sequences of the complexes with the X atom kept fixed. For such sequences the percentage of electrostatics decreases (and simultaneously the percentage of dispersion increases) while descending group 14 from M = C to M = Pb. This leads to the conclusion that the role of dispersion in the stabilization of the complexes is most significant for the weakest XB interaction. component does not change monotonously if M descends group 14. Moreover, the percentage share of in the total attraction between MH3X and HCN remains practically constant for all 10 complexes.

3.2. Assessment of Various SAPT Methods

In order to evaluate the performance of various SAPT methods in describing the XB interaction in the MH3X⋯NCH complexes, values predicted by these methods are compared with CCSD(T) reference results in a statistical manner. This statistical comparison is based on three metrics of errors in the SAPT values. The root mean-squared error (RMSE), the percentage relative root mean-squared error (rRMSE  (%)), and the mean signed error (MSE) are calculated using the following formulas:where is the number of the complexes studied (M = 10). RMSE is used to characterize the overall performance of each SAPT method in predicting . Because values are small themselves and the same applies to the RMSE values, then rRMSE  (%) should more conveniently and reliably illustrate the accuracy of the SAPT values relative to the corresponding CCSD(T) results. Finally, MSE provides information about systematic errors occurring in values obtained from each SAPT method. If MSE is negative (positive) for a given SAPT method, then this method tends to overestimate (underestimate) the strength of the XB interaction in the MH3X⋯NCH complexes.

The evaluation of SAPT accuracy has been performed for 12 SAPT methods. Table 2 lists these methods, together with the three mean errors determined for their values. Based on the RMSE values one can think that SAPT generally performs very well because RMSE is less than 1 kcal/mol for all methods. This is not necessarily true if rRMSE  (%) is analyzed as well. rRMSE  (%) is high for several SAPT methods, which indicates that the RMSE values are in fact quite large when compared to both the CCSD(T)/def2-QZVPP+mb and CCSD(T)/CBS values. It is justified to rank SAPT methods’ performances relative to the rRMSE  (%) values. From the magnitude of rRMSE  (%) one can conclude that the SAPT2+ method provides the best approximation of the CCSD(T) reference results. As with SAPT2+, the SAPT(DFT) method reproduces perfectly the CCSD(T)/def2-QZVPP+mb energies. A slightly larger rRMSE  (%) value is observed for SAPT2+3 but this value still does not exceed 10%. The values of MSE for the three methods indicate that SAPT2+ tends to slightly overestimate the strength of the XB interaction in MH3X⋯NCH, while SAPT(DFT) and SAPT2+3 show the opposite tendency (the mean underestimation of reaches a maximum of 0.14 kcal/mol for SAPT(DFT)). Unsurprisingly, two regular high-order SAPT methods SAPT2+ and SAPT2+3 yield with very good accuracy. It is rather puzzling, however, that SAPT2+3 performs slightly worse than SAPT2+ while the reverse ordering of these two methods could actually be expected. This means that the inclusion of all third-order interaction energy correction terms does not lead to any improvement in for the MH3X⋯NCH complexes. The SAPT2+ accuracy superior to that of SAPT2+3 is mostly due to the fact that the third-order induction and exchange-induction terms are poorly recovered by SAPT combined with ECPs [23] (see also Supplementary Section 2).

From the results in Table 2 the effect of the HF correction on can be established. This correction collects mainly higher-order induction and exchange-induction terms but it may also contain some spurious, unphysical terms. It is evident that the HF correction adversely affects the accuracies of all SAPT methods in describing the weak XB interaction in MH3X⋯NCH. Six SAPT methods utilizing this correction (the so-called hybrid approach) show many times larger RMSE values than the same methods without or . As illustrated by the negative values of MSE, the majority of SAPT methods tend to predict an overbinding between MH3X and HCN. Adding the HF correction results in further overestimation of the XB strength in the complexes. A similar overestimation of I⋯O XB strength was also detected in a previous SAPT study of the CH3I⋯OCH2 complex [57].

The deterioration of SAPT accuracy due to the incorporation of or into is likely caused by a spurious, unphysical part of the HF correction. It is known that spurious, unphysical effects contribute considerably to the HF correction if a molecular complex is bound primarily by dispersion [59]. This seems to be the case for the MH3X⋯NCH complexes for which the dispersion component of their is always significant or even dominant (the percentage of ; see Table 1).

The simplest SAPT method, that is, SAPT0, outperforms both SAPT2 and SAPT2+, although the latter methods incorporate more energy correction terms, and, therefore, they are computationally more demanding. This superior accuracy of SAPT0 seems to be due to a fairly good compensation of energy correction terms appearing in the SAPT0 expansion of . It has been demonstrated in another recent work [19] that SAPT0 yields in many cases not much worse interaction energies compared to more advanced SAPT methods which account for intramolecular electron correlation effects. The combination of SAPT0, SAPT2, or SAPT2+ with fails badly in describing the XB interaction in MH3X⋯NCH; rRMSE  (%) of the resulting hybrid levels exceeds 50% relative to the CCSD(T)/def2-QZVPP+mb reference results.

It is instructive to compare the performance of the SAPT methods with the performance of a series of popular wave function theory methods that compute within the supermolecular approach. Statistical metrics measuring the accuracy of supermolecular values obtained from 10 wave function theory methods against the CCSD(T) reference data are presented in Table 3. Of the 10 methods, SCS(MI)-MP2, MP2C, and MP2.5 perform best. They afford very small RMSE and rRMSE  (%) and these errors are even slightly smaller than those made by SAPT2+ and SAPT2+3. What is particularly important is that the SCS(MI)-MP2 calculations entail much lower computational cost (MP2 formally scales as with the system size). The good performance of SCS(MI)-MP2 has actually been expected because this method has been specially designed for weak intermolecular interactions. Similarly, the MP2C and MP2.5 methods were also developed for providing accurate interaction energies in various molecular complexes. For the MH3X⋯NCH complexes, these modifications of the original MP2 method do show a great improvement in the reproduction of over their ancestor. The performances of MP4 and CCSD are rather disappointing and these advanced methods yield less accurate values of than the regular high-order SAPT and SAPT(DFT) methods. According to the HF method, the interaction between MH3X and HCN destabilizes the MH3X⋯NCH complexes, which disqualifies this method as a computational tool for describing weak XB interactions reliably.

4. Conclusions

A description of the X⋯N XB interaction in a series of 10 MH3X⋯NCH complexes (where M = C–Pb and X = I, At) has been obtained from various SAPT methods combined with the def2-QZVPP basis set and midbond functions. This basis set assigns ECPs not only to the I or At atom directly participating in the XB interaction with HCN but also to the M atom when substituted with Sn or Pb. The accuracy of SAPT has been evaluated by comparing with CCSD(T) reference data. The main findings made in this work can be summarized as follows:(1)The CCSD(T)/CBS calculations of reveal that the complexes are stabilized by a weak X⋯N XB interaction whose values fall in the range between −0.55 and −2.36 kcal/mol. Astatine forms stronger XB interaction with HCN than iodine, by no more than 1 kcal/mol in each pair of the complexes with a given M. A relationship between and the kind of group 14 element substituting the M atom has been identified for two sequences of the complexes with the X atom kept fixed. The strength of the XB interaction decreases gradually when the M atom varies down group 14 from C to Pb. This is rationalized in terms of and MH3-group electronegativity.(2)Grouping the terms from the SAPT(DFT) expansion of into four principal components indicates the mixed electrostatic-dispersion nature of the XB interaction in MH3X⋯NCH. The XB weakening caused by the substitution of the M and X atoms clearly affects the percentage shares of and . This weakening is associated with the decrease of percentage share and the simultaneous increase of percentage share.(3)Regular high-order SAPT methods such as SAPT2+ and SAPT2+3 provide XB description very close to that obtained at the CCSD(T)/def2-QZVPP+mb and CCSD(T)/CBS levels of theory. The accuracy of obtained from the computationally favorable SAPT(DFT) method is comparable to that found for SAPT2+ interaction energies. The HF correction present in the hybrid approach leads to a significant deterioration in the accuracy of the SAPT description of the XB interaction in the MH3X⋯NCH complexes.

Competing Interests

The author declares no competing interests.

Acknowledgments

This work was partially supported by PL-Grid Infrastructure.

Supplementary Materials

The Supplementary Material attached to this article contains further details of computational methodologies used in the study (Supplemetary Sections 1-5) and additional tables with numerical data (Supplementary Section 6).

  1. Supplementary Material