Predicting the phase diagram of solid solutions is crucial for understanding their physical behavior under various conditions and at different compositions. However, conventional sampling methods face challenges in efficiently addressing configurational and vibrational disorder for substitutional solid solutions. Additionally, a persistent obstacle has been the lack of a robust theoretical approach for determining the free energy of interstitial solid solutions characterized by the fluid-like diffusion of interstitial species. The method presented in this paper overcomes both hurdles by coupling thermodynamic integration with hybrid Monte Carlo algorithms. We validate the accuracy of the method by computing the free energies of iron alloys described by a Lennard-Jones potential. We also showcase its efficiency by determining the phase diagram of the MgO-CaO system, described by a machine-learning interatomic potential. The high efficiencies achieved with this method pave the way to the determination of the free energies of solid solutions with ab initio accuracy.