are both evaluated at the equilibrium geometry (g 0 ). In the second term, which is introduced only in the case of molecules bearing flexible dihedrals δ , the weights W g are all initially set to the same value, ∆U g is the QM computed energy difference between the g and g 0 geometries and N geom is the number of the different geometrical arrangements sampled in the torsional scans for each considered δ in equation (9). Concretely, the JOYCE parameterization was usually carried out minimizing such an objective function through a two-step procedure, hereafter named Route I: i) a first step necessary to determine all harmonic parameters at once using the QM Hessian and minimizing the first term on the right hand side of the objective function (12) ii) only when flexible dihedrals are present in the target molecule, a second step consisting in constraining all harmonic parameters depending on r sti f f (i.e., E s , E b and E st in equation ( )), and parametrizing only the terms related to flexible dihedrals in E f t . These terms use reference QM torsional relaxed energy scans and the algorithm exploits the so-called Frozen Internal Rotation Approximation (FIRA). At this point, it is important to notice that in Route I, the FIRA in step (ii), inevitably introduces a geometry mismatch in equation (12). In fact, for each geometry g considered in the relaxed scans, the the QM and QMD-FF descriptors can be respectively written as: