This notebook demonstrates how you can calculate reaction energies using the Materials API and pymatgen.

Written using:

  • pymatgen==2018.3.13
In [1]:
from pymatgen.ext.matproj import MPRester
from pymatgen import Composition
from pymatgen.entries.computed_entries import ComputedEntry
from pymatgen.core.units import FloatWithUnit
from pymatgen.analysis.reaction_calculator import ComputedReaction

#This initializes the REST adaptor. Put your own API key in if necessary.
a = MPRester()

#This gets all entries belonging to the Ca-C-O system.
all_entries = a.get_entries_in_chemsys(['Ca', 'C', 'O'])

#This method simply gets the lowest energy entry for all entry with the same composition.
def get_most_stable_entry(formula):
    relevant_entries = [entry for entry in all_entries if entry.composition.reduced_formula == Composition(formula).reduced_formula]
    relevant_entries = sorted(relevant_entries, key=lambda e: e.energy_per_atom)
    return relevant_entries[0]

CaO = get_most_stable_entry("CaO")
CO2 = get_most_stable_entry("CO2")
CaCO3 = get_most_stable_entry("CaCO3")

reaction = ComputedReaction([CaO, CO2], [CaCO3])
energy = FloatWithUnit(reaction.calculated_reaction_energy, "eV atom^-1")

print("Reaction energy = {}".format("kJ mol^-1")))

# The following portions demonstrate how to get the experimental values as well.
exp_CaO = a.get_exp_entry("CaO")
exp_CaCO3 = a.get_exp_entry("CaCO3")

#Unfortunately, the Materials Project database does not have gas phase experimental entries. This is the value from NIST. We manually create the entry.
#Exp entries should be in kJ/mol.
exp_CO2 = ComputedEntry("CO2", -393.51)

exp_reaction = ComputedReaction([exp_CaO, exp_CO2], [exp_CaCO3])

print("Reaction energy = {} kJ mol^-1".format(exp_reaction.calculated_reaction_energy))
CaO + CO2 -> CaCO3
Reaction energy = -145.39165028828842 kJ mol^-1
CaO + CO2 -> CaCO3
Reaction energy = -178.30000000000064 kJ mol^-1