1) Enumerators#

Author: Matthew McDermott (U.C. Berkeley / Lawerence Berkeley National Laboratory)

Last Updated: 08/31/23

In this notebook, we cover reaction enumeration in the reaction-network package. The reaction enumerators are the core of the package and supply all reaction data to be used in reaction network construction and downstream analyses. This is the recommended starting point for working with the codebase. Please see the next notebook(s) for further applications (such as network construction).

If you use this package in your work, please consider citing the following paper:

McDermott, M. J., Dwaraknath, S. S., and Persson, K. A. (2021). A graph-based network for predicting chemical reaction pathways in solid-state materials synthesis. Nature Communications, 12(1). https://doi.org/10.1038/s41467-021-23339-x

1. Imports#

First, make sure the reaction-network package and Materials Project API package are installed. Run these commands:

pip install --upgrade reaction-network

pip install --upgrade mp-api

in your terminal. You may need to restart/reconfigure your Jupyter kernel after you do this.

[1]:

import logging
from pprint import pprint

from mp_api.client import MPRester
from pymatgen.core.periodic_table import Element

from rxn_network.core import Composition
from rxn_network.costs.functions import Softplus
from rxn_network.enumerators.basic import BasicEnumerator, BasicOpenEnumerator
from rxn_network.enumerators.minimize import MinimizeGibbsEnumerator, MinimizeGrandPotentialEnumerator
from rxn_network.entries.entry_set import GibbsEntrySet
from rxn_network.reactions.reaction_set import ReactionSet

# this is useful if you are editing your code locally!
%load_ext autoreload
%autoreload 2

2. Downloading and modifying entries#

We will work with an example chemical system: yttrium (Y), manganese (Mn), and oxygen (o).

First, we need to acquire thermodynamic data for phases in this system from the Materials Project (MP), a computed materials database containing calculations for 150,000+ materials.

[2]:

with MPRester() as mpr:  # insert your Materials Project API key if it's not stored as an environment variable
    entries = mpr.get_entries_in_chemsys("Y-Mn-O")

A unique feature of the reaction-network package is the GibbsEntrySet class.

This class allows us to automatically convert ComputedStructureEntry objects downloaded from the MP database into GibbsComputedEntry objects, where DFT-calculated energies have been converted to an AI-estimated equivalent values of the Gibbs free energies of formation, \(\Delta G_f\) for all entries at the specified temperature.

This will also include some experimental thermochemistry data (e.g., NIST-JANAF). For more information, check out the citation in the documentation for GibbsComputedEntry.

[3]:

temp = 900  # units: Kelvin
gibbs_entries = GibbsEntrySet.from_computed_entries(entries, temperature=temp, include_nist_data=True)

We can print the entries by calling .entries or .entries_list:

[4]:

gibbs_entries.entries_list[:5]  # the first five entries in the Y-Mn-O system (in alphabetical order)

[4]:

[GibbsComputedEntry | mp-723285-GGA | O8 (O2)
 Gibbs Energy (900 K) = 0.0000,
 GibbsComputedEntry | mp-1238773-GGA+U | Mn1 O1 (MnO)
 Gibbs Energy (900 K) = -2.6169,
 GibbsComputedEntry | mp-1238899-GGA+U | Mn1 O1 (MnO)
 Gibbs Energy (900 K) = -2.9922,
 GibbsComputedEntry | mp-25223-GGA+U | Mn1 O2 (MnO2)
 Gibbs Energy (900 K) = -3.7579,
 GibbsComputedEntry | mp-796077-GGA+U | Mn1 O2 (MnO2)
 Gibbs Energy (900 K) = 0.6017]

The GibbsEntrySet class has many helpful functions, such as the following filter_by_stability() function, which automatically removes entries which are a specified energy per atom above the convex hull of stability:

[5]:

filtered_entries = gibbs_entries.filter_by_stability(0.025)

You should now see a much shorter list of entries within the Y-Mn-O system (< 25 meV/atom below the hull):

[6]:

print(f"{len(gibbs_entries)} (unfiltered) -> {len(filtered_entries)} (filtered)")

108 (unfiltered) -> 13 (filtered)

Another useful function is the get_min_entry_by_formula function. This automatically finds the entry with the lowest energy matching the provided formula (composition). This is useful for querying the ground-state polymorph when there are many entries for a particular composition.

[7]:

gibbs_entries.get_min_entry_by_formula("Mn2O3")

[7]:

GibbsComputedEntry | mp-1172875-GGA+U | Mn32 O48 (Mn2O3)
Gibbs Energy (900 K) = -121.1202

3. Running enumerators#

Now that we’ve discussed creating a set of entries, we will learn how to enumerate reactions from those entries.

There are four distinct enumerator classes contained within rxn_network.enumerators: These are:

(basic) 1. BasicEnumerator: uses a combinatorial approach to identify all possible (closed) reactions within a set of entries. 2. BasicOpenEnumerator: uses a combinatorial approach to identify all open reactions within a set of entries and a list of specified open entries/elements.

(minimize)

  1. MinimizeGibbsEnumerator: uses a thermodynamic approach to identify all reactions within a set of entries that are predicted by minimizing Gibbs free energy between a set of two reacting phases touching at an interface.

  2. MinimizeGrandPotentialEnumerator: uses a thermodynamic approach to identify all reactions within a set of entries that are predicted by minimizing the grand potential energy between a set of two reacting phases touching at an interface with an open element at a specified chemical potential.

