Skip to content

Instantly share code, notes, and snippets.

@computron
Last active August 18, 2023 16:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save computron/0c0ed87f3c239409af13 to your computer and use it in GitHub Desktop.
Save computron/0c0ed87f3c239409af13 to your computer and use it in GitHub Desktop.
Check a new material's stability using the pymatgen codebase and Materials Project database
"""
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).')
@mxm2
Copy link

mxm2 commented Aug 25, 2018

You right, 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"

@computron
Copy link
Author

Thanks, I updated the code.

But now there is a simpler method - MPRester directly has a method called get_stability()

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