In the above equation, all the harmonic terms in equation ( ) vanished exploiting the FIRA and the remaining flexible torsions term E f t is computed over the same δ angle employed at QM level, thus limiting the mismatch due to the presence of the non-bonded intra-molecular term, E nb . In absence of strong coupling terms between the selected ICs, the latter term is not required and Route I allows for assigning the whole QM relaxed energy sampled along the scans ∆U g to the flexible dihedral term E f t (δ g ). Nonetheless, when dealing with large and very flexible molecules, where dihedral rotations might cause either specific internal interactions, as for instance intra-molecular HBs, or close contacts between far-lying atoms, it becomes necessary to include specific intramolecular non-bonded terms. This introduces a possible geometrical mismatch. One of the main consequences of such a mismatch, is the general underestimation of the r so f t distances evaluated at the unrelaxed g ′ FIRA geometry, which in turn leads to a significant overestimation of the intra-molecular LJ contribution in equation (11). To circumvent this issue, when including LJ intra-molecular terms, in previous work, Route I consisted in a rather time-consuming user-defined inclusion and tuning of the selected interacting atom pairs and their r so f t distances. This approach required the manual inclusion of specific pairs and their corresponding r so f t distances after the initial parameterization. This process hence relied on a trial-and-error methodology, where yet adjustments were still guided by comparing the modified QMD-FF with the original QM database.