After achieving convergence in the DFTB3 WTMD simulations, configurations were sampled from each trajectory following the protocol outlined in the Computational Methods. For each structure, the single-point energies were calculated at the PBE, B3LYP and RI-MP2 levels and used to reconstruct the corresponding FE curves by applying MetaFEP. The resulting profiles are reported in Figure with the reference DFTB3 FE profiles represented as dashed (blue or green) lines, and the PBE, B3LYP and RI-MP2 FE curves depicted as orange, red and magenta solid lines, respectively. The energy barriers associated with PET4 acylation (Figure ) obtained at DFTB3 and PBE levels exhibit remarkable agreement, with maxima positioned at ≈ 1.7 along the HLDA reaction coordinate and reaching heights of ≈ 21 and ≈ 22 kcal/mol, respectively. A slight difference is observed for the RI-MP2 and B3LYP FE barriers, which are nearly overlapping in shape and feature maxima of ≈ 24 kcal/mol at ≈ 1.3 along the CV axis. On the product side, the INTac-enz is formed at a relative FE of ≈ 4.8 kcal/mol according to DFTB3 calculation, while a plateau at ≈ 2 kcal/mol is observed in the PBE curve. Instead, the local minima obtained at RI-MP2 and B3LYP levels are positioned between those observed in the other two curves, with minimal values at ≈ 3 kcal/mol. Similar trends are observed for PEF4 acylation (Figure ), for which a first, sharp FE maximum ≈ 10 kcal/mol high is found at a CV value ≈ 0.7 for the DFTB3 and PBE curves, while the barrier observed at the RI-MP2 and B3LYP levels are broader and ≈ 1 kcal/mol higher. The second FE barrier ≈ 2 kcal/mol high is found at ≈ 2.2 Å along the DFTB3 curve, that is slightly shifted at lower values in the PBE (CV ≈ 2.0) and RI-MP2 and B3LYP (CV ≈ 1.8) profiles. Despite these minimal discrepancies, all calculations support the presence of a metastable intermediate, positioned at approximately the same coordinate along the HLDA CV axis (≈ 1.6-1.7) across all four FE curves. Differences are more pronounced in the region corresponding to the formation of the INTac-enz state, with the PBE, RI-MP2 and B3LYP calculations exhibiting reaction energies of ≈ 1.0 -1.5 kcal/mol, against the ≈ 2.5 kcal/mol estimated by DFTB3, suggesting that the process is slightly less endergonic than how predicted from WTMD simulations. Considering the deacylation stages, the results obtained for PET4 (Figure ) indicate a remarkable consistency in shape, height and position of the FE maxima across each method applied to perform the analysis, with a barrier peaking at ≈ 10-11 kcal/mol at a reaction coordinate of 0.2. However, notable differences are found in the reaction energy of the process. DFTB3 calculations yield minimal value of ≈ -1.7 kcal/mol, which is only ≈ 0.5 kcal/mol lower than that predicted at the PBE level. In contrast, the reaction energies obtained from RI-MP2 and B3LYP calculations are significantly lower, ≈ -5.0 kcal/mol, suggesting that the process is substantially more exergonic. Finally, the FE curves calculated for the PEF4 deacylation (Figure ) show some discrepancies among the four methods. Although each level of theory supports the single-step nature of the process, positioning a ≈ 5 kcal/mol high FE barrier at ≈ 0.24 along the reaction coordinate, the RI-MP2 and B3LYP energy walls show two shoulders, symmetrically positioned at CV values ≈ 0.17 and 0.30 relative to the central peak.