The most important experimental quantity describing the thermodynamics of metal-ion binding with various (in)organic ligands, or biomolecules, is the stability constant of the complex (β). In principle, it can be calculated as the free-energy change associated with the metal-ion complexation, i.e., its uptake from the solution under standard conditions. Because this process is associated with the interactions of charged species, large values of interaction and solvation energies are in general involved. Using the standard thermodynamic cycle (in vacuo complexation and solvation/desolvation of the reference state and of the resulting complexes), one usually subtracts values of several hundreds of kilocalories per mole to obtain final results on the order of units or tens of kilocalories per mole. In this work, we use density functional theory and Møller-Plesset second-order perturbation theory calculations together with the conductor-like screening model for realistic solvation to calculate the stability constants of selected complexes--[M(NH3)4](2+), [M(NH3)4(H2O)2](2+), [M(Imi)(H2O)5](2+), [M(H2O)3(His)](+), [M(H2O)4(Cys)], [M(H2O)3(Cys)], [M(CH3COO)(H2O)3](+), [M(CH3COO)(H2O)5](+), [M(SCH2COO)2](2-)--with eight divalent metal ions (Mn(2+), Fe(2+), Co(2+), Ni(2+), Cu(2+), Zn(2+), Cd(2+), and Hg(2+)). Using the currently available computational protocols, we show that it is possible to achieve a relative accuracy of 2-4 kcal·mol(-1) (1-3 orders of magnitude in β). However, because most of the computed values are affected by metal- and ligand-dependent systematic shifts, the accuracy of the "absolute" (uncorrected) values is generally lower. For metal-dependent systematic shifts, we propose the specific values to be used for the given metal ion and current protocol. At the same time, we argue that ligand-dependent shifts (which cannot be easily removed) do not influence the metal-ion selectivity of the particular site, and therefore it can be computed to within 2 kcal·mol(-1) average accuracy. Finally, a critical discussion is presented that aims at potential caveats that one may encounter in theoretical predictions of the stability constants and highlights the perspective that theoretical calculations may become both competitive and complementary tools to experimental measurements.