Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Create, plot, and analyze OPEN (grand canonical) Phase Diagrams using the pymatgen codebase and Materials Project database
This is a basic example of how to create, plot, and analyze OPEN Phase Diagrams using the pymatgen
codebase and Materials Project database. To run this example, you should:
* have pymatgen ( installed along with matplotlib
* obtain a Materials Project API key (
* paste that API key in the MAPI_KEY variable below, e.g. MAPI_KEY = "foobar1234"
For citation, see
For the accompanying comic book, see
from __future__ import print_function
from pymatgen import MPRester, Element
from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram, \
PhaseDiagram, PDPlotter
def plot_pd(pd, show_unstable=False):
plotter = PDPlotter(pd, show_unstable=show_unstable)
# plotter.write_image("{}.png".format('-'.join(system)), "png") # save figure
def analyze_pd(pd):
print('Stable Entries (formula, materials_id)\n--------')
for e in pd.stable_entries:
print(e.composition.reduced_formula, e.entry_id)
print('\nUnstable Entries (formula, materials_id, e_above_hull (eV/atom), decomposes_to)\n--------')
for e in pd.unstable_entries:
decomp, e_above_hull = pd.get_decomp_and_e_above_hull(e)
pretty_decomp = [("{}:{}".format(k.composition.reduced_formula, k.entry_id), round(v, 2)) for k, v in decomp.items()]
print(e.composition.reduced_formula, e.entry_id, "%.3f" % e_above_hull, pretty_decomp)
if __name__ == "__main__":
MAPI_KEY = None # You must change this to your Materials API key! (or set MAPI_KEY env variable)
system = ["Fe", "P", "O"] # system we want to get open PD for
# system = ["Li", "Fe", "P", "O"] # alternate system example
open_elements_specific = None # e.g., {Element("O"): 0} where 0 is the specific chemical potential
open_element_all = Element("O") # plot a series of open phase diagrams at critical chem pots with this element open
mpr = MPRester(MAPI_KEY) # object for connecting to MP Rest interface
# get data
entries = mpr.get_entries_in_chemsys(system, compatible_only=True)
if open_elements_specific:
gcpd = GrandPotentialPhaseDiagram(entries, open_elements_specific)
plot_pd(gcpd, False)
if open_element_all:
pd = PhaseDiagram(entries)
chempots = pd.get_transition_chempots(open_element_all)
all_gcpds = list()
for idx in range(len(chempots)):
if idx == len(chempots) - 1:
avgchempot = chempots[idx] - 0.1
avgchempot = 0.5 * (chempots[idx] + chempots[idx + 1])
gcpd = GrandPotentialPhaseDiagram(entries, {open_element_all: avgchempot}, pd.elements)
min_chempot = None if idx == len(chempots) - 1 else chempots[idx + 1]
max_chempot = chempots[idx]
print("Chempot range for diagram {} is: {} to {}".format(idx, min_chempot, max_chempot))
plot_pd(gcpd, False)

This comment has been minimized.

Copy link
Owner Author

@computron computron commented Mar 2, 2018

Updated 3/2/18 to support Python 3 and account for latest pymatgen (2018.03.02)


This comment has been minimized.

Copy link

@shyuep shyuep commented Mar 2, 2018

Thanks. Just a suggest that it might be helpful to move this to the matgenb repo as as notebook.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment