Quantifying the Intrinsic Strength of C–H⋯O Intermolecular Interactions

It has been recognized that the C–H⋯O structural motif can be present in destabilizing as well as highly stabilizing intermolecular environments. Thus, it should be of interest to describe the strength of the C–H⋯O hydrogen bond for constant structural factors so that this intrinsic strength can be quantified and compared to other types of interactions. This description is provided here for C2h-symmetric dimers of acrylic acid by means of the calculations that employ the coupled-cluster theory with singles, doubles, and perturbative triples [CCSD(T)] together with an extrapolation to the complete basis set (CBS) limit. Dimers featuring the C–H⋯O and O–H⋯O hydrogens bonds are carefully investigated in a wide range of intermolecular separations by the CCSD(T)/CBS approach, and also by the symmetry-adapted perturbation theory (SAPT) method, which is based on the density-functional theory (DFT) treatment of monomers. While the nature of these two types of hydrogen bonding is very similar according to the SAPT-DFT/CBS calculations and on the basis of a comparison of the intermolecular potential curves, the intrinsic strength of the C–H⋯O interaction is found to be about a quarter of its O–H⋯O counterpart that is less than one might anticipate.


Introduction
Hydrogen bonds (HBs) [1] form a particularly important class of noncovalent interactions [2]. From among the numerous types of HBs and related hydrogen-bonding (H-bonding) scenarios [3] (see reference [4] for the most recent survey of important experimental and theoretical research into H-bonding), it is the nature of the C-H· · · O intermolecular H-bonding that is perhaps the most fervently debated (see the important study [5] and work cited therein). This type of H-bonding is known to be quite important. Recent examples of its significance include studies of crystal packing (see the related debate in reference [6]), protein folding [7], structure of certain liquids [8], zeolites [9] and interfaces [10], formation of protein-ligand complexes [11], crosslinking of biomaterials [12], and even the molecular recognition of nerve agents [13]. Hence, in addition to numerous experimental investigations, the C-H· · · O intermolecular interactions were examined by a variety of theoretical approaches. The classic work of Scheiner and his coworkers [14] on the categorization of H-bonding should be mentioned, where the respective properties of the C-H· · · O and "conventional" HBs were compared and the C-H· · · O interactions were found to be "true" HBs (see also the subsequent investigation into the cooperativity of C-H· · · O and O-H· · · O H-bonding [15]). Some of these properties, which are not detailed here, were subsequently considered by a number of computational studies, with the work of Schaefer et al. [16] and Head-Gordon et al. [17] being the most notable examples. Importantly, in a classic paper describing a combination of theory and solid-state NMR experiments [18], it was shown that a distinction should be made between C-H· · · O "contacts" (cases of a small distance between some C-H group and an oxygen atom in

Interaction Energies
First, a computational scheme needs to be established for fully reliable predictions of the intermolecular interaction energy, ∆E. The focal-point approach to the CCSD(T)/CBS ∆E estimation was most recently tested [24]. When applied to the interaction energies of systems from the S22 dataset [25], this approach was shown to provide the ∆E values that differed only negligibly from their counterparts published by Sherrill et al. as the S22B collection [26]. Consequently, the same methodology as in reference [24] is used here. Its additional checks were performed for two important systems that feature C-H· · · O intermolecular bonding, namely, benzofuran:formaldehyde [27] and formaldehyde dimers [28]. In particular, the CCSD(T)/CBS ∆E data were obtained for two geometries of the global minimum of the potential-energy surface (PES) of benzofuran:formaldehyde complex, which is a highly challenging system [26]. One geometry was taken from reference [27] (where it is denoted as "semi-experimental" in Supporting Information), while the other, which is pictured in Figure S1, was obtained at the MP2/aTZ level (the second-order Møller-Plesset theory and the standard augmented correlation-consistent polarized valence triple-ζ basis set [29,30]; see Section 4 for further details of the computational methodology adopted in this work). It should be mentioned that rotational constants {A 0 , B 0 , C 0 } were measured in reference [27]. They are {1181, 1096, 788.3} MHz when rounded to four significant digits, and an inspection of Table 1 reveals that the "semi-experimental" structure reproduced these data better than the geometry provided by the MP2/aTZ optimization. However, structural differences led to only a small difference of less than one half of kJ/mol in ∆E values predicted for the two structures (see Table 1). Data in Table 1 also show that the interaction energy estimated in reference [27] by applying the cost-effective "jun-ChS" scheme [31] agrees well with its counterpart obtained in this work. Two additional methods were employed for benchmarking purposes, because these methods were also applied to much larger clusters that feature C-H· · · O and O-H· · · O intermolecular interactions (see Section 3). One of these methods is the DFT-based B2PLYP-D3/def2-QZVPPD approach (the double-hybrid B2-PLYP functional [32,33] combined with the D3 empirical dispersion correction [34] and applied together with the QZVPPD basis set [35]), which is expected to provide highly accurate interaction energies [36]. The other method is the DLPNO-CCSD(T)/CBS technique (the domain-based local pair natural orbital approximation to the CCSD(T) method with an extrapolation to the CBS limit; see Section 4 for pertinent references and computational details). In the following, the B2PLYP-D3/def2-QZVPPD and DLPNO-CCSD(T)/CBS interaction energies are denoted as ∆E(DFT) and ∆E(DLPNO), respectively. An inspection of Table 1 reveals that the interaction energies predicted by various methods for benzofuran:formaldehyde agree with each other very well.  [27] for details. 2 Obtained in reference [31] using "jun-ChS" scheme.
The MP2/aTZ method was further applied to two dimers of formaldehyde, namely, the adducts featuring C s and C 2h symmetry. The harmonic vibrational zero-point energies, ∆(ZPE) h , were computed and added to the corresponding ∆E values to estimate the related dissociation energies, (D 0 ) h , from (D 0 ) h = ∆E + ∆(ZPE) h . The (D 0 ) h data were compared to their counterparts from reference [36] (see Table 2). Expectedly, a close agreement was found, because a very similar computational methodology was adopted in reference [37]. Moreover, the anharmonic effects on zero-point energies were approximated by the VPT2 (vibrational second order perturbation) method at the MP2/aTZ level to obtain the ∆(ZPE) a values and combine them with relevant ∆E to assess the anharmonicitycorrected dissociation energies, (D 0 ) a , as (D 0 ) a = ∆E + ∆(ZPE) a . Results are summarized in Table 2 for the sake of comparison with analogous data presented in reference [38]. Importantly, the two (D 0 ) a estimates for the C s -symmetric structure, which is the global minimum of the PES of formamide dimer [28], agree within "spectroscopical accuracy" of one kJ/mol (see Table 2). It is noted that the ∆E(DFT) and ∆E(DLPNO) values for this structure are −19.3 and −18.7 kJ, respectively, while for the C 2h -symmetric minimum they amount to −15.4 and −15.8 kJ, respectively, and thus are all highly accurate. The dimers of ε-glycine (see Figure 1) are considered in the next paragraph.    Table S2 of reference [37]. 3 CCSD(T)-F12/a5Z value from Table 3 of reference [38]. 4 MP2/haTZ value from Table 2 of reference [37]. 5 The "semi-empirical" value from reference [38].

Figure 1.
The fragment of ε-glycine in its solid-phase structure [19] that is discussed in the text.
The presented procedure for a fully reliable estimation of the canonical CCSD(T)/CBS ∆ was applied to two dimers clipped out from the ND structure of the ε polymorph of glycine [19], which are pictured in Figure 1. The distances between atoms marked H and O2, and H and O1, amount to 218 and 220 pm, respectively. Based on these distances, the C-H•••O2 contact in the dimer shown in the left in Figure 1 would appear to be even stronger than a similar one in the other dimer, namely, the C-H•••O1 contact that is present in the dimer appearing in the right side of Figure 1. However, as already shown in reference [19], intermolecular interactions are destabilizing in the arrangement that includes the C-H•••O2 contact and are stabilizing in the dimer with the C-H•••O1 contact, which has a longer distance between pertinent proton and oxygen atoms. The ∆ values obtained for these two complexes are +51.8 and −41.1 kJ/mol, respectively, while it is noted that the PIXEL method [39] provided highly accurate values amounting to +51.2 and −40.4 kJ/mol, respectively (see Table 2 of reference [19], where also other geometrical parameters of the C-H•••O motifs are listed). It is also noted that the ∆ (DFT) and ∆ (DLPNO) results are +52.6 and +53.7 kJ/mol, respectively, for the repulsive arrangement (see references [40] and [41] for a discussion of repulsive intermolecular interactions), and they amount to −39.2 and −41.1 kJ/mol, respectively, in the stable dimer. It thus can be seen that all computational methods provided very similar interaction energies. Overall, these calculations suggest that the same structural environment needs to be included for a systematic comparison of the strength of various types of HBs. In the following part of the paper, such comparison is presented for the intermolecular C-H•••O and O-H•••O H-bonding in a simple model, namely, the AA dimers. At this point, it should be mentioned that the Supplementary Materials files "dimers.xlsx" and "DLPNO.xlsx" contain underlying data for the ∆ and ∆ (DLPNO), respectively, estimation of complexes described so far (benzofuran:formaldehyde and dimers of formamide and ε-glycine). The fragment of ε-glycine in its solid-phase structure [19] that is discussed in the text.
The presented procedure for a fully reliable estimation of the canonical CCSD(T)/CBS ∆E was applied to two dimers clipped out from the ND structure of the ε polymorph of glycine [19], which are pictured in Figure 1. The distances between atoms marked H and O2, and H and O1, amount to 218 and 220 pm, respectively. Based on these distances, the C-H· · · O2 contact in the dimer shown in the left in Figure 1 would appear to be even stronger than a similar one in the other dimer, namely, the C-H· · · O1 contact that is present in the dimer appearing in the right side of Figure 1. However, as already shown in reference [19], intermolecular interactions are destabilizing in the arrangement that includes the C-H· · · O2 contact and are stabilizing in the dimer with the C-H· · · O1 contact, which has a longer distance between pertinent proton and oxygen atoms. The ∆E values obtained for these two complexes are +51.8 and −41.1 kJ/mol, respectively, while it is noted that the PIXEL method [39] provided highly accurate values amounting to +51.2 and −40.4 kJ/mol, respectively (see Table 2 of reference [19], where also other geometrical parameters of the C-H· · · O motifs are listed). It is also noted that the ∆E(DFT) and ∆E(DLPNO) results are +52.6 and +53.7 kJ/mol, respectively, for the repulsive arrangement (see references [40,41] for a discussion of repulsive intermolecular interactions), and they amount to −39.2 and −41.1 kJ/mol, respectively, in the stable dimer. It thus can be seen that all computational methods provided very similar interaction energies. Overall, these calculations suggest that the same structural environment needs to be included for a systematic comparison of the strength of various types of HBs. In the following part of the paper, such comparison is presented for the intermolecular C-H· · · O and O-H· · · O H-bonding in a simple model, namely, the AA dimers. At this point, it should be mentioned that the Supplementary Materials files "dimers.xlsx" and "DLPNO.xlsx" contain underlying data for the ∆E and ∆E(DLPNO), respectively, estimation of complexes described so far (benzofuran:formaldehyde and dimers of formamide and ε-glycine).

Model Geometries
The acrylic acid dimer (AAD) is one of systems for which the double proton transfer was characterized experimentally (see the review [42]). In particular, the group of Caminati applied microwave spectroscopy measurements to the polar form of AAD and estimated a barrier for the exchange of its hydroxyl protons [43]. Here, nonpolar forms of AAD are considered. Namely, the model for an investigation of C-H· · · O H-bonding has both AA units with cis configuration of hydroxyl group relative to vinyl group (see Figure 2a), while that orientation is trans in the model describing O-H· · · O interactions ( Figure 2b). The full geometry optimization for these two models was carried out at the MP2/aTZ level as in the previous work on molecular complexes [36]. The distance, R, between atoms marked in Figure 2a with "C" and "O", and in Figure 2b with "O d " and "O a ", is 341 and 264 pm, respectively, in the pertinent MP2/aTZ minima. All other geometric parameters of these minima can be inferred from their coordinates provided in "min" sheets inside the Supplementary Materials files "CHO.xlsx" and "OHO.xlsx". It should be noted that in the single-crystal X-ray diffraction structure of AA that features layers of dimers with O-H· · · O hydrogen bonding, the corresponding value of R is 265 pm [44]. The two MP2/aTZ equilibrium geometries were used to vary the parameter R while fixing all other degrees of freedom in order to follow the PES curves, which are described in the next paragraph.

Model Geometries
The acrylic acid dimer (AAD) is one of systems for which the double proton transfer was characterized experimentally (see the review [42]). In particular, the group of Caminati applied microwave spectroscopy measurements to the polar form of AAD and estimated a barrier for the exchange of its hydroxyl protons [43]. Here, nonpolar forms of AAD are considered. Namely, the model for an investigation of C-H•••O H-bonding has both AA units with cis configuration of hydroxyl group relative to vinyl group (see Figure  2a), while that orientation is trans in the model describing O-H•••O interactions ( Figure  2b). The full geometry optimization for these two models was carried out at the MP2/aTZ level as in the previous work on molecular complexes [36]. The distance, , between atoms marked in Figure 2a with "C" and "O", and in Figure 2b with "Od" and "Oa", is 341 and 264 pm, respectively, in the pertinent MP2/aTZ minima. All other geometric parameters of these minima can be inferred from their coordinates provided in "min" sheets inside the Supplementary Materials files "CHO.xlsx" and "OHO.xlsx". It should be noted that in the single-crystal X-ray diffraction structure of AA that features layers of dimers with O-H•••O hydrogen bonding, the corresponding value of is 265 pm [44]. The two MP2/aTZ equilibrium geometries were used to vary the parameter while fixing all other degrees of freedom in order to follow the PES curves, which are described in the next paragraph.

The Interaction Energy Curves
For each type of AAD models introduced in Section 2.2, the CCSD(T)/CBS interaction energies were obtained at 20 points in a wide range of the respective intermonomer separation as expressed by the parameter . There are two equivalent hydrogen bonds in every investigated geometry due to the C2h symmetry of both types of AA dimers.  Table S1, and the actual geometries and raw energies in "CHO.xlsx" and "OHO.xlsx" spreadsheets). The two sets of ∆ values were accurately fitted to the same functional form, which is specified in Section 4 (further details are given in Table S2). Then, minima were found at distances

The Interaction Energy Curves
For each type of AAD models introduced in Section 2.2, the CCSD(T)/CBS interaction energies were obtained at 20 points in a wide range of the respective intermonomer separation as expressed by the parameter R. There are two equivalent hydrogen bonds in every investigated geometry due to the C 2h symmetry of both types of AA dimers. Figures 3 and 4 graphically present the results for C-H· · · O and O-H· · · O interactions, respectively (all R and ∆E values are collected in the Supplementary Materials Table S1, and the actual geometries and raw energies in "CHO.xlsx" and "OHO.xlsx" spreadsheets). The two sets of ∆E values were accurately fitted to the same functional form, which is specified in Section 4 (further details are given in Table S2). Then, minima were found at distances denoted as R min. . Values of R min. amount to 339 and 269 pm for C-H· · · O and O-H· · · O models, respectively, and the interaction energy at these points is −19.6 and −86.2 kJ/mol, respectively. These R min. values are thus in a good agreement with the aforementioned results of the unconstrained MP2/aTZ optimization (R min. = 341 and 264 pm accordingly). Moreover, the interaction energies at the minima of relaxed geometries are accordingly −19.9 and −85.6 kJ/mol and thus are quite close to their counterparts obtained from the fits of ∆E data. It is worth noting that the fits also enable to estimate a maximal distance at which the investigated contacts become repulsive, that is, a value of R for which ∆E = 0 kJ/mol. This distance is 294 and 226 pm for C-H· · · O and O-H· · · O interactions, respectively. Inflexion points of the fitted curves were located, too, and the corresponding distance between monomers is denoted as R infl. . Namely, R infl. of 382 and 305 pm, respectively, was found for the curves shown in Figures 3 and 4 (related ∆E amounts to −15.2 and −68.1 kJ/mol, respectively). Positions of these inflexion points are analyzed in terms of the scaled intermolecular distance R/R min. in Section 3 together with other parameters of the investigated PES. It is worth mentioning that results of the ∆E(DFT) and ∆E(DLPNO) calculations are −20.0 and −19.5 kJ/mol, respectively, for the structure fully optimized by the MP2/aTZ method and featuring C-H· · · O interactions. These values are −87.5 and −85.8 kJ/mol, respectively, for the MP2/aTZ optimized geometry with O-H· · · O H-bonding. It should also be mentioned that the accurately fitted dependencies of ∆E upon R might be of interest in a development and/or testing of force fields [45].

