[Bug]: Phase Diagram Failures with Mixed Entries

[kavanase]

Code snippet

from mp_api.client import MPRester
from pymatgen.analysis.phase_diagram import PhaseDiagram

with MPRester() as mpr:
    entries = mpr.get_entries_in_chemsys(["Cs", "Ti", "I"])  # default GGA_GGA+U_R2SCAN
pd = PhaseDiagram(entries)
host = next(e for e in entries if e.composition.reduced_formula == "Cs2TiI6")
print(pd.get_e_above_hull(host))  # ~4.574  (expected: 0, on the hull, see https://next-gen.materialsproject.org/phasediagram)

What happened?

I have found that with the (relatively-recent?) default of using the mixed GGA-R2SCAN entries as the return type with get_entries_in_chemsys, the functional-mixing energy correction is incorrectly applied, resulting in several phases being incorrectly predicted as unstable and vice versa.

For the example of Cs2TiI6, which is on the convex hull for GGA, R2SCAN and mixed on the MP database / phase diagrams (https://next-gen.materialsproject.org/phasediagram), this is correctly returned when specifying GGA or R2SCAN thermo types only:

from mp_api.client import MPRester
from pymatgen.analysis.phase_diagram import PhaseDiagram

with MPRester() as mpr:
    entries = mpr.get_entries_in_chemsys(["Cs", "Ti", "I"], additional_criteria={"thermo_types": ["GGA_GGA+U"]})
# or additional_criteria={"thermo_types": ["R2SCAN]}, same result
pd = PhaseDiagram(entries)
host = next(e for e in entries if e.composition.reduced_formula == "Cs2TiI6")
print(pd.get_e_above_hull(host))  # 0.0
-> 0

But with the default mixed thermo_type GGA_GGA+U_R2SCAN as shown in the MWE, it returns a spurious ~4.6 eV/atom above-the-hull energy.

This seems to be related to the cached mixing correction on the mixed-functional entries. When I run the following code which re-applies mixing corrections, it all works fine with the expected result:

from mp_api.client import MPRester
from pymatgen.analysis.phase_diagram import PhaseDiagram
from pymatgen.entries.mixing_scheme import MaterialsProjectDFTMixingScheme

with MPRester() as mpr:
    gga = mpr.get_entries_in_chemsys(
        ["Cs", "Ti", "I"], additional_criteria={"thermo_types": ["GGA_GGA+U"]},
    )
    r2 = mpr.get_entries_in_chemsys(
        ["Cs", "Ti", "I"], additional_criteria={"thermo_types": ["R2SCAN"]},
    )

# Re-apply the GGA(+U)/r2SCAN mixing scheme on the combined full-chemsys raw entries:
entries = MaterialsProjectDFTMixingScheme(run_type_2="r2SCAN").process_entries(gga + r2)

pd = PhaseDiagram(entries)
host = next(e for e in entries if e.composition.reduced_formula == "Cs2TiI6")
print(pd.get_e_above_hull(host))  # 0.0  (vs ~4.574 from the default chained call)

print(host.name)
print(host.energy_adjustments)
print(host.energy_per_atom)
~~~~~
Cs2TiI6
[]
-24.70879382444444

While printing the same after the MWE which fails shows:

Cs2TiI6
[ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=1.91e+02, uncertainty=0.0, description='Place r2SCAN energy onto the GGA(+U) hull', generated_by='MaterialsProjectDFTMixingScheme')]
-3.4718101466666678

For reference, printing the energy adjustments for all entries here after re-applying corrections gives:

for entry in entries:
    print(entry.name, entry.energy_adjustments)
~~~~~~
CsI [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-0.379, uncertainty=0.0055, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-47.3, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
TiI2 [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-0.758, uncertainty=0.011, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-50.7, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Cs [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-1.94e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
I [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-23.0, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Cs [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-96.9, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-5.04e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
CsI [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-1.52, uncertainty=0.022, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-1.89e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
I [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-23.0, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
I [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-23.0, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
TiI [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-0.758, uncertainty=0.011, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-55.8, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-5.04e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-5.04e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Cs [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-24.2, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-5.04e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-30.3, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti [ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-5.04e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
CsTi2I7 [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-10.6, uncertainty=0.154, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-7.78e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Cs3I [CompositionEnergyAdjustment(name='MP2020 anion correction (I)', value=-0.758, uncertainty=0.011, description='Composition-based energy adjustment', generated_by='MaterialsProject2020Compatibility'), ConstantEnergyAdjustment(name='MP GGA(+U)/r2SCAN mixing adjustment', value=-1.91e+02, uncertainty=0.0, description='Place GGA(+U) energy onto the r2SCAN hull', generated_by='MaterialsProjectDFTMixingScheme')]
Ti []
CsI []
Cs []
Cs []
TiI3 []
CsI3 []
Cs []
Cs []
Cs2TiI6 []
CsI []
Cs []
Ti []
CsTiI3 []
Cs []
I []
Ti []
Cs []
Cs []
CsI []
Cs []
TiI3 []
TiI4 []
CsI4 []
Cs []
Cs []
Cs []
Cs []
Cs []
I []
Ti []

Version

0.46.1

Which OS?

• MacOS
• Windows
• Linux

Log output

[kavanase]

Minor point to note btw; the docstring still states that GGA(+U) is the default thermo-type, but it's now actually mixed GGA-R2SCAN:

api/mp_api/client/mprester.py

Line 1062 in 64c0ae9

Note that by default this returns mixed GGA/GGA+U entries. For others,

[tsmathis]

Thanks for the report @kavanase, tangential but just noting that the second approach you listed of retrieving the gga and r2scan entries and applying the mixing scheme is functionally identical for what MP does to produce phase diagrams during the data release process. Using the phase diagram retrieval function from the client (and following your process above) gets me:

>>> with MPRester() as mpr:
...     mixed_pd = mpr.materials.thermo.get_phase_diagram_from_chemsys("Cs-Ti-I", thermo_type=ThermoType.GGA_GGA_U_R2SCAN)
...     host = next(e for e in mixed_pd.all_entries if e.composition.reduced_formula == "Cs2TiI6")
...     print(mixed_pd.get_e_above_hull(host))
0.0

So thankfully the phase diagrams served by the client and website should be sane.

That said, we have seen weird things with the mixing scheme from time to time, @esoteric-ephemera has been investigating more thoroughly so I'll leave the rest for him to add in when he has time.

[kavanase]

@tsmathis thanks for the response!

Yes I was assuming the MP website approach must be different given that it displays the correct energies, good to know it's functionally the same as the alternative re-calculation approach here.

[kavanase]

Another minor point just to flag from looking through some of this code; the default stated in the docstring for thermo_type doesn't seem to match the code for this function:

api/mp_api/client/routes/materials/thermo.py

Lines 145 to 153 in 64c0ae9

	def get_phase_diagram_from_chemsys(
	self, chemsys: str, thermo_type: ThermoType | str
	) -> PhaseDiagram:
	"""Get a pre-computed phase diagram for a given chemsys.
	Arguments:
	chemsys (str): A chemical system (e.g. Li-Fe-O)
	thermo_type (ThermoType): The thermo type for the phase diagram.
	Defaults to ThermoType.GGA_GGA_U.

[tsmathis]

Ah yes, good catch. We'll make sure to fix that along with the docstr issue you mentioned above for get_entries_in_chemsys for the next mp-api release.

[esoteric-ephemera]

Thanks @kavanase! Have added those docstr changes here which will be in the forthcoming release like @tsmathis said

Also noting I changed the default thermotype to be consistent across utilities - can revert if preferred