There is no “correct” enumerator to use; while the thermodynamic (i.e., minimize) enumerators may seem more logical, our thermodynamic data is not always exaclty correct – this means that some reactions will not be enumerated.

a) Basic enumerators#

We first create a basic enumerator object by initializing one from the BasicEnumerator class. Let’s initialize it with only default arguments:

[8]:

be = BasicEnumerator()

The BasicEnumerator class, as is true for all other enumerator classes, has many arguments that will customize its output. To view helpful documentation for these arguments, press Shift+Tab with your cursor inside the parentheses in the cell above.

A key argument to the basic enumerators is the maximum reactant/product cardinality, \(n\). The default is \(n=2\). This setting suffices for capturing many of the reactions that make up a solid-state reaction pathway according to the pairwise interface reaction hypothesis (most powder reactions proceed through reactions of interfacial pairs). Note that it is possible to set \(n=3\), however this dramatically increases the combinatorial complexity and can take a long time for more complex systems.

In general, the default arguments are good for generating a list of simple (unconstrained) reactions, as we might build for a reaction network. Run the following cell:

[9]:

all_rxns = be.enumerate(filtered_entries)

2023-08-31 11:22:12,439 INFO rxn_network.utils.ray Ray is not initialized. Checking for existing cluster...
2023-08-31 11:22:12,439 INFO rxn_network.utils.ray Could not identify existing Ray instance. Creating a new one...

2023-08-31 11:22:14,123 INFO worker.py:1627 -- Started a local Ray instance. View the dashboard at 127.0.0.1:8266 

2023-08-31 11:22:14,616 INFO rxn_network.utils.ray HOST: Matts-MBP.local, Num CPUs: 12.0, Total Memory: 37544088372.0

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1570.46it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:01<00:00,  1.54s/it]

You may notice that something called Ray is initialized. This is a python library for parallelizing functions and is used by the reaction-network code to parallelize enumeration.

The cell above should have completed somewhat quickly – ideally within a second or two. As a result, a list of 692 generated reactions will be stored within the all_rxns object. Note: this number may change in the future if the MP database changes…

[10]:

print(len(all_rxns))

692

Every time enumerate() is called, a ReactionSet will be returned. This is a memory-efficient object that can be used to store large sets of reactions.

The ReactionSet class stores sets of reactions as arrays. Note that the actual reaction objects can only be accessed by iterating through the reaction set. Lets print the “first” 10 reactions. These may be different on your side; the reactions are generated in no particualr order.

[11]:

for count, r in enumerate(all_rxns):
    print(r)
    if count>10:
        break

2 Y2Mn2O7 -> O2 + 4 YMnO3
O2 + 4 YMnO3 -> 2 Y2Mn2O7
Y2Mn2O7 -> Y2O3 + 2 MnO2
Y2O3 + 2 MnO2 -> Y2Mn2O7
2 YMnO3 -> Mn2O3 + Y2O3
Mn2O3 + Y2O3 -> 2 YMnO3
2 YMn2O5 -> Y2Mn2O7 + Mn2O3
Y2Mn2O7 + Mn2O3 -> 2 YMn2O5
YMn2O5 -> YMnO3 + MnO2
YMnO3 + MnO2 -> YMn2O5
6 Mn2O3 -> O2 + 4 Mn3O4
O2 + 4 Mn3O4 -> 6 Mn2O3

Looking at the list of reactions above, we see that all reactions are stoichiometrically balanced. If we look at any particular reaction object, we find that the reaction energy and uncertainty are properties which can be easily accessed:

[12]:

r = list(all_rxns)[0]
print(r)
print(f"{round(r.energy_per_atom, 4)} ± {round(r.energy_uncertainty_per_atom, 2)} eV/atom")

2 Y2Mn2O7 -> O2 + 4 YMnO3
0.038 ± 0.07 eV/atom

What if you want to enumerate all reactions from a known set of precursors? An inefficient way to do so might be to filter the previous reaction set by your phases of interest…

However, a more efficient solution has been provided. We can supply our precursor formulas when we initialize the BasicEnumerator object. This will reduce the number of calculations required significantly. Lets say we have Y2O3 as a precursor:

[13]:

be_precursors = BasicEnumerator(precursors=["Y2O3"])
y2o3_rxns_exclusive = be_precursors.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 1/1 [00:00<00:00, 1792.44it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00, 521.42it/s

[14]:

for r in y2o3_rxns_exclusive:
    print(r)

0.6667 Y2O3 -> O2 + 1.333 Y

Note that by default, this only produces reactions which exclusively have the provided precursor(s).

To include reactions that contain this precursor (and possibly others) set the exclusive_precursors=False:

[15]:

be_precursors = BasicEnumerator(precursors=["Y2O3"], exclusive_precursors=False)
y2o3_rxns = be_precursors.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 2344.82it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00,  4.25it/s]

We now get a much larger list of reactions, all of which contain Y2O3 as a precursor (and often another phase such as a manganese oxide):

[16]:

pprint(list(y2o3_rxns)[:10])  # a sample of just 10 reactions