SAPT-DFT Partitioning of the Interaction Energy
The SAPT-DFT method is a powerful tool for theoretical investigations of noncovalent bonding [46], while it is noted that also the energy decomposition analysis (EDA; see reference [47]) in the DFT framework can be usefully applied to noncovalent interactions [48,49]. The SAPT-DFT/CBS approach, which is specified in Section 4, enables accurate predictions of the intermolecular interaction energy (denoted as in the following).

SAPT-DFT Partitioning of the Interaction Energy
The SAPT-DFT method is a powerful tool for theoretical investigations of noncovalent bonding [46], while it is noted that also the energy decomposition analysis (EDA; see reference [47]) in the DFT framework can be usefully applied to noncovalent interactions [48,49]. The SAPT-DFT/CBS approach, which is specified in Section 4, enables accurate predictions of the intermolecular interaction energy (denoted as E total in the following). In particular, the root-mean-square deviation between the E total and corresponding CCSD(T)/CBS data (denoted as ∆E(CC) in the following) for a diverse set of 18 dimers was found to be as low as 0.84 kJ/mol [24]. At the same time, the SAPT-DFT technique partitions a E total value into well-defined contributions that stem from electrostatic, exchange, induction, and dispersion interactions [50]. These components of the total interaction energy are denoted here as E elst , E exch , E ind , and E disp , respectively, with E total = E elst + E exch + E ind + E disp (see Section 4). Their dependence upon the distance R in AA dimers was followed (all underlying values are collected in Table S3). Results are presented in Figures 5 and 6, where also pertinent ∆E(CC) data are included for comparison purposes. The key values are listed in Table 3 and reveal that the biggest absolute difference between the reference ∆E(CC) results and their E total counterparts is 2.8 kJ/mol and occurs in the model of O-H· · · O interactions with R = 247 pm (the relative difference of these interaction energies is 3.7%), while the highest relative difference amounts to 9.2% (in the case of the C-H· · · O model with R = 320 pm). The present SAPT-DFT/CBS data can thus be considered to be reliable throughout the investigated ranges of intermonomer separations and are analyzed below. Figures 5 and 6 show that the E elst energy, which arises in the first order of the perturbation theory of intermolecular interactions, is the highest stabilizing (that is, the most negative) contribution to all E total values. The E exch term is also the first-order contribution and surpasses E elst at the two shorter distances in models of both C-H· · · O and O-H· · · O interactions. Consequently, higher-order terms are responsible for binding in these four complexes [51]. At the two longer distances, however, the first-order terms alone would lead to an attractive interaction, as the sums (E elst + E exch ) are negative (see Figures 5 and 6, and also Table S3). In the computational approach adopted here, E disp component combines the second-order dispersion contributions, and E ind adds the second-order induction contributions and an estimate of all higher-order terms [52]. The E ind energy is a bit more important than E disp in AA dimers that represent O-H· · · O Hbonding, and conversely E disp contributes slightly more than E ind to a stabilization of the C-H· · · O models. Nevertheless, adding E ind energies to (E elst + E exch ) sums would result in a stabilization of all investigated dimers. This would not hold for E disp contribution, as it amounts to −64.0 kJ/mol in the O-H· · · O model with R = 247 pm and does not offset (E elst + E exch ) of 125.4 kJ/mol, contrary to E ind of −140.0 kJ/mol obtained in this case (see Table S3). The dispersion-to-electrostatics ratio provided by SAPT calculations is a parameter that is frequently employed for the categorization of noncovalently bound complexes [53]. These ratios obtained from the SAPT-DFT/CBS approach are listed in Table 3. They all fall into the range of values that are typical for HBs [54]. Overall, the electrostatics appear to be the most important contribution to the intermolecular bonding in both types of AA dimers considered here.
For the intermonomer separations considered in the previous paragraph, also the EDA calculations were performed. The implementation in the Amsterdam Modeling Suite [55] of the EDA approach from reference [56] was applied together with the B3LYP [57][58][59] combination of DFT functionals, the aforementioned D3 empirical dispersion correction, and the QZ4P basis set [60] (further details are given in Section 4). Resulting B3LYP-D3/QZ4P data are collected in Table S4 for both investigated H-bonding types. These data demonstrate the expected distance dependences of the interaction energy terms, which are graphically presented in Figure S2.  Figure 3. 2 Minimum of the curve shown in Figure 4.  For the intermonomer separations considered in the previous paragraph, also the EDA calculations were performed. The implementation in the Amsterdam Modeling Suite [55] of the EDA approach from reference [56] was applied together with the B3LYP [57][58][59] combination of DFT functionals, the aforementioned D3 empirical dispersion correction, and the QZ4P basis set [60] (further details are given in Section 4). Resulting B3LYP-D3/QZ4P data are collected in Table S4 for both investigated H-bonding types. These data demonstrate the expected distance dependences of the interaction energy terms, which are graphically presented in Figure S2.

