This notebook shows how to
Written using:
We use the Materials Project API to obtain energies of compounds.
In [1]:
from pymatgen.ext.matproj import MPRester
from pymatgen.analysis.interface_reactions import InterfacialReactivity
from pymatgen.analysis.phase_diagram import PhaseDiagram, GrandPotentialPhaseDiagram
from pymatgen import Composition, Element
%matplotlib inline
# Initialize the REST API interface. You may need to put your own API key in as an arg.
mpr = MPRester()
First, set the values of the two reactants. Optionally, simulate the case where the reaction system is in contact with an elemental reservoir.
Because the methodology here is to generate a pseudo-binary phase stability diagram of two reactants as a function of mixing ratio, the addition of an elemental reservoir implies construction of a so-called grand potential phase diagram.
In [2]:
# Chemical formulae for two solid reactants.
reactant1 = 'LiCoO2'
reactant2 = 'Li3PS4'
# Is the system open to an elemental reservoir?
grand = True
if grand:
# Element in the elemental reservoir.
open_el = 'Co'
# Relative chemical potential vs. pure substance. Must be non-positive.
relative_mu = -1
Now, compile the critical reaction information:
In [3]:
# Get the compositions of the reactants
comp1 = Composition(reactant1)
comp2 = Composition(reactant2)
# Gather all elements involved in the chemical system.
elements = [e.symbol for e in comp1.elements + comp2.elements]
if grand:
elements.append(open_el)
elements = list(set(elements)) # Remove duplicates
# Get all entries in the chemical system
entries = mpr.get_entries_in_chemsys(elements)
# Build a phase diagram using these entries.
pd = PhaseDiagram(entries)
# For an open system, include the grand potential phase diagram.
if grand:
# Get the chemical potential of the pure subtance.
mu = pd.get_transition_chempots(Element(open_el))[0]
# Set the chemical potential in the elemental reservoir.
chempots = {open_el: relative_mu + mu}
# Build the grand potential phase diagram
gpd = GrandPotentialPhaseDiagram(entries, chempots)
# Create InterfacialReactivity object.
interface = InterfacialReactivity(
comp1, comp2, gpd, norm=True, include_no_mixing_energy=True, pd_non_grand=pd, use_hull_energy=False)
else:
interface = InterfacialReactivity(
comp1, comp2, pd, norm=True, include_no_mixing_energy=False, pd_non_grand=None, use_hull_energy=False)
From here, you can plot reaction energy versus mixing ratio as below. Note that the mixing ratio is between the normalized compositions of the reactants.
In [4]:
plt = interface.plot()
In [5]:
from collections import OrderedDict
import pandas as pd
critical_rxns = [
OrderedDict([
("Atomic fraction", round(ratio, 3)),
("Reaction equation", rxn),
("E$_{rxt}$ per mol equation (kJ/mol)", round(rxn_energy, 1)),
("E$_{rxt}$ per reactant atom (eV/atom)", round(reactivity, 3)),
])
for _, ratio, reactivity, rxn, rxn_energy in interface.get_kinks()]
interface_reaction_table = pd.DataFrame(critical_rxns)
interface_reaction_table
Out[5]:
You can also obtain the mixing ratio between the original compositions, i.e. mol fraction of the first reactant.
In [6]:
import numpy as np
interface_reaction_table['Molar fraction'] = pd.Series(np.round(interface.get_critical_original_kink_ratio(), 3))
interface_reaction_table
Out[6]:
Note that the reaction equations are Reaction objects suitable for structured analysis:
In [7]:
rxn = critical_rxns[2]["Reaction equation"]
print(rxn)
print(type(rxn))
In [8]:
# Get interface reaction information for reactants LiCoO2 and Li3PS4 in open system to Co.
kinks_from_API = mpr.get_interface_reactions('LiCoO2','Li3PS4', open_el='Co', relative_mu=-1, use_hull_energy=False)
# Get inforamtion for the second critical reaction.
print(kinks_from_API[1])
The critical reaction information from REST API should be the same as in the previous table:
In [9]:
critical_rxns_from_API = [
OrderedDict([
("Atomic fraction", round(reaction['ratio_atomic'], 3)),
("Molar fraction", round(reaction['ratio_molar'], 3)),
("Reaction equation", reaction['rxn_str']),
("E$_{rxt}$ per mol equation (kJ/mol)", round(float(reaction['rxn_energy_sigdig_kJmol']),1)),
("E$_{rxt}$ per reactant atom (eV/atom)", round(reaction['energy'], 3))
])
for reaction in kinks_from_API]
pd.DataFrame(critical_rxns_from_API)
Out[9]:
In [ ]: