**2. Theoretical Methods**

Geometries of all the systems were fully optimized on the *ω*B97X-D/6-311++G(d,p) level of theory, i.e., using the *ω*B97X-D range-separated gradient- and dispersion-corrected hybrid exchange-correlation functional [99] of density functional theory [100–102] and the 6-311++G(d,p) basis set being of the triple-zeta type and containing both polarization and diffuse functions on all atoms [103–107]. It is worth noting that *ω*B97X-D was one of the best functionals out of 200 tested [108]. There were no imaginary frequencies showing that equilibrium structures were obtained each time. Both the geometry optimization and the frequency analysis were performed using the Gaussian 16 (Revision C.01) program [109].

Values of the electron density at the bond critical points (bcp) [110–112] of the C··· H-D hydrogen bonds and other accompanying interactions were computed using the AIMAll program [113]. Indeed, the C··· H-D hydrogen bond should also be investigated through quantum chemical topology, which, over the years, was shown to be an efficient theoretical approach [114,115].

The dissociation energy was calculated as the difference between the total energy of a dimer and the sum of total energies of isolated subsystems in their own fully optimized structures. Total energies were corrected for the zero-point vibrational energies (ZPVE). Dissociation energies are given as positive values. In order to determine the binding energy (in fact, the interaction energy [116]) of the individual interaction of interest, the formula suggested recently by Emamian et al. [116] was used:

$$E\_{\rm b}[\text{kcal/mol}] = -223.08 \cdot \rho\_{\rm bcp}[\text{a.u.}] + 0.7423 \tag{1}$$

This formula is based on the electron density value determined at the bond critical point of this interaction, which is the parameter that, as has been shown [116], is best correlated with binding energy among many wave function-based HB descriptors. As a result, Emamian et al. proposed using this equation for a quick estimate of the energy of hydrogen-bond-forming networks. It is worth emphasizing that the values of the

dissociation energies and the (sum) of the binding energies determined by Equation (1) are not comparable, since the former is a global quantity relating to the entire dimer, while the latter is a local parameter relating to an individual interaction. Besides, the dissociation energy takes into account deformation energies and ZPVEs.
