Magnesium (Mg2+) is an indispensable regulatory cation for numerous cellular enzymes and membrane proteins that participate in a variety of biological processes, including cell signaling, glycolysis, cell differentiation, and genome stability. Molecular mechanics (MM) simulations have the potential to provide detailed insights into the mechanisms of these Mg-dependent proteins. However, MM force fields continue to suffer from large inaccuracies when dealing with divalent cations. The inclusion of explicit polarization in MM models, such as in CHARMM Drude and AMOEBA models, has been promising, but despite substantial improvements over non-polarizable models, errors in interactions of proteins with divalent cations remain in excess of 10 kcal/mol. Such errors are present even with ‘nonbonded-fix (NB-fix)’ corrections. Modeling becomes even more problematic when dealing with enzymes like kinases and phosphatases, where Mg2+ ions interact directly with both proteins and nucleotides. Toward addressing this long-standing challenge, we recently presented a recalibration of the AMOEBA model (Delgado et al. J Chem. Info. Model. 2022). Notably, polarization terms of cation-coordinating groups were improved to respond better to the high electric fields present near cations. This systematically reduced errors in interactions of monovalent cations with proteins. Here we show that this revised model, along with our many body NB-fix (MB-NB-fix) corrections, reduces errors in Mg2+-protein interaction energies to 5 kcal/mol. MB-NB-fix corrections also improve predictions of Mg2+-ATP interactions substantially. Errors are computed with respect to estimates from vdW-inclusive density functional theory that we benchmark against "gold standard" CCSD(T) calculations and experiments. We will also present the ramifications of these improvements in local interactions on Mg2+-ATP binding free energies and binding modes in bulk water, and on complexation of Mg2+ and ATP with proteins.