Fast Quantum Approach for Evaluating the Energy of Non-Covalent Interactions in Molecular Crystals: The Case Study of Intermolecular H-Bonds in Crystalline Peroxosolvates

Energy/enthalpy of intermolecular hydrogen bonds (H-bonds) in crystals have been calculated in many papers. Most of the theoretical works used non-periodic models. Their applicability for describing intermolecular H-bonds in solids is not obvious since the crystal environment can strongly change H-bond geometry and energy in comparison with non-periodic models. Periodic DFT computations provide a reasonable description of a number of relevant properties of molecular crystals. However, these methods are quite cumbersome and time-consuming compared to non-periodic calculations. Here, we present a fast quantum approach for estimating the energy/enthalpy of intermolecular H-bonds in crystals. It has been tested on a family of crystalline peroxosolvates in which the H∙∙∙O bond set fills evenly (i.e., without significant gaps) the range of H∙∙∙O distances from ~1.5 to ~2.1 Å typical for strong, moderate, and weak H-bonds. Four of these two-component crystals (peroxosolvates of macrocyclic ethers and creatine) were obtained and structurally characterized for the first time. A critical comparison of the approaches for estimating the energy of intermolecular H-bonds in organic crystals is carried out, and various sources of errors are clarified.


Introduction
A quantitative description of intermolecular (non-covalent) interactions in the gas phase can be achieved using a variety of theoretical approaches: NBO analysis [1], SAPT [2], etc. [3][4][5]. There is a large amount of experimental data from which "reference" values of the energy/enthalpy of intermolecular interactions can be obtained [6][7][8][9][10]. Molecular structure, chemical bond theory, and "standard" quantum chemical approaches focus on non-periodic systems. Their "generalization" to molecular crystals is not straightforward because of the need to explicitly take into account the periodic boundary conditions [11,12]. The identification and quantitative description of intermolecular interactions in crystals is a non-trivial task, both experimentally and theoretically. An exception are molecular crystals in which the sublimation enthalpy/lattice energy is determined by one type of intermolecular (non-covalent) interactions, such as CF 4 [13,14]. A quantitative characteristic of intermolecular interactions in a crystal is the lattice energy or sublimation enthalpy, extrapolated to 0 K [15]. A number of theoretical approaches to the lattice energy assessment are reported in the literature. They mostly concern single-component crystals and use either quantum chemical modeling [16][17][18] or semi-empirical schemes [19][20][21]. In most of the papers cited above, short (strong) intermolecular H-bonds [22] are not considered since they are practically not realized in single-component crystals. On the other hand, such H-bonds determine the structure of multicomponent molecular crystals used in pharmaceuticals [23][24][25] and crystal engineering [26]. Intermolecular H-bonds of different types and strengths play a decisive role in the structure of crystal hydrates [27,28] and peroxosolvates [29,30]. Currently, there is a great demand for approaches that make it possible to quantitatively describe various intermolecular (non-covalent) interactions in organic crystals. In most cases, empirical approaches are used that relate the energy of the intermolecular interaction with one or another parameter of the electron density at the bond critical point [31]. In this case, the calculated values of the electron density, the values of the parameters retrieved from the precise X-ray diffraction data, and hybrid approaches are used [32]. The above approaches are often used to evaluate intermolecular interactions of different natures in various crystals [33][34][35], which gives rise to well-founded criticism [36,37]. Possible reasons for this criticism include the following. (1) The identification of a non-covalent interaction is based solely on the presence of the bond critical point [38]. (2) The limits of applicability of these approaches are not clear. For example, the approach developed in ref. [31] was used to estimate energies from~1.5 [33] to~600 kJ/mol [39]. (3) The dependence of the theoretical values of the electron density parameters at the bond critical point on the approximations used in the calculations by the DFT method has not been studied.
Peroxosolvates of organic molecules are convenient objects for a critical comparison of various DFT functionals and empirical approaches to assessing the energy of intermolecular interactions, that is, conventional H-bonds. At first, identification of the conventional H-bonds in crystals may be made without Bader analysis of periodic electron density [40]. Secondly, the enthalpy of H-bonds can be estimated by methods that do not use electron density characteristics, namely, changes in IR intensities [41], frequency shifts [40], and metric parameters of intermolecular H-bonds [42]. Third, organic crystals imply calculations with B3LYP and PBE-D3 functionals, which are widely used in calculations for periodic [12,14,17,18,34,43] and non-periodic [35] structures. Finally, H-bonds of various nature and strengths are realized in peroxosolvates [30].
In the present paper, we focus on the applicability of the B3LYP and PBE-D3 functionals for evaluating the geometric parameters of intermolecular H-bonds of different types and strengths in peroxosolvates of macrocyclic ethers and amino acids (Scheme 1). Four of these two-component crystals were obtained and structurally characterized for the first time. A critical comparison of several methods for estimating the energy of intermolecular H-bonds is carried out, and various sources of errors are clarified. Finally, a fast quantum approach for evaluating the energy of non-covalent interactions in molecular crystals is suggested. This approach is based on non-periodic (cluster) models. Its accuracy and limits of applicability are demonstrated. Scheme 1. Structures of the macrocyclic ethers (1-3) and creatine (4), the peroxosolvates I-IV of which were obtained and studied by X-ray diffraction.

Crystal Structures and Hydrogen Bonding
The structure of 18-crown-6 peroxosolvate C 12 H 24 O 6 •5H 2 O 2 (I) comprises one independent cyclic ether molecule lying on a crystallographic inversion center and three solvent peroxide molecules (Figure 1). Two of them occupy general positions, and the third, namely H8-O8-O8A-H8A, is sited on a crystallographic two-fold axis. The molecule H4-O4-O5-H5 is disordered over two sites with an occupancy ratio of 0.927(3)/0.083(3) (Figure 1). In I, the ether molecule adopts an ordinary crown conformation with O1, O3, and O2A atoms lying on one side of the molecule and O2, O1A, and O3A on the other. The C-O bond lengths vary within the normal range (1.4209(12)-1.4288(11) Å) for 18-crown-6 derivatives (CSD ver. 5.43 November 2021 [44]). All oxygen atoms of organic coformer are involved in one hydrogen bond with neighboring H 2 O 2 molecules. In the solvent molecules, O-O distances lie within the normal range for peroxosolvates (1.441(10)-1.4692(10) Å) [45]. The torsion H-O-O-H angles (80.8(18)-139.4(15) • ) are also typical for solvent peroxide molecules. As in the structures of all known peroxosolvates, H 2 O 2 molecules in I participate in two donor hydrogen bonds [30]. Additionally, two of three peroxide molecules are involved in acceptor H-bonding. In the case of H6-O6-O7-H7 molecules, this hydrogen bonding results in the unprecedented formation of 2-fold cyclic dimers with R 2 2 (6) topology ( Figure S1). Of interest, atom H7 is involved in bifurcated H-bonding, and I represents the second example of a structure containing an H 2 O 2 molecule with a bifurcated donor hydrogen bond [28]. Recently, it was supposed that the usage of organic coformers without donor H-atoms results in the targeted formation of peroxosolvates with infinite 1D-hydrogen peroxide clusters [46,47]. The structure of I is direct evidence of this rule ( Figure 2). Within the chains, central network-forming H 2 O 2 molecules are hydrogen bonded only with adjacent peroxides, while lateral molecules are additionally linked to crown ether coformers by donor HOOH···OR 2 bonds. To the best of our knowledge, this untrivial chain topology was not observed previously either for peroxosolvates or for hydrate structures.   Figure 3). All molecules lie on general positions. In the organic molecule, all bond lengths and angles have ordinary values [44]. The coformer adopts "OO" conformation [49] with all ether oxygen atoms directed toward the center of the molecule. Obviously, this compact conformation prevents ether oxygen atoms from hydrogen bonding and denies the inclusion of solvents inside the cryptand skeleton. All solvent peroxide molecules possess skew conformation  In the crystal, the neighboring molecules are combined in N + -O − ···(HOOH)n··· − O-N + (n = 1, 2) synthons. These fragments were shown to be structure-forming in peroxosolvates of organic N-oxides [50]. Finally, H-bonded layers perpendicular to ab-diagonal are formed in the structure II ( Figure 4; Figure 5). To the best of our knowledge, II is the second example of a crystal structure containing [2.2.2]-cryptand di-N-oxide molecule [51].  The structure of dibenzo-18-crown-6 (db18c6) peroxosolvate C 20 H 24 O 6 •4H 2 O 2 (III) comprises one independent cyclic ether molecule and four solvent peroxide molecules. All molecules lie on general positions ( Figure 6). The metric parameters of the macrocyclic ether are ordinary with bond lengths close to those found in I. The organic coformer adopts an unusual "planar" conformation with approximately parallel aromatic rings (dihedral angle 1.5(4) • ) and with both alkyl oxygen atoms O2 and O5 deviating in opposite directions from the mean plane of the molecule. In contrast, db18c6 possesses a bent conformation in the unsolvated structure [52], which coincides with the most stable ground state conformation in the gaseous phase [53]. According to CSD, a bent conformation with a dihedral angle between the aromatics of greater than 30 • is also preferable for metal and H-bonded complexes of db18c6 (290 of total 298 entries). The remaining eight entries (BOYVED, DIHZEA, NOKVIU, OYAHAZ02, RAXJAB, RORXOL, VADTEZ, and XIGPEJ) correspond to planar conformation with a dihedral angle less than 5 • , most of them possessing a crystallographically imposed central symmetry. Five ether oxygen atoms participate in hydrogen bonds with hydrogen peroxide. Only catecholate atom O3 is not involved in H-bonding. The peroxide molecules possess skew conformation (77(7)-117(7) • ) with O-O bond lengths ranging from 1.427(6) to 1.465(6) Å. All of them participate in two donor hydrogen bonds and one inter-solvent acceptor HOOH···O 2 H 2 linkage. As expected, the latter results in formation of hydrogen-bonded peroxide chains with only lateral peroxide molecules H11-O11-O12-H12 and H41-O41-O42-H42 interacting with the organic coformer ( Figure 7).  The asymmetric unit in the structure of creatine peroxosolvate C 4 H 9 N 3 O 2 •2H 2 O 2 (IV) contains one independent zwitterionic molecule of natural amino acid and two H 2 O 2 molecules ( Figure 8). All geometric parameters of the organic molecule (including conformation) are close to those previously found for the structures of creatine monohydrate and three polymorphs of unsolvated creatine [54][55][56].
In these four structures, the characteristic torsion angle C-N(Me)-CH 2 -C between two planar fragments (guanidinium and amino acid) varies from 68 to 76 • , while in IV, this angle is equal to 81.9(1) • . In IV, all creatine hydrogen atoms are involved in H-bonding ( Figure 8). Of interest, all four sp 2 -hybridized lone electron pairs of the carboxylate group participate in hydrogen bonds with neighboring peroxide molecules resulting in [SA; SA] linkage mode [57].  The crystal structure of phenylserine peroxosolvate (V) ( Figure S5) was discussed in detail previously [58].