[Y2O3 + 2 MnO2 -> Y2Mn2O7,
 Mn2O3 + Y2O3 -> 2 YMnO3,
 0.6667 Y2O3 -> O2 + 1.333 Y,
 Mn5O8 + 0.5 Y2O3 -> 0.5 Y2Mn2O7 + 2 Mn2O3,
 Mn5O8 + 2.5 Y2O3 -> 0.5 Y2Mn2O7 + 4 YMnO3,
 Y2O3 + 2 YMn2O5 -> Y2Mn2O7 + 2 YMnO3,
 Mn2O3 + 0.25 Y2O3 -> 0.25 Y2Mn2O7 + 0.5 Mn3O4,
 Mn5O8 + Y2O3 -> Y2Mn2O7 + Mn3O4,
 Y2O3 + 8 YMn2O5 -> 5 Y2Mn2O7 + 2 Mn3O4,
 Mn2O3 + 0.75 Y2O3 -> 0.75 Y2Mn2O7 + 0.5 Mn]

This same approach can be used for the target phase(s) as well.

With exclusive_targets=True:

[17]:

be_target = BasicEnumerator(targets=["YMnO3"], exclusive_targets=True)
ymno3_rxns_exclusive = be_target.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 2158.95it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00, 28.38it/s]

[18]:

for r in ymno3_rxns_exclusive:
    print(r)

Mn2O3 + Y2O3 -> 2 YMnO3

Due to the fact that the right side of the reaction can only contain YMnO3, we get just one reaction above!

Now, with exclusive_targets=False (the default):

[19]:

be_target = BasicEnumerator(targets=["YMnO3"])
ymno3_rxns = be_target.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 2322.75it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00,  4.16it/s]

[20]:

pprint(list(ymno3_rxns)[:10])

[2 Y2Mn2O7 -> O2 + 4 YMnO3,
 Mn2O3 + Y2O3 -> 2 YMnO3,
 YMn2O5 -> YMnO3 + MnO2,
 0.5 Y2Mn2O7 + 0.5 Mn2O3 -> YMnO3 + MnO2,
 0.5 Y2Mn2O7 + 2.5 Mn2O3 -> Mn5O8 + YMnO3,
 Mn5O8 + 2.5 Y2O3 -> 0.5 Y2Mn2O7 + 4 YMnO3,
 Y2O3 + 2 YMn2O5 -> Y2Mn2O7 + 2 YMnO3,
 0.3333 Y2Mn2O7 + 0.6667 Mn3O4 -> Mn2O3 + 0.6667 YMnO3,
 0.5 Y2Mn2O7 + 0.25 Mn3O4 -> YMnO3 + 0.75 MnO2,
 1.333 Y2Mn2O7 + 1.667 Mn3O4 -> Mn5O8 + 2.667 YMnO3]

And finally, with multiple targets specified (e.g., YMnO3 and O2):

[21]:

be_targets = BasicEnumerator(targets=["YMnO3", "O2"], exclusive_targets=True)
ymno3_rxns_o2 = be_targets.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 2457.12it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00, 13.44it/s]

[22]:

for r in ymno3_rxns_o2:
    print(r)

2 Y2Mn2O7 -> O2 + 4 YMnO3
Mn2O3 + Y2O3 -> 2 YMnO3
Y2O3 + 2 MnO2 -> 0.5 O2 + 2 YMnO3
Mn5O8 + 2.5 Y2O3 -> 0.25 O2 + 5 YMnO3
Y2O3 + 2 YMn2O5 -> 0.5 O2 + 4 YMnO3

What about open gases, liquids, etc.??#

In the previous cell, we showed that it was possible to specify YMnO\(_3\) as a target, along with O\(_2\). However, because O\(_2\) is a gas, it is often desirable to include it as an open entry in addition to the 2 possible precursors or targets. For example, we may want to specify a reaction as follows:

\[A + B ~ +~\textrm{O}_2 \rightarrow C + D\]
\[\textrm{or}\]
\[A + B \rightarrow C + D ~ +~\textrm{O}_2\]

To do this, we use the BasicOpenEnumerator class, which is an extension to the original basic enumerator. All of the lessons learned above also apply to this class, although now a list of open entry formulas (e.g., ["O2"]) must be specified:

[23]:

be_target_open = BasicOpenEnumerator(open_phases=["O2"], targets=["YMnO3"])
ymno3_rxns_open = be_target_open.enumerate(filtered_entries)

