Exploring a metal-involved biochemical process at a molecular level often requires a reliable description of metal properties in aqueous solution by classical nonbonded models. An additional C 4 term for considering ion-induced dipole interactions was previously proposed to supplement the widely used Lennard-Jones 12-6 potential (known as the 12-6-4 LJ-type model) with good accuracy. Here, we demonstrate an alternative to modeling divalent metal cations (M 2+ ) with the traditional 12-6 LJ potential by developing nonbonded point charge models for use with 11 water models: TIP3P, SPC/E, SPC/E b , TIP4P-Ew, TIP4P-D, and TIP4P/2005 and the more recent OPC3, TIP3P-FB, OPC, TIP4P-FB, and a99SB-disp. Our designed models simultaneously reproduce the experimental hydration free energy, ion−oxygen distance, and coordination number in the first hydration shell accurately for most of the metal cations, an accuracy equivalent to that of the complex 12-6-4 LJ-type and double exponential potential models. A systematic comparison with the existing M 2+ models is presented as well in terms of effective ion radii, diffusion constants, water exchange rates, and ion−water interactions. Molecular dynamics simulations of metal substitution in Escherichia coli glyoxalase I variants show the great potential of our new models for metalloproteins.