H-Bond Enthalpies/Energies in Peroxosolvates of Macrocyclic Ethers and Amino Acids
It was noted that there are several empirical approaches to assessing the enthalpy ∆H HB of the intermolecular H-bond in a solid. The most accurate of these is Iogansen's approach to exploring changes in IR intensities [41]. In contrast to the liquid state [59], the experimental determination of the integral intensities of OH stretching vibrations in crystals is a non-trivial problem. Periodic DFT calculations of IR intensities are usually performed in the "double harmonic approximation" [60]. This has limited applicability for estimating the frequencies and IR intensities of stretching vibrations of OH/NH groups involved in the formation of short (strong) H-bonds due to significant mechanical and electro-optical anharmonicity. A theoretical description of the considered bands beyond the "double harmonic approximation" is a time-consuming procedure [61][62][63]. On the other hand, the Iogansen approach exploring changes in OH stretching frequencies to estimate ∆H HB has limited applicability to ionic H-bonds [30] since some + N-H···O − bonds are not realized in the gas phase. Moreover, the application of this approach to bifurcated H-bonds is not obvious. This is why we use the approaches described below.
The ∆H HB values are estimated using the Rozenberg approach [42]: where the R(H···O) is the H···O distance (nm). Equation (1) describes −∆H HB varying from 10 to 80 kJ/mol [42]. The R(H···O) values are obtained as a result of geometry optimization. The energy of intermolecular H-bonds E HB in the considered crystals is evaluated according to ref. [64] as: where G b is the positively-defined local electronic kinetic energy density at the H···O bond critical point [38]. Equation (2) will be referred to as the G b -E HB approach below. We start with crystal I (crown Figure S1) since it has weak H-bonds (−∆H HB < 25 kJ/mol, R(H···O) > 1.80 Å and ρ b < 0.04 a.u.). Approaches (1) and (2) lead to similar results for these H-bonds (Table S1a). Moderate H-bonds (25 < −∆H HB < 50 kJ/mol, 1.80 > R(H···O) > 1.60 Å and 0.04 < ρ b < 0.06 a.u.) exist in crystalline II (cryptand), III (benzocrown), and creatine diperoxosolvate IV, (Figures S2-S4). The differences between the two approaches can reach 10 kJ/mol (Tables S2a-S4a), and the absolute enthalpy values always turn out to be lower than the calculated values of the energy of the moderate H-bonds. The obtained results agree with the literature data [65][66][67][68].
The −∆H HB and E HB values evaluated for the considered molecular crystals versus the calculated H···O distances are given in Figure 10. With a decrease in the H···O distance (that is, an increase in the strength of H-bonds), the difference between the energy values calculated using Equations (1) and (2) increases. Since the accuracy of the Rozenberg approach for the considered H-bonds is several kJ/mol, we conclude that Equation (2) can significantly overestimate E HB for short (strong) H-bonds in molecular crystals.
The specific feature of Equations (1) and (2) is the possibility of using the G b and R(H···O) values obtained from experiments (high precision X-ray diffraction, neutron diffraction) or the periodic DFT computations followed by the Bader analysis of the periodic electronic density. From a theoretical point of view, Equation (1) presents the simplest case because the R(H···O) values are obtained as a result of geometry optimization. From the experimental point of view, the applicability of this approach is not straightforward due to the usual lack of accuracy that X-ray structures have on the H atoms. It should be noted that neutron diffraction data on H-bonded crystals are very limited [69]. Initially, Equation (2) was developed for gas-phase systems with intermolecular H-bonds with energies ranging from 3 to 90 kJ/mol [31,64]. Currently, approaches that use the characteristics of the electron density at the point (3, −1) are used to estimate the energies of non-covalent interactions in various systems in gas and condensed media, and the E HB values vary within wide limits [33,39]. From a theoretical point of view, errors in the energy values obtained using Equation (2)   To estimate the error in determining the −∆H HB value caused by the use of B3LYP, the R(H···O) values were computed at the PBE-D3/6-31G** level (Tables S6-S8). The −∆H HB values calculated using PBE-D3 agree nicely with those obtained using B3LYP. In accordance with the literature ( [28] and Table S1a in [71]), the PBE-D3 calculations slightly overestimate the H-bonded energy compared to the B3LYP calculations.

Fast Quantum Approach for Evaluating the Energy of Intermolecular H-Bonds in Molecular Crystals
The cluster approximation is widely used to describe the structure and properties of organic crystals: the energy of the intermolecular interactions [72,73], the mobility of charges in organic semiconductors [74,75], and the chemical shift of the nucleus [76]. This requires knowledge of the structure of the considered crystal, that is, the cif file. Usually, a cluster consisting of several molecules is cut out of a crystal, and the geometry is partially optimized, i.e., the position of only hydrogen atoms is varied [77]. In contrast to C-H bonds whose length can be uniformly normalized [78], the length of OH and NH depends on the strength of the H-bonds formed by these groups ( Table 2 in [22]). Note that the complete optimization of the cluster under non-periodic conditions can significantly change the mutual orientation of the molecules, that is, the number and strength of the studied H-bonds [79].
The cluster size is determined by the condition that all "unique" H-bonds that occur in the crystal under study are taken into account. In some cases, it can be quite large. For all the crystals considered, except for crystal III, the cluster includes no more than 100 atoms. Calculations of such clusters are an order of magnitude faster than calculations of crystals (geometry optimization). The total enthalpies/energies of all unique H-bonds in each of the studied crystals, calculated using periodic DFT, are given in Tables SXa in Supplementary Materials, where X = 1, 2, . . . 5. The corresponding values obtained in the cluster approximation are given in Tables SXb in Supplementary Materials, where X = 1, 2, . . . 5. It follows that the total enthalpies/energies calculated using the cluster approximation are in good agreement with the values obtained from periodic DFT calculations. This justifies the use of the cluster approximation for estimating the energies of H-bonds in molecular crystals. The G b -E HB approach gives reasonable results for energies from 2.0 [79] to 50 kJ/mol and is applicable to describe bifurcated H-bonds formed by the OH group. This has been called the fast quantum approach.
It should be emphasized that the analysis of normal vibrations of the cluster structure, obtained as a result of partial optimization usually leads to imaginary frequencies. Therefore, this approximation is hardly applicable for calculating the experimentally observed characteristics of molecular crystals, namely: sublimation enthalpies [17,18,43], IR [80,81] and Raman spectra [68,71].
Unlike Equation (1), the G b -E HB approach can be applied to intramolecular H-bonds [82,83], which play a structure-directed role in pharmaceutical multicomponent organic crystals [84,85]. It should be noted that empirical approaches use different geometric parameters to estimate the energies of inter-and intramolecular H-bonds in organic crystals, compare [42] with [86].

