We now analyze the impact of these two corrections on the overall hydrogen bond network of the material when the large size system (300 molecules of each molecular compound) is simulated. We do it through the combined distribution functions (CDF) where the RDF between atoms forming the H-bond (x-axis) and the angular distribution function that contains these two atoms (y-axis) are shown in combination. Figure shows the CDFs of the H-bonds for which the parameter has been optimized. In terms of occurrence values, these three type of H-bonds emerge as the strongest ones throughout the simulations (see Figure for comparison with the remaining type of H-bonds). Notably, FF1 underestimates the interaction with LACOOH and LAOH and slightly overestimates the interaction with choline (see Figure ). The agreement between AIMD with the FF2 and FF3 force fields regarding the H-bonds shown in Figure is highly satisfactory, despite the large di↵erences on the size of the system and simulation timescale employed at each level of theory. Importantly, the observed similarities extend beyond RDFs, which were the ones employed during the fitting process. The di↵erent ADFs also exhibit large agreement. This supports the use of small size systems containing 16 molecules of each chemical compound for the fitting process.