Building chunks...: 100%|████████████████████| 3/3 [00:00<00:00, 1811.79it/s]
Enumerating reactions (BasicOpenEnumerator): 100%|█| 2/2 [00:01<00:00,  1.82i

The BasicOpenEnumerator generally requires more time to run as it considers a larger combinatorial space.

You will now see that the reactions from before are included in this new list, along with many other open-O\(_2\) reactions – some which even have a total of 3 reactants or 3 products. In other words, the open entry/entries do not count towards the specified maximum cardinality, \(n\).

[24]:

for r in list(ymno3_rxns_open)[:10]:  # first 10 reactions
    print(r)

2 Y2Mn2O7 -> O2 + 4 YMnO3
O2 + 0.6667 Mn2O3 + 1.333 Y -> 1.333 YMnO3
O2 + 4 Mn3O4 + 6 Y2O3 -> 12 YMnO3
O2 + 0.4 Mn3O4 + 1.2 Y -> 1.2 YMnO3
O2 + 1.333 Mn + 0.6667 Y2O3 -> 1.333 YMnO3
O2 + 0.6667 Mn + 0.6667 Y -> 0.6667 YMnO3
O2 + 4 MnO + 2 Y2O3 -> 4 YMnO3
O2 + MnO + Y -> YMnO3
O2 + 2 MnO2 + 2 Y -> 2 YMnO3
O2 + 0.2857 Mn5O8 + 1.429 Y -> 1.429 YMnO3

It is even possible to specify multiple open entries, allowing for more complex reactions. For example, specifiyng Y2O3 and O2 as both being open:

[25]:

be_target_open2 = BasicOpenEnumerator(open_phases=["Y2O3", "O2"],targets=["YMnO3"])
ymno3_rxns_open2 = be_target_open2.enumerate(filtered_entries)

Building chunks...: 100%|██████████████████████| 2/2 [00:00<00:00, 28.43it/s]
Enumerating reactions (BasicOpenEnumerator): 100%|█| 6/6 [00:01<00:00,  3.15i

[26]:

for r in list(ymno3_rxns_open2)[:10]:  # first 10 reactions
    print(r)

0.5 Y2Mn2O7 + 0.3333 Y -> YMnO3 + 0.1667 Y2O3
O2 + 0.1026 YMn12 -> 0.5641 Mn2O3 + 0.1026 YMnO3
Y2O3 + 0.6667 Mn3O4 + 0.1667 O2 -> 2 YMnO3
Mn + 0.5 Y2O3 + 0.75 O2 -> YMnO3
Y2O3 + 2 MnO + 0.5 O2 -> 2 YMnO3
Y2O3 + 0.1818 YMn12 + 1.773 O2 -> 2.182 YMnO3
Y2O3 + 0.6667 Mn3O4 + 0.3333 Y2Mn2O7 -> 2.667 YMnO3
Mn + 0.5 Y2O3 + 1.5 Y2Mn2O7 -> 4 YMnO3
Y2O3 + 2 MnO + Y2Mn2O7 -> 4 YMnO3
Y2O3 + 0.1818 YMn12 + 3.545 Y2Mn2O7 -> 9.273 YMnO3

This may be useful for modeling systems with 2 or more gaseous, liquid, or molten phases, such as \(O_2\) and a molten salt (e.g. \(LiCl\)).

If the system is open to a particular element such as \(O_2\), then it is generally a goood idea to model the thermodynamics of the system with the grand potential energy, \(\Phi\). This can be easily done by providing the open element and chemical potential to the ReactionSet class. In this case, we will reinitialize the reactions:

[27]:

grand_potential_rxns = ymno3_rxns_open2.set_chempot(open_el=Element("O"), chempot=0.0)
pprint(list(grand_potential_rxns)[:10])

[0.5 Y2Mn2O7 + 0.3333 Y -> YMnO3 + 0.1667 Y2O3 (mu_O=0.0),
 O2 + 0.1026 YMn12 -> 0.5641 Mn2O3 + 0.1026 YMnO3 (mu_O=0.0),
 Y2O3 + 0.6667 Mn3O4 + 0.1667 O2 -> 2 YMnO3 (mu_O=0.0),
 Mn + 0.5 Y2O3 + 0.75 O2 -> YMnO3 (mu_O=0.0),
 Y2O3 + 2 MnO + 0.5 O2 -> 2 YMnO3 (mu_O=0.0),
 Y2O3 + 0.1818 YMn12 + 1.773 O2 -> 2.182 YMnO3 (mu_O=0.0),
 Y2O3 + 0.6667 Mn3O4 + 0.3333 Y2Mn2O7 -> 2.667 YMnO3 (mu_O=0.0),
 Mn + 0.5 Y2O3 + 1.5 Y2Mn2O7 -> 4 YMnO3 (mu_O=0.0),
 Y2O3 + 2 MnO + Y2Mn2O7 -> 4 YMnO3 (mu_O=0.0),
 Y2O3 + 0.1818 YMn12 + 3.545 Y2Mn2O7 -> 9.273 YMnO3 (mu_O=0.0)]

Each of these reactions is a new class: OpenComputedReaction and its energy (per atom) is calculated assumning oxygen is an open element with the defined chemical potential. We will discuss this further later in the notebook.

b) Minimize enumerators#

The “minimize” enumerators produce reactions via a thermodynamic free energy minimization approach, rather than a purely combinatorial one. This means that reactions are produced directly from the compositional phase diagram, where a new convex hull is drawn connecting two compositions within a closed (Gibbs) or open (Grand Potential) system. See the InterfacialReactivity class within the pymatgen package for more information.

It is important to note that reactions produced with the minimize enumerators may overlap some with the basic enumerators, but the minimize enumerators have the restriction that all reactions they originally produce must have a negative reaction energy and result in a set of product phases which are stable with respect to each other (i.e. they share a facet of the phase diagram). That being said, these enumerators (when unrestricted) will also supply the reverse (i.e. positive energy) reaction as well.

In general, this means that the basic enumerators and minimize enumerators do not perfectly overlap in their outputs, and may be used in tandem. This is recommended! Let’s initialize the MinimizeGibbsEnumerator with all default arguments (press Shift+Tab to view documentation).

[28]:

mge = MinimizeGibbsEnumerator()

The default arguments, as before, help produce all reactions in a set of entries, given minimal constraints.

[29]:

rxns = mge.enumerate(filtered_entries)

Building phase diagrams (MinimizeGibbsEnumerator): 100%|█| 4/4 [00:00<00:00,
Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1240.73it/s]
Enumerating reactions (MinimizeGibbsEnumerator): 100%|█| 1/1 [00:00<00:00,  4

[30]:

for r in list(rxns)[0:10]:  # first 10 reactions
    print(r)

0.5 Mn2O3 + 0.5 Y2Mn2O7 -> YMn2O5
0.5 Mn5O8 -> Mn2O3 + 0.5 MnO2
0.5 Mn2O3 + 0.5 Y2O3 -> YMnO3
0.5 O2 + 2 YMnO3 -> Y2Mn2O7
MnO2 + YMnO3 -> YMn2O5
0.5 Mn5O8 -> Mn2O3 + 0.5 MnO2
2 MnO2 + Y2O3 -> Y2Mn2O7
0.5 Mn5O8 -> Mn2O3 + 0.5 MnO2
0.5 Mn5O8 -> Mn2O3 + 0.5 MnO2
0.5 Mn5O8 -> Mn2O3 + 0.5 MnO2

And as before, we can specify various combinations of precursors and targets, as well as whether or not they should be “exclusive”. Here we specify Y2O3 as a required, but non-exclusive precursor):

[31]:

mge_precursors = MinimizeGibbsEnumerator(precursors=["Y2O3"], exclusive_precursors=False)
rxns = mge_precursors.enumerate(filtered_entries)

Building phase diagrams (MinimizeGibbsEnumerator): 100%|█| 4/4 [00:00<00:00,
Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1344.54it/s]
Enumerating reactions (MinimizeGibbsEnumerator): 100%|█| 1/1 [00:00<00:00, 34

[32]:

for r in rxns:
    print(r)

2.5 Mn2O3 + 0.5 Y2O3 -> YMn2O5 + Mn3O4
0.5 Mn3O4 + 0.5 Y2O3 -> YMnO3 + 0.5 MnO
2 MnO2 + 0.5 Y2O3 -> YMn2O5 + 0.25 O2
0.3333 Mn5O8 + 0.6667 Y2O3 -> YMnO3 + 0.3333 YMn2O5
0.625 Mn5O8 + 0.5 Y2O3 -> YMn2O5 + 0.375 Mn3O4
0.6667 Mn5O8 + 0.3333 Y2O3 -> Mn2O3 + 0.6667 YMn2O5
0.5 Mn2O3 + 0.5 Y2O3 -> YMnO3
2 MnO2 + Y2O3 -> Y2Mn2O7

And here we specify both Y2O3 and Mn2O3 as the exclusive precursors of the reaction:

[33]:

mge_precursors = MinimizeGibbsEnumerator(precursors=["Y2O3", "Mn2O3"])
rxns = mge_precursors.enumerate(filtered_entries)

Building phase diagrams (MinimizeGibbsEnumerator): 100%|█| 4/4 [00:00<00:00,
Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1110.27it/s]
Enumerating reactions (MinimizeGibbsEnumerator): 100%|█| 1/1 [00:00<00:00, 81

[34]:

for r in rxns:
    print(r)

2.5 Mn2O3 + 0.5 Y2O3 -> YMn2O5 + Mn3O4
0.5 Mn2O3 + 0.5 Y2O3 -> YMnO3

Now we specify YMnO3 as the exclusive target of the reaction.

[35]:

mge_targets = MinimizeGibbsEnumerator(targets=["YMnO3"], exclusive_targets=True)
rxns = mge_targets.enumerate(filtered_entries)

Building phase diagrams (MinimizeGibbsEnumerator): 100%|█| 4/4 [00:00<00:00,
Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1311.03it/s]
Enumerating reactions (MinimizeGibbsEnumerator): 100%|█| 1/1 [00:00<00:00,  5

This only identifies one reaction!

[36]:

for r in rxns:
    print(r)

0.5 Mn2O3 + 0.5 Y2O3 -> YMnO3

Lastly, we identify YMnO3 as a required target, with any number of byproducts (i.e., non-exclusive target):

[37]:

mge_targets = MinimizeGibbsEnumerator(targets=["YMnO3"], exclusive_targets=False)
rxns = mge_targets.enumerate(filtered_entries)

Building phase diagrams (MinimizeGibbsEnumerator): 100%|█| 4/4 [00:00<00:00,
Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1345.41it/s]
Enumerating reactions (MinimizeGibbsEnumerator): 100%|█| 1/1 [00:00<00:00,  4

[38]:

for r in list(rxns)[:10]:  # first 10 reactions
    print(r, r.energy_per_atom)

Mn3O4 + 2 Y2Mn2O7 -> YMnO3 + 3 YMn2O5 -0.039591749392732606
0.5 Mn + 0.5 Y2Mn2O7 -> YMnO3 + 0.5 MnO -0.22628212150322882
0.375 Mn + 0.5 Y2Mn2O7 -> YMnO3 + 0.125 Mn3O4 -0.19916896885873445
0.3333 Mn + 0.6667 Y2Mn2O7 -> YMnO3 + 0.3333 YMn2O5 -0.15230582998188094
1.5 MnO + 0.5 Y2Mn2O7 -> YMnO3 + 0.5 Mn3O4 -0.07145795072025463
MnO + Y2Mn2O7 -> YMnO3 + YMn2O5 -0.06058835627265508
0.3333 Y + 0.5 Y2Mn2O7 -> YMnO3 + 0.1667 Y2O3 -0.4614723709416988
0.03448 YMn12 + 0.4828 Y2Mn2O7 -> YMnO3 + 0.3793 MnO -0.2475633372715506
0.02752 YMn12 + 0.4862 Y2Mn2O7 -> YMnO3 + 0.1009 Mn3O4 -0.22088724880664687
0.025 YMn12 + 0.625 Y2Mn2O7 -> YMnO3 + 0.275 YMn2O5 -0.17363615371992666

Notice that there are many more complex reactions here – all but the previous one we identified feature a byproduct phase that is required to balance the reaction.

Open entries#

Once again, as before, we can do all the same analysis with open entries. As previously mentioned, the grand potential is used as the thermodynamic free energy which is minimized:

\[\Phi = G - \mu_iN_i\]

Where \(i\) is the open species with chemical potential \(\mu_i\). The amount of the open species in a particular phase, \(N_i\), is used to adjust the energy of that phase.

Let’s create a MinimizeGrandPotentialEnumerator with default arguments and oxygen as the open element with chemical potential \(\mu=0\). Note that this is equivalent to \(\mu=\mu_O^0\) due to the fact that we are working with Gibbs free energies of formation:

[39]:

mgpe = MinimizeGrandPotentialEnumerator(open_elem=Element("O"), mu=0)
open_rxns = mgpe.enumerate(filtered_entries)

Building phase diagrams (MinimizeGrandPotentialEnumerator): 100%|█| 3/3 [00:0
Building chunks...: 100%|█████████████████████| 3/3 [00:00<00:00, 902.78it/s]
Enumerating reactions (MinimizeGrandPotentialEnumerator): 100%|█| 1/1 [00:00<

[40]:

for r in list(open_rxns)[:10]:  # first 10 reactions
    print(r)

2 YMnO3 + 0.5 O2 -> Y2Mn2O7
2 YMnO3 + 0.5 O2 -> Y2Mn2O7
0.5 Mn2O3 + 0.25 O2 -> MnO2
2 YMnO3 + 0.5 O2 -> Y2Mn2O7
0.3333 Mn3O4 + 0.3333 O2 -> MnO2
Mn + O2 -> MnO2
MnO + 0.5 O2 -> MnO2
0.2 Mn5O8 + 0.2 O2 -> MnO2
2 Y + 1.5 O2 -> Y2O3
0.5 Mn2O3 + 0.25 O2 -> MnO2

Note that this enumerator often returns duplicates particularly with terminal reactions.

To ensure that duplicates are not returned, we can initialize the enumerator with a different flag:

[41]:

mgpe = MinimizeGrandPotentialEnumerator(open_elem=Element("O"), mu=0, filter_duplicates=True)
open_rxns = mgpe.enumerate(filtered_entries)

for r in list(open_rxns)[:10]:  # first 10 reactions
    print(r)

Building phase diagrams (MinimizeGrandPotentialEnumerator): 100%|█| 3/3 [00:0
Building chunks...: 100%|█████████████████████| 3/3 [00:00<00:00, 894.44it/s]
Enumerating reactions (MinimizeGrandPotentialEnumerator): 100%|█| 1/1 [00:00<
Filtering duplicates: 100%|████████████████████| 8/8 [00:00<00:00, 22.67it/s]

0.3333 Mn3O4 + 0.3333 O2 -> MnO2
0.5 Mn2O3 + 0.25 O2 -> MnO2
0.2 Mn5O8 + 0.2 O2 -> MnO2
2 Y + 1.5 O2 -> Y2O3
Mn + O2 -> MnO2
MnO + 0.5 O2 -> MnO2
2 YMnO3 + 0.5 O2 -> Y2Mn2O7
2 MnO2 + 2 Y + 1.5 O2 -> Y2Mn2O7
2 MnO2 + Y + 0.5 O2 -> YMn2O5
4 MnO2 + 2 Y2Mn2O7 -> O2 + 4 YMn2O5

Also note that the reaction objects returned are NOT (by default) configured to report their energy as a change in the grand potential. To configure this, we need to transform these reactions to OpenComputedReaction objects. This class allows for easy specification of reactions where one of the elements is assigned a chemical potential.

We did this earlier through a convenience constructor method:

[42]:

new_open_rxns = open_rxns.set_chempot(Element("O"), chempot=0.0)

We should now see that (mu_O=0.0) appears in the repr of the object:

[43]:

sample_open_rxn = list(new_open_rxns)[0]
print(sample_open_rxn.__class__.__name__)

sample_open_rxn

OpenComputedReaction

[43]:

0.3333 Mn3O4 + 0.3333 O2 -> MnO2 (mu_O=0.0)

We can also, as before, customize which precursors and targets are specified:

[44]:

mgpe_precursors = MinimizeGrandPotentialEnumerator(open_elem=Element("O"), mu=0,
                                                   precursors=["Y2O3"], exclusive_precursors=False, filter_duplicates=True)
open_rxns_precursors = mgpe_precursors.enumerate(filtered_entries)

Building phase diagrams (MinimizeGrandPotentialEnumerator): 100%|█| 3/3 [00:0
Building chunks...: 100%|█████████████████████| 3/3 [00:00<00:00, 883.38it/s]
Enumerating reactions (MinimizeGrandPotentialEnumerator): 100%|█| 1/1 [00:00<
Filtering duplicates: 0it [00:00, ?it/s]

[45]:

for r in open_rxns_precursors:
    print(r)

8 MnO2 + 2 Y2O3 -> O2 + 4 YMn2O5
2 MnO + Y2O3 + O2 -> Y2Mn2O7
2 MnO + 0.5 Y2O3 + 0.75 O2 -> YMn2O5
0.4 Mn5O8 + Y2O3 + 0.4 O2 -> Y2Mn2O7
0.4 Mn5O8 + 0.5 Y2O3 + 0.15 O2 -> YMn2O5
0.6667 Mn3O4 + Y2O3 + 0.6667 O2 -> Y2Mn2O7
0.6667 Mn3O4 + 0.5 Y2O3 + 0.4167 O2 -> YMn2O5
2 Mn + Y2O3 + 2 O2 -> Y2Mn2O7
2 Mn + 0.5 Y2O3 + 1.75 O2 -> YMn2O5
Mn2O3 + Y2O3 + 0.5 O2 -> Y2Mn2O7
Mn2O3 + 0.5 Y2O3 + 0.25 O2 -> YMn2O5
0.1667 YMn12 + 0.9167 Y2O3 + 2.125 O2 -> Y2Mn2O7
0.1667 YMn12 + 0.4167 Y2O3 + 1.875 O2 -> YMn2O5
0.5 Y2O3 + YMn2O5 + 0.25 O2 -> Y2Mn2O7

[46]:

mgpe_precursors = MinimizeGrandPotentialEnumerator(open_elem=Element("O"), mu=0,
                                                   targets=["Y2Mn2O7"], filter_duplicates=True)
open_rxns_targets = mgpe_precursors.enumerate(filtered_entries)

Building phase diagrams (MinimizeGrandPotentialEnumerator): 100%|█| 3/3 [00:0
Building chunks...: 100%|█████████████████████| 3/3 [00:00<00:00, 915.79it/s]
Enumerating reactions (MinimizeGrandPotentialEnumerator): 100%|█| 1/1 [00:00<
Filtering duplicates: 100%|██████████████████| 1/1 [00:00<00:00, 1811.79it/s]

[47]:

for r in open_rxns_targets:
    print(r)

2 YMnO3 + 0.5 O2 -> Y2Mn2O7
2 MnO2 + 2 Y + 1.5 O2 -> Y2Mn2O7
2 MnO + 2 Y + 2.5 O2 -> Y2Mn2O7
2 MnO + Y2O3 + O2 -> Y2Mn2O7
0.4 Mn5O8 + 2 Y + 1.9 O2 -> Y2Mn2O7
0.4 Mn5O8 + Y2O3 + 0.4 O2 -> Y2Mn2O7
0.6667 Mn3O4 + 2 Y + 2.167 O2 -> Y2Mn2O7
0.6667 Mn3O4 + Y2O3 + 0.6667 O2 -> Y2Mn2O7
2 Mn + 2 Y + 3.5 O2 -> Y2Mn2O7
2 Mn + Y2O3 + 2 O2 -> Y2Mn2O7
Mn2O3 + 2 Y + 2 O2 -> Y2Mn2O7
Mn2O3 + Y2O3 + 0.5 O2 -> Y2Mn2O7
0.1667 YMn12 + 1.833 Y + 3.5 O2 -> Y2Mn2O7
0.1667 YMn12 + 0.9167 Y2O3 + 2.125 O2 -> Y2Mn2O7
Y + YMn2O5 + O2 -> Y2Mn2O7
0.5 Y2O3 + YMn2O5 + 0.25 O2 -> Y2Mn2O7

Note that setting the chemical potential to a value outside of the range of stability of the target (e.g., mu_O = -3, causes the enumerator to yield no reactions:

[48]:

mgpe_precursors = MinimizeGrandPotentialEnumerator(open_elem=Element("O"), mu=-3,
                                                   targets=["Y2Mn2O7"])
open_rxns_targets = mgpe_precursors.enumerate(filtered_entries)

Building phase diagrams (MinimizeGrandPotentialEnumerator): 100%|█| 3/3 [00:0
Building chunks...: 100%|█████████████████████| 3/3 [00:00<00:00, 929.04it/s]
Enumerating reactions (MinimizeGrandPotentialEnumerator): 100%|█| 1/1 [00:00<

[49]:

print(len(open_rxns_targets))

0

This concludes the introduction of the four enumerator classes. In general, the configurations we showed are all that you will really need to use the reaction-network package. You are free to explore the other enumerator arguments; note that many of these often do not need to be changed.

It is encouraged for you to try out what you have learned above on your own system of interest!

4. Running enumerators with the jobflow package#

Running reaction enumeration calculations in high-throughput can be challenging. Often, we want to be able to enumerate reactions as part of some workflow and keep track of them using a database structure.

For this reason, we turn to the jobflow package, which is a fantastic toolit for developing and running computational workflows. The reaction-network package has several jobs (and flows) written using jobflow that can be either run directly or strung together as part of custom workflows.

Let’s run one of the jobs via ReactionEnumerationMaker locally on your computer. This job does exactly what we were doing above – run a list of enumerators on a provided entry set and collect the reactions in a ReactionSet. First, we import the job maker and run function:

[50]:

from jobflow.managers.local import run_locally
from rxn_network.jobs.core import ReactionEnumerationMaker

Now, all we need to do is initialize the Maker and use the make() function with a list of enumerators and our entry set. The enumerators can be customized when they are created and this will be passed to the jbo:

[51]:

be = BasicEnumerator()  # put your desired settings here
boe = BasicOpenEnumerator(open_phases=["O2"])

maker = ReactionEnumerationMaker()
job = maker.make([be, boe], filtered_entries)

The job can now be run either locally (as shown here) or launched on a remote workstation using the fireworks package.

Ti configure launching these jobs on other machines, or how the data is stored in various databases, then please see the jobflow documentation here: https://materialsproject.github.io/jobflow/.

WARNING: You must have jobflow properly configured in order to run this successfully (see link above). This means configuring your jobflow.yaml file and having a working connection to the database specified. For this particular job, you should also specify an additional data store called “rxns” in your jobflow.yaml file. I typically make this a GridFSStore as the outputs can get quite large.

[52]:

output = run_locally(job)

2023-08-31 11:22:25,150 INFO Started executing jobs locally
2023-08-31 11:22:25,152 INFO Starting job - enumerate_reactions (9471b6cb-ace1-40db-bb96-75dc155c00bd)
2023-08-31 11:22:25,154 INFO rxn_network.jobs.core Running enumerators...
2023-08-31 11:22:25,154 INFO rxn_network.enumerators.utils Running BasicEnumerator

Building chunks...: 100%|████████████████████| 4/4 [00:00<00:00, 1892.95it/s]
Enumerating reactions (BasicEnumerator): 100%|█| 1/1 [00:00<00:00,  1.96it/s]

2023-08-31 11:22:25,677 INFO rxn_network.enumerators.utils Adding 692 reactions to reaction set
2023-08-31 11:22:25,678 INFO rxn_network.enumerators.utils Running BasicOpenEnumerator


Building chunks...: 100%|████████████████████| 3/3 [00:00<00:00, 1213.63it/s]
Enumerating reactions (BasicOpenEnumerator): 100%|█| 2/2 [00:00<00:00,  2.16i

2023-08-31 11:22:26,627 INFO rxn_network.enumerators.utils Adding 204 reactions to reaction set
2023-08-31 11:22:26,628 INFO rxn_network.enumerators.utils Completed reaction enumeration. Filtering duplicates...


Filtering duplicates: 100%|██████████████████| 22/22 [00:00<00:00, 61.66it/s]

2023-08-31 11:22:26,999 INFO rxn_network.enumerators.utils Completed duplicate filtering.
2023-08-31 11:22:26,999 INFO rxn_network.jobs.core Building task document...



2023-08-31 11:22:27,501 INFO Finished job - enumerate_reactions (9471b6cb-ace1-40db-bb96-75dc155c00bd)
2023-08-31 11:22:27,502 INFO Finished executing jobs locally

When the job completes, it will store its data in the JobStore database. However, the output reactions can also be accessed locally here in the notebook. To do this, access the rxns property of the output EnumeratorTaskDocument:

[53]:

response = output[job.uuid][1].output
response.dict()

[53]:

{'task_label': 'enumerate_reactions',
 'last_updated': '2023-08-31 18:22:26.999828',
 'rxns': <rxn_network.reactions.reaction_set.ReactionSet at 0x34617bdf0>,
 'targets': [],
 'elements': [Element Y, Element Mn, Element O],
 'chemsys': 'Y-Mn-O',
 'added_elements': None,
 'added_chemsys': None,
 'enumerators': [<rxn_network.enumerators.basic.BasicEnumerator at 0x3410e38e0>,
  <rxn_network.enumerators.basic.BasicOpenEnumerator at 0x3410e23e0>]}

[54]:

for r in list(response.rxns)[0:10]:
    print(r)

0.5 Mn2O3 -> MnO + 0.25 O2
MnO + 0.25 O2 -> 0.5 Mn2O3
0.6 Mn5O8 -> Mn3O4 + 0.4 O2
Mn3O4 + 0.4 O2 -> 0.6 Mn5O8
Mn + 0.8 O2 -> 0.2 Mn5O8
0.2 Mn5O8 -> Mn + 0.8 O2
0.4 Mn5O8 -> Mn2O3 + 0.1 O2
Mn2O3 + 0.1 O2 -> 0.4 Mn5O8
0.5 Y2Mn2O7 -> YMnO3 + 0.25 O2
YMnO3 + 0.25 O2 -> 0.5 Y2Mn2O7

5. Conclusion#

That’s it! We hope this notebook was helpful in introducing the enumerator classes found in the reaction-network package. If any significant errors are encountered, please first double-check that your settings are configured properly (e.g., proper installation of all dependencies, configuration of Materials Project API key, etc.).

If the error persists, then please raise an Issue here: https://github.com/materialsproject/reaction-network/issues.

Happy enumerating!