To bridge this gap, we have developed a set of MM parameters tailored for electrostatic embedding QM/MM simulations of Zr-MOFs. This novel parameter set accurately reproduces the adsorption energies of a diverse array of molecules, achieving a root mean square error (RMSE) of 1.1 kcal/mol and a mean absolute error (MAE) of 1.0 kcal/mol. Moreover, it demonstrates robust compatibility across multiple functionals (ωB97X-D, ωB97X-V, B97-D, and M062X), with the MAE maintained within 2 kcal/mol. To showcase the practical applications of our QM/MM model, we present glucose isomerization and epimerization in UiO-66 and MOF-808 as case studies, highlighting the capacity of the model for simulating catalytic reactions within Zr-MOFs. Zr-MOF geometries. The geometries of the Zr-MOFs were modeled using cluster representations based on the crystallographic structures of UiO-66 and MOF-808. As shown in Fig. , both UiO-66 and MOF-808 are represented by cluster models consisting of 31 secondary building units (SBUs). These SBUs are interconnected by 1,4-benzenedicarboxylate (BDC) linkers in UiO-66 and 1,3,5-benzenetricarboxylate (BTC) linkers in MOF-808. The choice of 31 SBUs for each model is driven by the dual goals of ensuring convergence in the calculated adsorption energies and adhering to the principle of centering the QM region within the model. To further refine the UiO-66 model for catalytic activity studies, BDC linkers were Please do not adjust margins Please do not adjust margins strategically removed to introduce defects. This modification is based on evidence from prior research indicating that defective UiO-66 exhibits enhanced catalytic performance. QM/MM computations. Implementation of the QM/MM model in this work followed the electrostatically embedded scheme proposed previously, the details of which can be found elsewhere. As shown in Fig. , adsorbates and an SBU encompassing the active site were described by QM, whereas the other atoms in the MOF structure were described by MM with a standard force field. At the boundary of the QM cluster, a single hydrogen atom was used to replace a C-C bond to divide the QM and MM subsystems. These hydrogen link atoms were positioned along each of the terminal C-C bonds (ca. 154 pm) at a distance of 0.709 R(C-C), which is close to the C-H bond length (ca. 110 pm). All geometry optimizations were conducted with the relaxation of only the QM atoms, keeping all the MM atoms frozen at their crystallographic positions. Because the MM atoms are fixed in space, the only relevant terms in the force field are the two MM-QM interactions, i.e., the electrostatic force and the Lenard-Jones interaction. The electrostatic interaction is given by