Calculation of critical chemical potentials

I found the algorithm (pymatgen.phasediagram.analyzer.PDAnalyzer.get_transition_chempots) in pymatgen calculates chemical potential ranges from Gibbs phase diagram but use them to construct grand phase diagram. However, these two kind of phase diagrams seem not to have the same set of equilibrium transitions. How to understand that?


Hi @peikai

Can you give us more info on what you’re doing? Which chemical system? Maybe some code snippets of what you’re comparing?

Hi, @shyamd,

I’m going to find out how to calculate critical potentials in principle. In my opinion, critical potentials are calculated from each of equilibria in Gibbs phase diagrams in the pymatgen with a principle below, with a reference of the method _get_facet_chempots() in class PhaseDiagram, (

In the equilibrium Li-Be-BeO of the Li-Be-O system, we can get \mu_{Li} = -1.909 eV/atom from the equation set,




x_{Be}^{BeO} is the atomic fraction of element Be in phase BeO, that is 0.5

E is the Gibbs free energy per atom (also is the corrected final energy per atom)

In this way, we can get \mu_{Li} for all equilibria of the phase diagram

equilibrium \mu_{Li} eV/atom
Li-BeO-Be -1.909
LiO8-BeO-O2 -5.559
Li2O-Li-BeO -1.909
Li2O2-LiO8-BeO -5.161
Li2O2-Li2O-BeO -4.807

After sorting, critical chemical potentials and corresponding ranges are [-1.909, -4.807, -5.161, -5.559], as the same as the webapp displayed.

What puzzles me is why chemical potential ranges calculated from the Gibbs phase diagrams are used to construct the grand potential phase diagrams.

Now I realize that two kind of phase diagrams should have the same set of critical chemical potentials. Actually, chemical potentials are hidden info of equilibrium in both the Gibbs phase diagrams and the grand phase diagrams. Different equilibria in the Gibbs phase diagrams may have different chemical potentials, and different the grand phase diagrams contains different chemical potentials. The critical chemical potentials should be able to obtained from both the Gibbs phase diagrams and the grand phase diagrams, but how to calculate critical chemical potentials directly from the grand phase diagrams may still be a problem. Just correct me if anything wrong.


The grand potential phase diagram (PD) is just the regular PD with some of those chemical potentials fixed to whatever the user supplies. If you look at the code in pymatgen, all it does is change the energies to subtract the chemical potential of fixed elements: pymatgen/pymatgen/analysis/ at 676e33da0be7c099af4ecb59d963592d4d56c813 · materialsproject/pymatgen · GitHub

The critical chemical potentials are the transition points, the places where chemical stability changes. We use the middle’s from the regular PD since these should be regions of stable phase equilibria. That way all the stable phases can be identified.

Yes, it is \Phi = G - \mu N.

It is exactly my question, I think it is correct only if the \Phi PD have the same set of equilibrium transitions with the G PD, but how to explain that? Does \Phi = G - \mu N not change transition points of equilibria?

It will if you only pick the transition points along one tie line. Let’s say we make the PD open with respect to Li. The set of transition chemical potentials are not just the critical potentials for Li, but are the Li chemical potential associated with ANY critical point in the PD.

Note how the pymatgen code checks every facet:

The transition chemical potentials are when the PD stability changes. Then EACH segment between transition chemical potentials is a region of constant stability. That’s why we can choose the middle of even 25% alone the way.

I agree with your statement, but what pymatgen did is use segments of chemical potentials derived from facets of Gibbs phase diagram (namely regular PD).

    for facet in self.facets:
        chempots = self._get_facet_chempots(facet)
    def _get_facet_chempots(self, facet):
        complist = [self.qhull_entries[i].composition for i in facet]
        energylist = [self.qhull_entries[i].energy_per_atom for i in facet]
        m = [[c.get_atomic_fraction(e) for e in self.elements] for c in complist]
        chempots = np.linalg.solve(m, energylist)
        return dict(zip(self.elements, chempots))

It is the same with what I wrote above,

In the equilibrium Li-Be-BeO of the Li-Be-O system, we can get \mu_{Li} = -1.909 eV/atom from the equation set,

How to clearly explain the reason of this operation, i.e, regarding transition points of the Gibbs PD (regular PD) as that of the grand phase diagram. In my opinion, these two kind of phase diagrams may have the same set of critical chemical potentials and segmentations.

I’m very confused as to what you’re confused about.
Are you confused that:

  1. The PD and GrandPD have the same transition chemical potentials?
  2. The PD and GrandPD have different transition chemical potentials?
  3. The PD and GrandPD use the same method to get transition chemical potentials?
  4. ???

pymatgen uses 1 to implement codes, I confuse why