The accurate determination of stability constants for metal-ligand complexes is central to understanding metal speciation in aqueous environments, yet it remains a significant challenge for theoretical modelling. Density-functional theory (DFT) calculations often struggle to produce reliable results for these species due to the complexities of modelling the solvation of charged metal complexes. In this work, we compare a purely implicit solvation model with a hybrid implicit/explicit approach to calculate the solvation energies of a series of metal cations, as well as the stability constants of their corresponding metal-cyanide and metal-glycinate complexes. By employing thermodynamic cycles to incorporate discrete water molecules into the first coordination sphere, we achieve a more physically realistic description of short-range solute-solvent interactions. While purely implicit models produce significant errors in solvation energies, the hybrid approach substantially improves accuracy, leading to calculated stability constants that showed strong linear correlations with experimental data across different metal ligand systems. This refined methodology provides a robust predictive framework for determining stability constants in cases where experimental data are unavailable and offers a practical pathway for computational screening of metal-ligand systems in aqueous environments.