Hi,
I have a few questions regarding pymatgen’s pourbaix_diagram module. First, I don’t know what the meaning of nPhi vs nH2O is in a PourbaixEntry. For example, when I use find_stable_entry(pH,V), the output is something like:
Multiple Pourbaix Entry: energy = 4.0981, npH = -3.3999999999999977, nPhi = -3.199999999999998, nH2O = 1.6999999999999988, entry_id = [‘mp-1180536’, ‘ion-11’, ‘mp-1094139’], species: NaO8(s) + FeO4[2-] + Ni2O5(s)
According to the documentation here, nPhi and nH2O both represent the number of H2O. Since these values are different, I don’t understand what nPhi and nH2O are actually telling me.
Additionally, in the above output from find_stable_entry(pH,V), it gives me an energy. I’m having trouble understanding how this energy is calculated and how exactly to interpret the energy. If it is a free energy, and this is a ‘stable’ entry, I would expect the energy to be zero? So this then begs the question, what defines ‘stability’ in a stable entry?
The larger context for the problem I am trying to solve is: I have a list of materials and their associated mp-id. Out of this list, I would like to find materials that do not interact with H2O at a given potential and pH. I’m still not quite sure how to accomplish this with Pourbaix, so any tips are welcome! Thank you for any responses.
Hi @mdresser , find_stable_entry tells you which materials are stable in a given Pourbaix domain, based on the Pourbaix energy (which is based on a Legendre transformation of the free energy of formation for each material). In this case, you get multiple stable phases, indicating that all three may coexist or compete under those conditions. You can find more about the theory behind the Pourbaix energy in the Methods section of this work
As to nH2O vs nPhi, it looks like the documentation of nPhi is inaccurate (thanks for reporting!). If you look at the source code, nPhi is defined as
self.npH - self.charge
Which you could think of as the number of excess protons (not associated with H2O) minus the overall charge of the material, or in other words, the number of electrons involved in a redox reaction with that species.
Thank you so much for your response!
That paper was quite helpful and made me realize in order to accomplish my task I most certainly care about the electrochemical stability (which seems to be an additional metric that can be added to the Pourbaix diagram). I see the method plot_entry_stability is used to add the electrochemical stability of a specific entry to a Pourbaix plot. However, I want to perform a more high throughput search by simply extracting deltaGpbx (in the language of the paper you linked) for a specific material at a given pH and voltage rather than trying to look at the plot for each material. It isn’t clear to me from the documentation what method could be used to accomplish this. For example, I see the method normalized_energy_at_conditions(pH , V) but is this the “Pourbaix Energy” you mentioned above, or is it something like the deltaGpbx the paper focuses on? Or are these energies equivalent and perhaps I just misunderstood the paper? Thank you for any clarification!
I believe G_pbx is the same thing as the Pourbaix energy, and you can use normalized_energy_at_conditions to get it.
plot_entry_stability should show you the relative stability of a given entry relative to the stable phases at every location on the Pourbaix diagram. So if it’s stability is zero, then the entry you’re examining is the stable phase, but if it’s a positive number that means it’s unstable (by that amount) relative to the stable phases.
Hope this helps!
Hi, thank you for your response! After looking further at the paper and documentation, I think I need to use get_decomposition_energy(entry , pH , V). However, when I try to use this I get an error. Here is my code:
from pymatgen import MPRester
import pymatgen.analysis.pourbaix_diagram as pb
els = ['O', 'Na', 'Si', 'Al']
mp_id = "mp-560334"
comp_dict = {'Na': 0.14285714285714285, 'Al': 0.14285714285714285, 'Si': 0.14285714285714285, 'O': 0.5714285714285714}
with MPRester(api_key) as m:
pe = m.get_pourbaix_entries(els)
dia = pb.PourbaixDiagram(pe, comp_dict)
mpid_entry = [entry for entry in pe if entry.entry_id==mp_id][0]
decomp_e = dia.get_decomposition_energy(mpid_entry, pH, V)
When I run this, I get the following error:
ValueError: Composition of stability entry does not match Pourbaix Diagram
Upon looking further, I see that this particular material is included in the “unprocessed” entries:
unproc = [entry for entry in dia.unprocessed_entries if entry.entry_id==mp_id][0]
print("unproccessed",unproc)
This indeed gives the output:
unproccessed Pourbaix Entry : Na8 Al8 Si8 O32 with energy = -88.2080, npH = -64.0, nPhi = -64.0, nH2O = 32.0, entry_id = mp-560334
Presumably this error occurs because this entry is “unprocessed”? Can you give me any information on how I could fix this error and why this entry is “unprocessed”? Thank you!
Hi @mdresser , I’m unfamiliar with this error, but after a bit of code diving I think the problem is that you have included O in your comp_dict. Since O and H are intrinsic parts of all Pourbaix diagrams, these elements are not considered when the “composition” of the Pourbaix diagram is defined. See this section of code. In this case, the code checks that comp_dict has the same composition as the diagram as a whole, which is not the case since O is excluded from the diagram’s composition.
I’d have to think more about whether this behavior makes sense or is just a necessary workaround based on how we have to transform the energies. In any case, give it a try by defining your comp_dict based on the non-OH elements only.
Thank you! This fixed the issue.