Defect formation energies in materials generally depend on chemical potentials determined by a chemical equilibrium condition. In particular, an aqueous solution environment is important for biomaterials such as hydroxyapatite studied here. Therefore, a methodology to obtain ionic chemical potentials under chemical equilibrium between solid and aqueous solution was introduced, and was applied to substitutional divalent cations formed via ion exchange with Ca(2+) in hydroxyapatite. The calculated ranking of the stability of substitutional cations in HAp was in good agreement with the experimentally observed trend. The present theoretical approach would be useful to explore the thermodynamic stability of defects in materials subjected to an aqueous solution environment.