To address fundamental questions in bioinorganic chemistry, such as metal ion selectivity, accurate computational protocols for both the gas-phase association of metal ligand complexes and solvation/desolvation energies of the species involved are needed. In this work, we attempt to critically evaluate the performance of the ab initio and DFT electronic structure methods available and recent solvation models in calculations of the energetics associated with metal ion complexation. On the example of five model complexes ([M(II)(CH(3)S)(H(2)O)](+), [M(II)(H(2)O)(2-) (H(2)S)(NH(3))](2+), [M(II)(CH(3)S)(NH(3))(H(2)O)(CH(3)COO)] [M(II)(H(2)O)(3)(SH) (CH(3)COO)(Im)], [M(II)(H(2)S)(H(2)O)(CH(3)COO)-(PhOH)(Im)(+) in typical coordination geometries) and four metal ions (Fe(2+), Cu(2+), Zn(2+), and Cd(2+); representing open- and closed-shell and the first- and second-row transition metal elements), we provide reference values for the gas-phase complexation energies, as presumably obtained using the CCSD(T)/aug-cc-pVTZ method, and compare them with cheaper methods, such as DFT and RI-MP2, that can be used for large-scale calculations. We also discuss two possible definitions of interaction energies underlying the theoretically predicted metal-ion selectivity and the effect of geometry optimization on these values. Finally, popular solvation models, such as COSMO-RS and SMD, are used to demonstrate whether quantum chemical calculations can provide the overall free enthalpy (Delta G) changes in the range of the expected experimental values for the model complexes or match the experimental stability constants in the case of three complexes for which the experimental data exist. The data presented highlight several intricacies in the theoretical predictions of the experimental stability constants: the covalent character of some metal ligand bonds (e.g., Cu(II) thiolate) causing larger errors in the gas-phase complexation energies, inaccuracies in the treatment of solvation of the charged species, and difficulties in the definition of the reference state for Jahn-Teller unstable systems (e.g., [Cu(H(2)O)(6)](2+)). Although the agreement between the experimental (as derived from the stability constants) and calculated values is often within 5 kcal.mol(-1), in more complicated cases, it may exceed 15 kcal.mol(-1). Therefore, extreme caution must be exercised in assessing the subtle issues of metal ion selectivity quantitatively.