All QM/MM calculations were carried out using LICHEM , for interfacing Gaussian16 67 and TINKER . The MM region of both the systems were represented with the ff14SB forcefield. The QM region of the Mg 2+ system, consisting of the two metals, incoming nucleotide, terminal 3', conserved carboxylates and coordinating waters of both metals and tri-phosphates, were treated at the ωB97X-D/6-31G(d,p) level. For the Mn 2+ system, the QM region (similar atoms to Mg 2+ system) was treated at the ωB97X-D level with LANL2DZ for the Mn 2+ ions and 6-31G(d,p) for all other atoms, for a preliminary optimization for the sake of efficiency. Once converged, further optimizations were performed with ωB97X-D/6-31G(d,p) for all QM atoms, including Mn 2+ ions. The metals in the active site for the Mn 2+ system were modeled in a high spin state (ferromagnetic, multiplicity of 11), based on previous studies on similar systems . The long-range electrostatic correction (LREC) with a 12 Å cutoff was used for the QM calculations, whereas PME was used for the MM calculations. For proper treatment of the covalent boundaries between the QM and MM atoms, the pseudobond approach was used . The number of QM atoms for Mg 2+ system is 174, with 3536 active MM atoms, whereas for the Mn 2+ system, the number of QM atoms is 171, with 3579 active MM atoms.