The appeal of hybrid QM/MM models for catalytic reactions in MOFs lies in their optimal balance between simulation accuracy and computational costs. Choomwattana et al. examined carbonyl-ene reactions between formaldehyde and propylene in MOF-11 using the QM/MM method. Their findings revealed that incorporating the extended framework of MOF significantly lowers the reaction barrier compared to the uncatalyzed system. Similarly, Yadnum et al. applied QM/MM calculations to investigate the Mukaiyama aldol reaction between formaldehyde and silyl ether on Cu-MOF-505 and Cu-ZSM-5. These simulations indicated that the inclusion of the microporous framework stabilizes reacting species, resulting in lower energy barriers compared to those in models excluding the framework. Hirao et al. conducted a series of studies on catalytic mechanisms in various systems using the QM/MM methodology. For example, they explored the asymmetric ringopening reaction of cyclohexane oxide with aniline in a twodimensional MOF containing a Cu-paddlewheel (Cu-PDW) unit. Their results underscored the critical role of CH- interactions between aniline and naphthol groups in dictating enantioselectivity. Additionally, QM/MM simulations were used to investigate the catalytic mechanisms of CO2 reacting with styrene oxide in Co-MOF-74 and Mg-MOF-74, revealing distinct rate-determining steps in these systems. Moreover, the hydroxylation of ethane through a high-spin rebound mechanism in Fe-MOF-74 was studied using QM/MM calculations, showcasing the broad applicability and efficacy of the method. However, accurate modeling of interactions between QM and MM regions is pivotal for maintaining the efficacy of QM/MM simulations. In QM/MM, the electrostatic interactions between the QM and MM regions can be modeled by mechanical or electrostatic embedding schemes. The mechanical embedding approach, while commonly employed in numerous studies, overlooks the polarization of the QM region by the MM atoms, as it calculates interactions between the QM and MM regions at the MM level. This simplification may not adequately capture the complex interplay within MOFs, as demonstrated by Wu et al., who highlighted the crucial role of polarization effects in accurately depicting the catalytic behavior of MOF catalysts. They advocate for electrostatic embedding, where electrostatic interactions are treated at the QM level, to better reflect the complex interactions between the QM and MM regions. Moreover, the accuracy of QM/MM calculations significantly hinges on the MM parameters employed. Despite the availability of several force fields designed for MOFs, the selection and validation of MM parameters for QM/MM simulations demand rigorous scrutiny. Parameters optimized for classical force field simulations may not seamlessly translate to the hybrid computational framework of QM/MM methods. Currently, there appears to be a lack of an electrostatic embedding QM/MM model specifically for Zr-MOFs, a class of materials renowned for their utility in gas storage, separations, and catalysis. This gap highlights a significant opportunity for advancing the field by developing models that can more accurately simulate adsorptions and catalytic reactions occurring in Zr-MOFs.