Discussion
The most recent review article [61] stressed the importance of a proper denomination of weak interactions in various kinds of aggregates. For a quantification of multitude of intermolecular contacts, it should be beneficial to apply high-level quantum chemical calculations to models that avoid competing structural effects [62]. This way the intrinsic strength of the respective types of noncovalent bonding can be obtained. Here, the CCSD(T)/CBS approach was used to predict the intermolecular interaction energy as a

Discussion
The most recent review article [61] stressed the importance of a proper denomination of weak interactions in various kinds of aggregates. For a quantification of multitude of intermolecular contacts, it should be beneficial to apply high-level quantum chemical calculations to models that avoid competing structural effects [62]. This way the intrinsic strength of the respective types of noncovalent bonding can be obtained. Here, the CCSD(T)/CBS approach was used to predict the intermolecular interaction energy as a function of the C-H· · · O or O-H· · · O H-bonding distances in analogous environments provided by the topology of acrylic acid dimers (see Sections 2.2 and 2.3). Accurate fits of resulting ∆E to the same functional form were obtained and enable a side-by-side comparison of the two interactions. It is convenient for such comparison to employ the relative distance, r, r = R/R min. with the respective R min. values specified in Section 2.3. Figure 7 depicts the dependence of two sets of interaction energies upon the same r data ranging from about 0.86 up to 1.75. In relative terms, the two curves are very similar (see Figure 7). In particular, inflexion points of the C-H· · · O and O-H· · · O models lie at r = 1.128 and 1.138, respectively. The two ∆E sets, however, span completely different intervals of absolute values in the investigate range of r. For instance, ∆E = −15.0 kJ/mol for the C-H· · · O model is reached at around r = 0.924 in the descending part of the dissociation curve, while at this relative distance the O-H· · · O model features ∆E as high as −76.1 kJ/mol. Additionally, at R min. and R infl. distances, the interaction energy of the O-H· · · O arrangement is over four times higher than exhibited by its C-H· · · O counterpart (see Section 2.3 for details). Consequently, the intrinsic strength of the C-H· · · O interaction is much lower than could be anticipated on the basis of, for example, experimental estimates of interaction energies in fragments of crystalline austdiol (C 12 H 12 O 5 ) featuring either C-H· · · O or O-H· · · O HBs [63], and the ∆E data obtained for N-methylacetamide:dimethylformamide dimers containing either C-H· · · O or classical N-H· · · O HBs [64]. This finding is in line with conclusions of the most recent computational study of host-guest interactions, where the C-H· · · O H-bonding was found to be weak, yet very important for the thermodynamic stabilization of certain complexes [65].  It should be of interest to reproduce the aforementioned finding in some more complex systems. Consequently, using the coordinates of β-maltose from the ND study [66], the cluster was prepared in order to study C-H  It should be of interest to reproduce the aforementioned finding in some more complex systems. Consequently, using the coordinates of β-maltose from the ND study [66], the cluster was prepared in order to study C-H· · · O and O-H· · · O H-bonding separately, but in the same structural environment (see Figure 8). Two dimers are thus considered that both contain one β-maltose and one water molecule. Respective water molecules replace the neighboring β-maltose featuring either C-H· · · O interaction (in this case the C-O distance is 366 ppm) or O-H· · · O interaction (the distance between pertinent oxygen atoms is 277 pm). Due to the size of these dimers (48 atoms [67]. Hence, clusters that contain this water molecule and the pertinent L-asparagine were prepared and the ∆ , ∆ (DLPNO), and ∆ (DFT) values were established while employing the ND geometry. Results are summarized in Table 4 (the underlying absolute energies are available from related Excel spreadsheets included in the Supplementary Materials). As expected, the interaction energy values computed by the three high-level methods agree with each for both types of H-bonding. Data in Table 4     Furthermore, the present theoretical approach was used in order to compare C-H· · · O and O-H· · · O H-bonding arrangements that involve a crystalline water molecule. The ND structure of L-asparagine monohydrate [67] is suitable for such a comparison. As shown in Figure 9, both C-H· · · O and O-H· · · O HBs are formed by the water that bridges L-asparagine molecules, which are positioned approximately in a direction of the crystallographic a axis of the structure from reference [67]. Hence, clusters that contain this water molecule and the pertinent L-asparagine were prepared and the ∆E, ∆E(DLPNO), and ∆E(DFT) values were established while employing the ND geometry. Results are summarized in Table 4 (the underlying absolute energies are available from related Excel spreadsheets included in the Supplementary Materials). As expected, the interaction energy values computed by the three high-level methods agree with each for both types of H-bonding. Data in Table 4 reconfirm the finding presented above, as they show the O-H· · · O interaction to be about four times stronger than its C-H· · · O counterpart in essentially the same structural environment.  [67] is suitable for such a comparison. As shown in Figure 9, both C-H•••O and O-H•••O HBs are formed by the water that bridges L-asparagine molecules, which are positioned approximately in a direction of the crystallographic a axis of the structure from reference [67]. Hence, clusters that contain this water molecule and the pertinent L-asparagine were prepared and the ∆ , ∆ (DLPNO), and ∆ (DFT) values were established while employing the ND geometry. Results are summarized in Table 4 (the underlying absolute energies are available from related Excel spreadsheets included in the Supplementary Materials). As expected, the interaction energy values computed by the three high-level methods agree with each for both types of H-bonding. Data in Table 4

Materials and Methods
All the MP2/aTZ geometry optimizations and estimations of vibrational zero-point energies were performed using the Gaussian 16, revision C.01 suite of codes [68] with default settings. This version of the Gaussian program package was also used to obtain the B2PLYP-D3(BJ)/def2-QZVPD interaction energies.
The CCSD(T)/CBS interaction energies were obtained by the focal-point method expressed by Equation (1) (see reference [24] for further details): where subscripts denote the respective energy terms, namely, the total Hartree-Fock energy (HF), the MP2 correlation energy (MP2), and the higher-order correlation energy (post-MP2), and superscripts specify the basis set used to compute the respective term. The MP2/a5Z correlation energies were obtained in the resolution-of-the-identity integral approximation [69,70] while using the relevant auxiliary basis sets [70]. Calculations of the HF/a5Z and MP2/a5Z energies were performed in Turbomole, version 7.1 [71].
Calculations of the canonical CCSD(T)/aTZ and MP2/aTZ correlation energies were carried out in Molpro 2021.2 [72]. Additionally, in Molpro 2021.2, the SAPT-DFT/CBS interaction energies were estimated using the same procedures as in our most recent work [36]. The E elst , E exch , E disp , and E ind contributions to the total interaction energy, E total , from Section 2.4 are related to the underlying interaction energy terms as follows: E elst and E exch are the polarization and exchange energy contributions, respectively, arising in the first order of the perturbation theory of intermolecular interactions [73]; E disp is the dispersion energy contribution obtained as a sum of the second order terms E SAPT (2) disp. and E SAPT (2) disp.−exch. [74]; and E ind is the induction energy contribution approximated by a sum of the second order terms E SAPT (2) ind.−exch. and E SAPT (2) ind. [75] and the correction term E SAPT δ(HF) , which is computed at the HF level [51]. DLPNO-CCSD(T)/CBS interaction energies were estimated by the focal-point method from reference [24], which applies Equation (2) (the notation is as in Equation (1), and right arrow is used to indicate an application of the two-point extrapolation formula from reference [76]): while the CCSD(T) and MP2 correlation energies were obtained in the DLPNO approximation [77][78][79][80]. The ORCA 5.0.3 program package [81] was used with the "TightPNO" set of parameters for the truncation of the electron-correlation space and with the default method of the orbital localization. The least-squares fits of ∆E(R) data employed the following functional form: ∆E(R; r e , a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , V e ) = a 0 ξ 2 1 + a 1 ξ + a 2 ξ 2 + a 3 ξ 3 + a 4 ξ 4 + a 5 ξ 5 + a 6 ξ 6 + V e with ξ = (R − r e )/R and R as defined in Section 2.2. The Levenberg-Marquardt algorithm from the "lsqcurvefit" function of MATLAB ® Optimization Toolbox™ was applied.
The B3LYP-D3/QZ4P EDA calculations were performed in the AMS version 2022.103 [55]. In these calculations, core orbitals were not frozen, and the "Good" convergence thresholds were applied.

Conclusions
More than 20 years ago, it was found that in a certain context the C-H· · · O H-bonding could be even stronger than the O-H· · · O interaction [82]. However, a computational procedure proposed here, which is based on an analysis of the CCSD(T)/CBS data obtained for accurate geometries, shows that C-H· · · O HBs are intrinsically weaker than O-H· · · O HBs by a factor of about four. This procedure should be applied to systematically study other noncovalent interactions to assess their role in various types of intermolecular association processes [83,84].