Cocrystal Preparation
Colorless crystals of I-IV were obtained by cooling to −21 • C saturated solutions (rt) of respective compounds in 98% hydrogen peroxide. Compounds I and III are stable just for a few minutes without mother liquor, while crystal IV is stable for several hours. In contrast, the peroxosolvate of cryptand II demonstrates high stability and does not exhibit noticeable decay for several months.

Single Crystal X-ray Diffraction Experiments
The crystals were extracted from the mother liquor and covered immediately with inert oil to prevent contact with atmospheric moisture. Experimental data sets were collected on a Bruker SMART APEX II diffractometer (for I, III, and IV) and a Bruker SMART Photon II diffractometer (for II). Absorption corrections based on measurements of equivalent reflections were applied [90]. The structures were solved by direct methods and refined by full-matrix least-squares on F 2 with anisotropic thermal parameters for all non-hydrogen atoms [91]. The studied thin plate crystal of II was found to be twinned by the law (1-10)  with an approximate domain ratio of 0.79/0.21. Both domains were indexed separately, and full data were integrated and scaled for absorption with two matrices (APEX III software). The data for the main domain were used for calculations. In I, one of the peroxide molecules was disordered over two sites with an occupancy ratio of 0.927(3)/0.083(3). All hydrogen atoms in I and IV were found from difference Fourier synthesis and refined with isotropic thermal parameters (except for a minor part of the disordered H 2 O 2 molecule). In II and III, all carbon H atoms were placed in calculated positions and refined using a riding model, while "active" hydrogen atoms were found from difference Fourier synthesis, and their positional parameters were refined. In all structures, the partial substitutional disorder of hydrogen peroxide by water molecules [92][93][94][95] was not observed since no residual peaks with an intensity more than 0.17 e/A 3 were seen in the hydrogen peroxide molecule regions. Crystal data and details of X-ray analysis are given in Table S5a. CCDC 2158207-2158210 contains the supplementary crystallographic data for I-IV. These data can be obtained free of charge via http://www.ccdc.cam.ac.uk/conts/ retrieving.html (accessed on 18 May 2022) (or from the CCDC, 12 Union Road, Cambridge CB2 1EZ, UK; Fax: +44-1223-336033; E-mail: deposit@ccdc.cam.ac.uk).

Periodic (Solid-State) DFT Computations
B3LYP and PBE-D3 are the most popular functionals in periodic (solid-state) DFT calculations of organic crystals [12,17,18,81]. These functionals provide a grounded tradeoff between the accuracy and the rate of calculations of experimentally observed properties of multicomponent organic crystals [28,43,68,71]. Computations with all-electron Gaussiantype localized orbital basis 6-31G** were conducted using the CRYSTAL17 package [96]. We employed B3LYP [97,98] and PBE [99]. The London dispersion interactions were taken into account by introducing the D3 correction with Becke-Jones damping (PBE-D3) developed by Grimme et al. [100]. The space groups and the unit cell parameters of the considered twocomponent crystals (I, II, III, IV, and V) obtained from the X-ray diffraction experiment were fixed, and the structural relaxations were limited to the positional parameters of the atoms (AtomOnly). The model structure of I was obtained from the crystal structure, which has one crystallographically independent hydrogen peroxide molecule disordered over two sites with an occupancy ratio of 0.927(3)/0.073(3). Only the higher occupancy was considered in our computations. Tolerance on energy controlling the self-consistent field convergence for geometry optimizations and frequencies computations is set to 10 −10 and 10 −11 hartree, respectively. The shrinking factor of the reciprocal space net is set to 3. All structures optimized at the B3LYP/6-31G** level are found to correspond to the minimum point on the potential energy surface. Further details of calculations are given in the Supplementary Materials.

