One of the more powerful tools in materials screening is the computational phase stability diagram. Unfortunately, it is only utilized at the moment by a few research groups (although I do see its usage increasing), and I thought that a comic book about them might improve the situation.
So here’s that comic book! In addition, this post contains Python examples to create, plot, and analyze phase diagrams using the pymatgen library and Materials Project database. You can now do what earlier took a month of research (computing and generating an entire ternary or quaternary phase diagram) in a few seconds!
This post has three parts:
- The comic!
- Interactive phase diagram examples
- Further resources
The comic!
Click here to download the full high quality PDF version (19MB) of the Phase Diagram comic.
There’s also a small file size version (3MB) for slower connections.
Interactive phase diagram examples
Materials Project Phase Diagram App
The Materials Project phase diagram app allows one to access a database of tens of thousands of DFT calculations and construct interactive computational phase diagrams. You can build binary, ternary, and quaternary diagrams as well as open -element diagrams. No programming required!
Python code example: Creating a phase diagram
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| This is a basic example of how to create, plot, and analyze Phase Diagrams using the pymatgen | |
| codebase and Materials Project database. To run this example, you should: | |
| * have pymatgen (www.pymatgen.org) installed along with matplotlib | |
| * obtain a Materials Project API key (https://www.materialsproject.org/open) | |
| * paste that API key in the MAPI_KEY variable below, e.g. MAPI_KEY = "foobar1234" | |
| For citation, see https://www.materialsproject.org/citing | |
| For the accompanying comic book, see http://www.hackingmaterials.com/pdcomic | |
| """ | |
| from pymatgen import MPRester | |
| from pymatgen.entries.compatibility import MaterialsProjectCompatibility | |
| from pymatgen.phasediagram.pdanalyzer import PDAnalyzer | |
| from pymatgen.phasediagram.pdmaker import PhaseDiagram | |
| from pymatgen.phasediagram.plotter import PDPlotter | |
| if __name__ == "__main__": | |
| MAPI_KEY = None # You must change this to your Materials API key! (or set MAPI_KEY env variable) | |
| system = ["Fe", "P"] # system we want to get PD for | |
| # system = ["Fe", "P", "O"] # alternate system example: ternary | |
| mpr = MPRester(MAPI_KEY) # object for connecting to MP Rest interface | |
| compat = MaterialsProjectCompatibility() # sets energy corrections and +U/pseudopotential choice | |
| # Create phase diagram! | |
| unprocessed_entries = mpr.get_entries_in_chemsys(system) | |
| processed_entries = compat.process_entries(unprocessed_entries) # filter and add energy corrections | |
| pd = PhaseDiagram(processed_entries) | |
| # Plot! | |
| plotter = PDPlotter(pd, show_unstable=False) # you can also try show_unstable=True | |
| plotter.show() | |
| # plotter.write_image("{}.png".format('-'.join(system)), "png") # save figure | |
| # Analyze phase diagram! | |
| pda = PDAnalyzer(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 = pda.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.iteritems()] | |
| print e.composition.reduced_formula, e.entry_id, "%.3f" % e_above_hull, pretty_decomp |
Python code example: Creating a grand canonical phase diagram
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| 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 (www.pymatgen.org) installed along with matplotlib | |
| * obtain a Materials Project API key (https://www.materialsproject.org/open) | |
| * paste that API key in the MAPI_KEY variable below, e.g. MAPI_KEY = "foobar1234" | |
| For citation, see https://www.materialsproject.org/citing | |
| For the accompanying comic book, see http://www.hackingmaterials.com/pdcomic | |
| """ | |
| 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.show() | |
| # 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) | |
| analyze_pd(gcpd) | |
| 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 | |
| else: | |
| 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) | |
| analyze_pd(gcpd) |
Python code example: Checking to see if your materials is stable with respect to compounds in the MP database
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| This is a basic example of how to check a new material for stability using the pymatgen | |
| codebase and Materials Project database. It assumes that you have already calculated | |
| the energy of your predicted compound with DFT. | |
| It is very important to note that this code does NOT do a detailed check to make sure | |
| your computed energies are compatible with those in the Materials Project database | |
| (i.e., uses the same functional, +U corrections, k-point mesh, etc). You must ensure | |
| that the energies you provide are directly comparable with those in Materials Project, | |
| otherwise this code will give you false results! | |
| To run this example, you should: | |
| * have pymatgen (www.pymatgen.org) installed | |
| * obtain a Materials Project API key (https://www.materialsproject.org/open) | |
| * paste that API key in the MAPI_KEY variable below, e.g. MAPI_KEY = "foobar1234" | |
| as well as: | |
| * update MY_COMPOSITION with the composition of your prediction | |
| * update MY_ENERGY_PER_ATOM with the energy per atom (in eV) of your prediction (via DFT) | |
| * update MY_PARAMETERS with the functional and POTCAR code, as well as any hubbards applied | |
| For citation, see https://www.materialsproject.org/citing | |
| For the accompanying comic book, see http://www.hackingmaterials.com/pdcomic | |
| A similar script that directly reads VASP output files and checks stability is at: | |
| https://gist.github.com/shyuep/3570304 | |
| """ | |
| from pymatgen import MPRester, Composition | |
| from pymatgen.entries.compatibility import MaterialsProjectCompatibility | |
| from pymatgen.entries.computed_entries import ComputedEntry | |
| from pymatgen.analysis.phase_diagram import PhaseDiagram | |
| if __name__ == "__main__": | |
| MAPI_KEY = None # You must change this to your Materials API key! (or set MAPI_KEY env variable) | |
| MY_COMPOSITION = Composition("Fe5O8") # put the formula of your prediction | |
| MY_ENERGY_PER_ATOM = –6.5 # note: this must be energy normalized PER ATOM in eV! | |
| MY_PARAMETERS = {"potcar_symbols": ['pbe Fe_pv', 'pbe O'], "hubbards": {"Fe": 5.3}} | |
| mpr = MPRester(MAPI_KEY) # object for connecting to MP Rest interface | |
| compat = MaterialsProjectCompatibility() # sets energy corrections and +U/pseudopotential choice | |
| # Create phase diagram! | |
| unprocessed_entries = mpr.get_entries_in_chemsys([e.symbol for e in MY_COMPOSITION.elements]) | |
| processed_entries = compat.process_entries(unprocessed_entries) # filter and add energy corrections | |
| # define user entry and add corrections | |
| my_entry = ComputedEntry(MY_COMPOSITION, MY_ENERGY_PER_ATOM * MY_COMPOSITION.num_atoms, | |
| parameters=MY_PARAMETERS) | |
| corrections_dict = compat.get_corrections_dict(my_entry) | |
| pretty_corrections = ["{}:{}".format(k, round(v, 3)) for k, v in corrections_dict.items()] | |
| print("We are applying the following corrections (eV) to the user entry: {}".format(pretty_corrections)) | |
| my_entry.correction = sum(corrections_dict.values()) | |
| processed_entries.append(my_entry) | |
| pd = PhaseDiagram(processed_entries) | |
| if my_entry in pd.stable_entries: | |
| print('Your entry {} with energy/atom {} is stable! Its "inverse hull energy" is {} eV/atom.'.format(MY_COMPOSITION, my_entry.energy_per_atom, "%.3f" % pd.get_equilibrium_reaction_energy(my_entry))) | |
| else: | |
| decomp, e_above_hull = pd.get_decomp_and_e_above_hull(my_entry) | |
| pretty_decomp = [("{}:{}".format(k.composition.reduced_formula, k.entry_id), round(v, 2)) for k, v in decomp.items()] | |
| print('Your entry {} with energy/atom {} eV/atom is unstable. Its "decomposition energy" is {} eV/atom. Its predicted decomposition is to: {}'.format(MY_COMPOSITION, "%.3f" % my_entry.energy_per_atom, "%.3f" % e_above_hull, pretty_decomp)) | |
| print('A rule of thumb is decomposition energy less than about 0.100 eV/atom might be considered "metastability", although this is not a quantitative statement.') | |
| print('For more details, see Sun, W., Dacek, S. T., Ong, S. P., Hautier, G., Jain, A., Richards, W. D., Gamst, A. C., Persson, K. A. & Ceder, G. The thermodynamic scale of inorganic crystalline metastability. (2016).') | |
Further resources
“Accuracy of ab initio methods in predicting the crystal structures of metals: A review of 80 binary alloys” by Curtarolo et al.
This (somewhat epic!) paper contains data for 80 binary convex hulls computed with density functional theory. The results are compared with known experimental data and it is determined that the degree of agreement between computational and experimental methods is between 90-97%.
“A Computational Investigation of Li9M3(P2O7)3(PO4)2 (M = V, Mo) as Cathodes for Li Ion Batteries” by Jain et al.
The endpoints of a binary convex hull need not be elements. For example, in the Li ion battery field one often searches for stable intermediate phases that form at certain compositions as lithium is inserted into a framework structure. The paper above is just one example of many computational Li ion battery papers that use such “pseudo-binary” convex hulls.
“Configurational Electronic Entropy and the Phase Diagram of Mixed-Valence Oxides: The Case of LixFePO4” by Zhou et al.
Incorporating temperature into first-principles convex hulls is often possible, but not always straightforward or easy to do. Here is one example of this process that focuses on electronic entropy.
Wikipedia article on Ternary plots
Ternary plots are not only for phase diagrams (the most creative usage I’ve ever seen is in Scott McCloud’s Understanding Comics, where it is used to explain the language of art and comics). Wikipedia does a good job of explaining the basics of how to read and interpret compositions on ternary diagrams.
“Li-Fe-P-O2 phase diagram from first principles calculations” by Ong et al.
Here is a nice example of the computation of a quaternary phase diagram – sliced into ternary sections – from first principles calculations.
“Accuracy of density functional theory in predicting formation energies of ternary oxides from binary oxides and its implication on phase stability” by Hautier et al.
How accurate are computational phase diagrams? The correct answer, like always, is “it’s complicated”. But based on results from this paper and some experience, colleagues of mine and I have found that an error bar of 25 meV/atom is usually a good estimate. We usually double that to 50 meV/atom when searching for materials to synthesize by conventional methods.
“Formation enthalpies by mixing GGA and GGA + U calculations” by Jain et al.
In an ideal world, first principles calculations would live up to their name and require no adjustable parameters. In practice, however, DFT errors do not always cancel when comparing energies of compounds with different types of electronic states. This paper shows how one can mix two DFT approximations along with some experimental data in order to produce a correct phase diagram across a changing landscape of electronic states.
“First-Principles Determination of Multicomponent Hydride Phase Diagrams: Application to the Li-Mg-N-H System” by Akbarzadeh et al.
An alternate (but equivalent) approach to the convex hull algorithm for determining phase diagrams is to use a linear programming approach. This is demonstrated by Akbarzadeh et al. in the search for H2 sorbents.
“Thermal stabilities of delithiated olivine MPO4 (M = Fe, Mn) cathodes investigated using first principles calculations” by Ong et al.
If Li ion battery cathode materials (generally oxygen-containing compounds) release O2 gas from their lattice, it can lead to runaway electrolyte reactions that cause fire. Thus, a safe cathode material resists O2 release even under extreme conditions. Stated another way, safety is the “price point” (inverse O2 chemical potential) at which a cathode material will give up its oxygen. The higher the price point, the more stable the compound. This paper compares the critical chemical potential for O2 release between MnPO4 and FePO4 cathode materials, finding that similar chemistry and structure doesn’t necessarily imply similar safety.
“CO2 capture properties of M–C–O–H (M.Li, Na, K) systems: A combined density functional theory and lattice phonon dynamics study” by Duan et al.
The CO2 capture problem is to find a compound that absorbs CO2 from an open environment at chemical potentials found in industrial processes, and then releases the CO2 back into some other open environment under sequestration conditions. This paper constructs multi-dimensional phase diagrams to predict how different chemical systems will react with CO2 under different conditions.
HACKING MATERIALS 
I just stumbled upon this blog during my research and I think it is really awesome! I’m very impressed.
I had a quick question about phase diagrams (I don’t know if this is the right place to ask this). When we plot the phase diagram on using the materials project api, this returns the phase diagram at 0 K and 0 atm, am I correct? This is different from the “Li−Fe−P−O2 Phase Diagram from First Principles Calculations” paper where they use 0.1 MPa as the reference pressure?
So what I would obtain on the materials project would be somewhat shifted in chemical potential when compared to the paper?
Hi Prateek
Glad you found it and thanks!
For the phase diagram question, it’s probably best to contact Materials Project support with such inquiries. But the Materials Project does the same thing as in the Li-Fe-P-O2 paper – in both instances, all solid phases are calculated without any pressure (0 atm), but the gas phase O2 energy is fit (not necessarily calculated) using experimental data for standard conditions (1 atm, 293K). So even in the paper, it is not as if the solid phases are being adjusted for pressure.
For the solid phases, having 1 atm of pressure versus zero pressure really doesn’t change anything at all – meaning if we (or the paper) did re-calculate all the solid phases at 1 atm, we wouldn’t expect anything different. So solid phase calculations performed at 0 atm (and to a certain extent, 0K) can for practical purposes be assumed to be close to “standard temperature and pressure”.
Dear ANUBHAV, I’ve analyzed your “check_stability.py” code and interest how it can be modified:
1. what about multi component composition consisting 5-6 and even more elements, is it able?
2. is it possible for your code to simply provide the compounds of ground states generated by the convex-hull construction, so that user is able to compute them by himself by his own DFT software and thus analyze the stability by his own DFT software. This can be important because some users can not have VASP or not prefer to use it
Hi Maxim,
1. Yes, it should be able to do any number of elements without any modifications.
2. Yes, just skip the code about adding your own calculation entry. You can still call pd.stable_entries to get all the stable entries.
There are many other functions available in the code if you look around the PhaseDiagram and PDAnalyzer classes in pymatgen. If you do give it a try and get stuck, you can always post a question to:
https://groups.google.com/forum/#!forum/pymatgen
However, you should give it a try first and send a report/question if you’re unable to get it working after giving it your best shot.
Thanks a lot
I modified it a little for my python and pymatgen and it works fine!
The main issues were related to inconsistent names in import section etc.
Ok! Glad you got it to work. Yes, the underlying libraries like pymatgen have changed since I first posted the code example. If you want, you could submit a pull request to Github with the updated import statements for future users.
Let me send it to your mail so you post it more conveniently
Thank you again)
my email is: ars21031960@gmail.com
Hi Maxim – thanks! I got your email and updated the code snippet based on the parts that changed (import statements and Py3 compatible print statements)
Dear Anubhav,
Thank you a lot; using your code I was able to publish in a good journal
http://pubs.acs.org/doi/abs/10.1021/acs.jpcc.7b01575
My question is to the code “check_stability.py” and phrase in it: “decomposition energy less than about 0.100 eV/atom might be considered “metastability””
Could you give a reference for this phrase, if to write it to the “print” code it will be useful. The users will be able to read the corresponding paper or will be able to give a reference in their own papers.
The second question is to the same 0.100 eV/atom value. For the Li2MnSiO4 compound there is difficulty in obtaining single-phase product (several phases coexist) and the difference in energy between them is within 0.010 eV/atom.
For ZrO2 only single phase exist in this 0.010 eV/atom range and indeed it is single phase monoclinic form.
I updated the gist, see Sun, W., Dacek, S. T., Ong, S. P., Hautier, G., Jain, A., Richards, W. D., Gamst, A. C., Persson, K. A. & Ceder, G. The thermodynamic scale of inorganic crystalline metastability. (2016).
Thank you again! Good luck to your project!
Hi Anubhav,
The code “check_stability.py” is outdated and gives many errors, I believe this is due upgrades in pymatgen
Hi Anubhav,
I’ve fixed it:
replace “from pymatgen.phasediagram.analyzer import PDAnalyzer
from pymatgen.phasediagram.maker import PhaseDiagram”
with
“from pymatgen.analysis.phase_diagram import PhaseDiagram”,
delete ” pda = PDAnalyzer(pd)”
replace “pda” with “pd”
Ok I updated the code – discussion migrated to the gist itself (rather than the blog)
Super!