*2.2. Aromaticity*

Aromaticity is a multi-factorial concept which is thought to modify and even to determine the structural, energetic, electronic, and magnetic properties of some molecules. Due to lack of an inherent Dirac observable defining it, aromaticity is usually described in terms of its effects on conjugated systems, such as enhanced thermodynamic stability or structural rigidity. Although the idea of aromaticity was conceived solely upon the interpretation of experimental results, the birth of quantum and computational chemistry motivated the development of multiple tools and techniques aiming at its quantitative analysis. One of the most common approaches to study and measure aromaticity is the nucleus-independent chemical shifts (NICS) [46,47] method, as originally proposed in 1996 by Schleyer and coworkers [46]. The NICS approach has been used for several decades to study numerous *π* skeletons in a variety of fields. Nevertheless, some results obtained through the NICS descriptor have turned to be highly questionable [48–50], even contradicting, in some cases, other aromaticity measures based on reactivity [51]. Among many other criteria exploited to quantify aromaticity we can find (i) structural indices, which evaluate bond equalisation [52], (ii) energy decomposition analyses which require a reference molecule [53], and (iii) electronic descriptors which evaluate the amount of electron delocalisation among the atoms forming a cyclic structure. The last-mentioned set includes a number of methods that have been developed relying on the partition of the electronic density offered by QTAIM, such as the Para Delocalisation Index (PDI) [54], the Aromatic Fluctuation Index (FLU) [55], or the Multicenter Index (MCI) [56]. The PDI and FLU approaches provide an estimate of the aromaticity of a system in terms of the electron delocalisation within the cyclic skeleton, whereas the MCI method arises from a generalised population analysis leading to a many-atoms bond index. In the present work, we have made use of some electronic delocalisation-based descriptors, such as the MCI or the FLU indices, to quantify the changes in aromaticity and antiaromaticity of each molecule upon the formation of the corresponding dimer. We have chosen the FLU and MCI indicators to account for the changes in aromaticity and antiaromaticity upon interaction of the monomers under consideration given their proven accuracy and reliability [57]. Furthermore, these methods are fully compatible with the rest of the QCT analyses performed in this report. Unfortunately, the PDI method is not applicable to some of the examined systems herein, because it can only be used for six-membered rings. Finally, we also used the IQA partitioning to study in detail the energetic changes accompanying

the generation of the molecular clusters shown in Figure 2, with a particular emphasis on their role in HB formation.

### **3. Computational Details**

The structures of the hydrogen-bonded dimers in Figure 2 were optimised in the gas phase and the resultant approximate wavefunctions and electron densities were afterwards dumped for further analysis. All geometry optimisations were performed using the ORCA quantum chemistry package version 5.0.3 [58] using the PBE0 hybrid functional [59] along with the Def2-TZVP basis set [60] and the atom-pairwise dispersion correction with the Becke–Johnson damping scheme [61,62]. For the sake of computational efficiency, the Resolution of Identity (RI) approximation was used for the Coulomb integrals with the default COSX grid for HF exchange, as implemented in ORCA [58]. On the other hand, the auxiliary Def2 basis set was used for the RI-J approximation. The combination of such an exchange-correlation functional and basis set has proven [32] suitable for the characterisation of the systems under study. Moreover, DNLPO-CCSD(T)/def2-QZVPP single point calculations were performed on the optimised geometries of the monomers and dimers in order to test the accuracy of our DFT results. An extrapolation to the complete basis set limit was performed through the def2-SVP/def2-TZVP scheme as implemented in ORCA [58] so as to ameliorate the Basis Set Superposition Error (BSSE). The nature of the stationary points (corresponding to local minima of the potential energy surface) was characterised through the computation of the corresponding harmonic frequencies. QTAIM and IQA calculations were performed using the AIMALL [63] and PROMOLDEN codes [64]. The exchange-correlation energy was partitioned as indicated in reference [65]. Finally, all aromaticity indices discussed along this manuscript were computed using the ESI-3D code [66]. We have denoted the dimers with (i) the hydrogen bond Acceptor Contained within the Ring and (ii) the hydrogen bond Donor Contained within the Ring, as ACR and DCR, respectively. For the ACR azet-2-1H-one (AZH) dimer, we performed a geometry-constrained optimisation to ensure the attainment of the right tautomer of the constituting monomers. On the other hand we observed for the DCR-AZH dimer a structure considerably deviated from planarity as opposed to the rest of the systems. Additional information such as optimised structures, electronic energies, and a more complete survey of IQA and QTAIM (e.g., IQA energy of different groups as well as other QTAIM descriptors) can be found in the electronic supporting information.
