Inconsistency in Ion Energies for C-Containing Entries in Pymatgen Pourbaix Diagrams

While constructing Pourbaix diagrams (Ti-C system) using Pymatgen, I noticed something strange about the ion energies in the Pymatgen dictionary. The way ion energies are calculated in the code (pymatgen/pymatgen/ext/ at a1d508844dadbae6d930b66ab88bb13df52bf538 · materialsproject/pymatgen · GitHub) seems logical and consistent for some ions like Ti²⁺ and Ti³⁺. However, when it comes to ions containing carbon (C), the values in the dictionary look unusual.

I understand that ion energy is calculated as following.

Ion_energy = (ion_energy_exp ) + [ref_solid_energy_DFT - ref_solid_energy_exp]* (cf) where cf is ratio of major element in ion ion to major element in ref solid.

This implies that for the single C-containing entries (such as CH4, CO3-2, CH3OH etc.,) the difference between the dictionary value and the experimental value (available at MPContribs) for each entry is supposed to be a constant. The correction is applied for any C-containing ions using Na2CO3 as the reference solid. So, as long as the ions contain the same amount of C in their formula, the energy difference between the dictionary and experimental values should be the same. But it is not the case.

Could you please highlight the possible reasons for this discrepancy? Thanks.
(Energy in eV/formula)



Please let me know what the “Dict Value” means (which you said is available at MPContribs).
I’m asking this as the only value available at the MPContribs is the experimental value for this ion and its reference solid. Also, how did you get the DFT values for these Ions?

Please refer to this paper on how these ion energies are calculated : Phys. Rev. B 85, 235438 (2012) - Prediction of solid-aqueous equilibria: Scheme to combine first-principles calculations of solids with experimental aqueous states

Hi @Arpit_Bhardwaj , thanks for the response. I will elaborate the issue here in more detail. If I run the code given below (code to makes Pourbaix Diagram) and print “entries”, I get a list of entries to make Pourbaix diagram (dict values). Here are a few values.

[PourbaixEntry(Ti1 with energy = -4.2458, npH = 0, nPhi = -2.0, nH2O = 0, entry_id = 'ion-0'),

PourbaixEntry(C1 H4 with energy = -1.3085, npH = 4.0, nPhi = 4.0, nH2O = 0, entry_id = 'ion-15'),

PourbaixEntry(C1 O3 with energy = 0.9529, npH = -6.0, nPhi = -4.0, nH2O = 3.0, entry_id = 'ion-8')

Now referring to the section “G. Aqueous ions” of the article you shared, energy values should have been calculated by following formula.

Ion_energy = (ion_energy_exp ) + [ref_solid_energy_DFT - ref_solid_energy_exp]

The values for the first and third term of left-hand side of the equation can be obtained from MPContribs whereas value for second term comes from

I made a few calculations (screenshot attached), for Ti2+ (ion-0 here) calculated values was quite close to what code printed, but it was not the case for C containing ions.

> from mp_api.client import MPRester
> from pymatgen.analysis.pourbaix_diagram import PourbaixDiagram, PourbaixPlotter, PDEntry, PourbaixEntry
> %matplotlib inline
> # Initialize the MP Rester
> mpr = MPRester(MP_API_KEY)
> # Get all pourbaix entries corresponding to the chemical system.
> entries = mpr.get_pourbaix_entries(["Ti","C"])
> print (entries)

Hi @nihal , I’m not sure I follow. I would not say that your calculated values for Ti+2 are close - -3.96 vs. -4.24 eV is a big difference. If you are following the procedure correctly you should have exactly the right number there.

I think the issue is with your ref_solid_energy_DFT The .energy attribute of the PourbaixEntry is the formation energy of your entire structure (Ti, CH4, or CO3), in eV per formula unit, corrected for concentration (see here). I suspect that in your manual calculations you haven’t accounted for the concentration correction.

You might want to examine the source of get_ion_entries here to understand the steps.

If you go through the above and still think there is a bug, please let us know!