Plotting a pourbaix diagram using the Materials Project API

If you use this infrastructure, please consider citing the following work:

K. A. Persson, B. Waldwick, P. Lazic and G. Ceder, Phys. Rev. B, 2012, 85, 235438.

A. K. Singh, L. Zhou, A. Shinde, S. K. Suram, J. H. Montoya, D. Winston, J. M. Gregoire and K. A. Persson, Chem. Mater., 2017, DOI: acs.chemmater.7b03980.

Notebook Author: Joseph Montoya

---

The Materials Project REST interface includes functionality to construct pourbaix diagrams from computed entries. Note that the Pourbaix infrastructure is still undergoing revision, but now includes a simplified interface that enables MP and pymatgen users to fetch entries that have been processed according to the Materials Project Aqueous Compatibility scheme. Thus, users can reproduce web Pourbaix diagrams in two or three steps in pymatgen.

# Uncomment the subsequent lines in this cell to install dependencies for Google Colab.
# !pip install pymatgen==2022.7.19

# Import necessary tools from pymatgen
from pymatgen.analysis.pourbaix_diagram import PourbaixDiagram, PourbaixPlotter
from pymatgen.ext.matproj import MPRester

%matplotlib inline

# Initialize the MP Rester
mpr = MPRester()

To retrieve entries necessary to construct a Pourbaix Diagram use MPRester.get_pourbaix_entries(LIST_OF_ENTRIES) with a list of entries comprising your chemical system. It is not necessary to include 'O' and 'H' in your list, as they are added automatically. This function also contains all of necessary preprocessing to ensure the PourbaixEntries are compatible with the pymatgen PourbaixDiagram constructor.

# Get all pourbaix entries corresponding to the Cu-O-H chemical system.
entries = mpr.get_pourbaix_entries(["Cu"])

Pourbaix diagrams can be constructed using PourbaixDiagram(RETRIEVED_ENTRIES) as below. Note that a comp_dict keyword argument may also be supplied to the PourbaixDiagram constructor if a fixed composition for a multi-element pourbaix diagram is desired.

# Construct the PourbaixDiagram object
pbx = PourbaixDiagram(entries)

The PourbaixAnalyzer includes a number of useful functions for determining stable species and stability of entries relative to a given pourbaix facet (i.e. as a function of pH and V).

# Get an entry stability as a function of pH and V
entry = [e for e in entries if e.entry_id == "mp-1692"][0]
print(
    "CuO's potential energy per atom relative to the most",
    "stable decomposition product is {:0.2f} eV/atom".format(
        pbx.get_decomposition_energy(entry, pH=7, V=-0.2)
    ),
)

CuO's potential energy per atom relative to the most stable decomposition product is 0.44 eV/atom

This suggests that CuO, for example, has a large driving force for decomposition at neutral pH and mildly reducing conditions.

To see this in more detail, we can plot the pourbaix diagram. The PourbaixPlotter object is also initialized using an instance of the PourbaixDiagram object.

plotter = PourbaixPlotter(pbx)

plotter.get_pourbaix_plot().show()

The PourbaixPlotter object can also plot the relative stability of an entry across the pourbaix diagram. To do this, use the PourbaixPlotter.plot_entry_stability method.

plt = plotter.plot_entry_stability(entry)
plt.show()

Plotting k-nary systems

Pymatgen also supports binary/ternary pourbaix diagrams with fixed compositions of the non-H or O elements. This is achieved by finding all possible combinations of entries that fulfill the composition constraint and treating them as individual entries in pourbaix space. Note that you can supply a composition dictionary and to further tune the pourbaix diagram.

# Get all pourbaix entries corresponding to the Fe-O-H chemical system.
entries = mpr.get_pourbaix_entries(["Bi", "V"])
# Construct the PourbaixDiagram object
pbx = PourbaixDiagram(
    entries,
    comp_dict={"Bi": 0.5, "V": 0.5},
    conc_dict={"Bi": 1e-8, "V": 1e-8},
    filter_solids=True,
)

Note that the filter_solids keyword argument in the PourbaixDiagram instantiation above tells the constructor whether to filter solids by phase stability on the compositional phase diagram. Note that Pourbaix Diagrams generated with and without this argument may look significantly different in the OER and HER regions, since highly oxidized materials (e. g. Bi$_2$O$_5$) or highly reduced materials (e. g. most hydrides) may not be stable on the compositional phase diagram. The filtering process significantly reduces the time it takes to generate all of the combined entries for a binary or ternary pourbaix diagram though, so it may be prudent to use in those cases.

# Construct the pourbaix analyzer
plotter = PourbaixPlotter(pbx)
plt = plotter.get_pourbaix_plot()
plt.show()

Getting the heatmaps for a solid entry in these cases is a bit more involved, because many of the regions of the pourbaix diagram include multiphase entries. Here's an example for this case.

bivo4_entry = [entry for entry in entries if entry.entry_id == "mp-613172"][0]
plt = plotter.plot_entry_stability(bivo4_entry)

If a compound is provided that doesn't meet the constraints of the composition (in this case, equal parts Bi and V), the plotter will attempt to find the most stable entry containing that solid and any spectating ions. In this case, it estimates the stability of BiO$_2$(s) + VO$_4$-.

bio2_entry = [entry for entry in entries if entry.entry_id == "mp-557993"][0]
plt = plotter.plot_entry_stability(bio2_entry)