Conclusions
Four new crystalline peroxosolvates of 18-crown-6, [2.2.2]-cryptand, dibenzo-18crown-6, and creatine were crystallized from the solution of the corresponding coformer in 98% hydrogen peroxide. H 2 O 2 participate in the formation of moderate hydrogen bonds with organic coformer and adjacent molecules of hydrogen peroxide.
Hydrogen bond enthalpies/energies in peroxosolvates of macrocyclic ethers and amino acids are evaluated using the approaches based on the H···O distance (the R(H···O)-∆H HB approach) and the local electronic kinetic G b densities at the hydrogen bond critical point (the G b -E HB approach). In the case of weak and moderate hydrogen bonds, both approaches give similar results. For short (strong) hydrogen bonds (−∆H HB > 50 kJ/mol, R(H···O) < 1.60 Å and ρ b > 0.60 a.u.), the differences −∆H HB and E HB can be 15 kJ/mol or more. We conclude that the G b -E HB approach should be used with caution for describing short (strong) hydrogen bonds in organic crystals.
Reasonable agreement between the hydrogen bonds enthalpies/energies obtained using the H···O distances and the local electronic kinetic G b densities at the hydrogen bond critical point in crystalline peroxosolvates indicates that the use of the Bader approach for localization and characterization of intermolecular hydrogen bonds in organic crystals is fully justified.
It is shown that the cluster approximation can be used to evaluate intermolecular interactions using the R(H···O)-∆H HB and G b -E HB approaches. The latter approach gives reasonable results for energies from 2.0 to 50 kJ/mol and is applicable to describe bifurcated hydrogen bonds formed by the OH group. It has been called the fast quantum approach. Due to the use of partial optimization, the cluster approximation is hardly applicable for calculating sublimation enthalpies, IR and Raman spectra, and thermodynamic and elastic properties of molecular crystals.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27134082/s1, References [38,75,[101][102][103][104][105] are cited in Supplementary Materials; Section S1, Details of periodic (solid-state) DFT computations; Section S2, Non-periodic DFT computations; Figure Figure S4; Part of the H-bond network in creatine peroxosolvate IV. Displacement ellipsoids are shown at the 50% probability level. H atoms of macrocyclic ether omitted for clarity. The H-bonds are given by the dashed lines. Symmetry operation: (i) x, 1 − y, −0.5 + z; (ii) x, −1 + y, z; (iii) 0.5 − x, 1.5 − y, 1 − z; (iv) 1 − x, 1 − y, 1 − z; Figure S5: Part of the H-bond network in 3-phenylserineperoxosolvate V. Displacement ellipsoids are shown at the 50% probability level. Only «active» hydrogen atoms are shown. The H-bonds are given by the dashed lines. Symmetry operation: (i) −x, y, −z; (ii) 0.5 − x, 0.5 − y, 0.5 + z; (iii) x, y, −1 + z; Table S1a: B3LYP/6-31G** values of the (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the H···O bond critical point in I. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S1b: The values of (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point obtained using partial B3LYP/6-31G** optimization of the I cluster. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns.; Table S2a: B3LYP/6-31G** values of the electron density, ρ b , and the local electronic kinetic energy density, G b , at the H···O bond critical point in II. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns.; Table S2b: The values of (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point obtained using partial B3LYP/6-31G** optimization of the II cluster. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S3a; B3LYP/6-31G** values of the (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point in III. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S3b: The values of (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point obtained using partial B3LYP/6-31G** optimization of the III cluster. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S4a; B3LYP/6-31G** values of the (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point in creatine peroxosolvate IV. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S4b: The values of (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point obtained using partial B3LYP/6-31G** optimization of the IV cluster. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S5a: B3LYP/6-31G** values of the (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point in 3-phenylserine peroxosolvate V (refcode VILGAB). The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S5b: The values of (H···O) distance, R(H···O), the electron density, ρ b , and the local electronic kinetic energy density, G b , at the O···O bond critical point obtained using partial B3LYP/6-31G** optimization of the V cluster. The H-bond energy E HB and enthalpy ∆H HB evaluated using Equations (2) and (1), respectively, are given in the last two columns; Table S6: Distances between the atoms involved in the formation of intermolecular H-bonds R(O···O) and R(H···O) in I a) , obtained using periodic DFT computations (PBE-D3/6-31G** and B3LYP/6-31G**). Theoretical values of the enthalpy, ∆H HB of intermolecular H-bonds evaluated using Equation (1) are given in the last two columns; Table S7: Distances between the atoms involved in the formation of intermolecular H-bonds R(O···O) and R(H···O) in II, obtained using periodic DFT computations (PBE-D3/6-31G** and B3LYP/6-31G**). Theoretical values of the enthalpy, ∆H HB of intermolecular H-bonds evaluated using Equation (1) are given in the last two columns; Table S8: Distances between the atoms involved in the formation of intermolecular H-bonds R(X···O) and R(H···O) in creatine diperoxosolvate IV, where X = N or O, obtained using periodic DFT computations (PBE-D3/6-31G** and B3LYP/6-31G**). Theoretical values of the enthalpy, ∆H HB of intermolecular H-bonds evaluated using Equation (1) are given in the last two columns; Table S9: X-ray structure determination summary.