rn.build() on MGPE-enumerated reactions

[ancarnevali]

Hi!
I am trying to construct a reaction network after enumerating reactions with MinimizeGrandPotentialEnumerator, but am encountering some probelms when running find_pathways.

with MPRester(key) as mpr:
    entries = mpr.get_entries_in_chemsys("Fe-O-Nb")
temp = 1073
gibbs_entries = GibbsEntrySet.from_computed_entries(entries, temp)
tresh = 0.005
entry_set = gibbs_entries.filter_by_stability(tresh)
###
mgpe = minimizeGrandPotentialEnumerator(
    open_elem=Element("O"), mu=0, filter_duplicates=True
)
rxns = mgpe.enumerate(entry_set)
grand_potential_rxns = rxns.set_chempot(open_el=Element("O"), chempot=0.0)
###
cf = Softplus(temp)
rn = ReactionNetwork(grand_potential_rxns, cf) # should I be passing non `set_chempot`ted reactions instead?
rn.build()
rn.set_precursors(["Fe2O3", "NbO"])
rn.set_target("Nb2FeO6")
paths = rn.find_pathways(["Nb2FeO6"])

This returns no paths for the target compound.
On a sidenote, plotting the rn shows a rather expected graph, too. This might indicate that, if there is an issue, this might be in the find_pathways method?

Running the same script with the basic open enumerator be = BasicOpenEnumerator(open_phases=["O2"], filter_duplicates=True) has no issues though, and returns reaction paths as expected.

Is there something I am missing about the use of MGPE for building the reaction network, or what do you think is going on here?

Thank you!

[mattmcdermott]

Hi @ancarnevali

The issue you're running into with MinimizeGrandPotentialEnumerator is because the Nb2FeO6 target is not a stable phase at $\mu_O = 0$ eV. To confirm this, you an look at the grand potential phase diagram:

from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram
GrandPotentialPhaseDiagram(entry_set, {Element("O"):0})

Which shows:
Nb-Fe GrandPotentialPhaseDiagram with chempots = 'mu_O = 0.0000'3 stable phases: NbFeO4, Fe2O3, Nb2O5

Using a chemical potential around -2 eV should stabilize the phase and make the rest of your code produce results.

In general, for open systems, I recommend using all four enumerators together: BasicEnumerator, BasicOpenEnumerator, MinimizeGibbsEnumerator and MinimizeGrandPotentialEnumerator. Combine their resulting reactions into a new ReactionSet and then set the chemical potential to the one of interest. This tends to be more comprehensive in capturing all possible reactions. Hope this helps!

[ancarnevali]

I see, thank you for your help!
Actually, at this point I'm not sure how the mu_0 parameter is to be interpreted though (and how it differs from the mu that is to be passed to the MGPE instantiation). For example, what would it mean if a compound is found stable at mu_0 = -5, but not at a less negative value? What in the opposite case?
I am guessing it is the chemical potential of the open phase, but how does it affect the costruction of the GrandPotentialPhaseDiagram? Also in the A graph-based network for predicting chemicalreaction pathways in solid-state materialssynthesis paper, I am struggling to find an explicit explanation of (although I totally might just be not finding it).

More in general, it seems like a good way to use the package would be to first find the chemical potential at which the target compound is stable at and then proceed with enumerating and constructing the network. Is there a proper way to do so, or would you just iterate over a range of chempots and inspect the stable phases?