Bounding region from a set of chemical potentials

Is there any way currently to get some representation of a stability region for variable chemical potentials? (specifically thinking about the case for defect formation energy)
Example Ca in GaN needs Ga to form and Ca3N2 to not form:

3\mu_{Ca} 2\mu_N + < H(Ca3N2)
\mu_{\rm Ga} + \mu_N > H(Al2O3)

These inequalities along with the fact that the chemical potential is negative should allow me to create a bounding region that defines chemical possible chemical potentials.

I think this can be done with something like cddlib but I just want to be sure this is not already implemented somewhere.

There is an implementation of this in pymatgen in the PhaseDiagram.getmu_vertices_stability_phase() function:

https://pymatgen.org/pymatgen.analysis.phase_diagram#pymatgen.analysis.phase_diagram.BasePhaseDiagram.getmu_vertices_stability_phase

In the future, there will likely be some additional features added to pymatgen related to this – particularly the implementation of a more robust algorithm for calculating chemical potential domains, as well as the plotting of chemical potential diagrams in 3D. Please contact me if you (or anyone reading this) needs these features before then.