Following the equilibration simulations, a series of 12 ns production runs were performed for each 𝜆 window under isothermal-isobaric ensemble conditions using the Parrinello-Rahman pressure coupling scheme. The simulations were performed with and without enhanced sampling using the Hamiltonian replica-exchange MD (H-REMD) approach. All simulations were performed in triplicate to estimate the uncertainty in the calculations and to check for convergence in the results. For each protein-ligand model system, an aggregate total of ≃ 5.5 µs in sampling time was achieved. The results of the simulations were analyzed using the alchemical_analysis.py python tool, and were calculated using multistate Bennet acceptance ratio (MBAR) method after discarding the first 2 ns from each 𝜆 window as equilibration. The reported binding free energy results are the averages and standard deviation of the three independent repeats. ONIOM(QM:MM) Methodology: To estimate the energy leading to the formation of the covalent bond between the target Cys145 of the M pro and the reactive center of the inhibitors, we performed QM/MM calculations on the complexes using the ONIOM 59 approach implemented in Gaussian 16. Explicit solvent water molecules were excluded from the model system. For the QM/MM ONIOM model, we partitioned the system into two layers: (1) a High-level QM region consisting of the full inhibitor structure, Cys145 and His41 side chains (2) Low-level molecular mechanical region which accounts for the remaining residues of the system. On average, our high-level QM region consisted of 74-82 atoms, while the low-level MM region contained 4599-4649 atoms. We treated the high-level region using the M06-2X functional with the def2-TZVP basis set-which performs well for modeling the thermochemistry and reaction energies of chemical reactions. The lowlevel region on the other hand, was treated with the amber molecular mechanical force field. Missing ligand atom parameters were obtained from the generalized Amber force field (GAFF). Hydrogen link atoms were used to cap the QM-MM boundary region. Figure shows a representative snapshot of the atoms comprising the QM region for a selected model system. The ONIOM model structures were geometry optimized using the M06-2X/def2-TZVP:AMBER method within an electronic embedding scheme. We also performed frequency analyses on the optimized geometries using the same level of theory to verify that the optimized structures were minima and to calculate the Gibbs energy.