E_hull prediction and E_form prediction with code build_phase_diagram.py

[RyotaroOKabe]

I would like to use Matbench-discovery to evaluate the stability of materials. I am trying to calculate E_hull, but I have trouble doing that.

I saw build_phase_diagram.py present how PatchedPhaseDiagram object ppd_mp is defined, and formation energy is calculated.
I tried the source code, but the predicted E_form was not accurate as written in the preprint. I attached the scatter plot and the code I used. I used almost the same code as the default build_phase_diagram.py, but I commented out several lines (I marked the lines with "#!" with the error message.). I pasted my current code of build_phase_diagram.py below. It would be nice if any of you knew how I can get the correct prediction of E_hull and E_form.

"""Build a PatchedPhaseDiagram from all MP ComputedStructureEntries for calculating
DFT-ground truth convex hull energies.
"""


# %%
import gzip
import json
import os
import pickle

import pandas as pd
import pymatviz
from pymatgen.analysis.phase_diagram import PatchedPhaseDiagram
from pymatgen.entries.compatibility import MaterialsProject2020Compatibility
from pymatgen.entries.computed_entries import ComputedEntry, ComputedStructureEntry
from pymatgen.ext.matproj import MPRester
from tqdm import tqdm

from matbench_discovery import MP_DIR, ROOT, Key, today
from matbench_discovery.data import DATA_FILES
from matbench_discovery.energy import get_e_form_per_atom, get_elemental_ref_entries

module_dir = os.path.dirname(__file__)


# %% run on 2022-09-16 and again on 2023-02-07
all_mp_computed_structure_entries = MPRester().get_entries("")

# save all ComputedStructureEntries to disk
# mp-15590 appears twice so we drop_duplicates()
df = pd.DataFrame(all_mp_computed_structure_entries, columns=["entry"])
df.index.name = Key.mat_id
df.index = [e.entry_id for e in df.entry]
df.reset_index().to_json(
    f"{module_dir}/{today}-mp-computed-structure-entries.json.gz",
    default_handler=lambda x: x.as_dict(),
)


# %%
# data_path = f"{module_dir}/2023-02-07-mp-computed-structure-entries.json.gz"
data_path = f"{module_dir}/2024-03-05-mp-computed-structure-entries.json.gz"    #! specify the date in the filename.
# df = pd.read_json(data_path).set_index(Key.mat_id)
df = pd.read_json(data_path)    #.set_index(Key.mat_id) #!  KeyError: "None of [<Key.mat_id: 'material_id'>] are in the columns"

# drop the structure, just load ComputedEntry, makes the PPD faster to build and load
mp_computed_entries = [ComputedEntry.from_dict(dct) for dct in tqdm(df.entry)]

print(f"{len(mp_computed_entries) = :,} on {today}")
# len(mp_computed_entries) = 146,323 on 2022-09-16
# len(mp_computed_entries) = 154,719 on 2023-02-07


# %% build phase diagram with MP entries only
ppd_mp = PatchedPhaseDiagram(mp_computed_entries, verbose=True)
print(f"{ppd_mp} on {today}")
# prints:
# PatchedPhaseDiagram covering 44805 sub-spaces on 2022-09-16
# PatchedPhaseDiagram covering 46216 sub-spaces on 2023-02-07

# save MP PPD to disk
# with gzip.open(f"{module_dir}/{today}-ppd-mp.pkl.gz", "wb") as zip_file:  #! AttributeError: Can't pickle local object 'cached_class.<locals>._decorated'
#     pickle.dump(ppd_mp, zip_file)


# %% build phase diagram with both MP entries + WBM entries
df_wbm = pd.read_json(DATA_FILES.wbm_computed_structure_entries).set_index(Key.mat_id)

# using ComputedStructureEntry vs ComputedEntry here is important as CSEs receive
# more accurate energy corrections that take into account peroxide/superoxide nature
# of materials (and same for sulfides) based on atomic distances in the structure
# wbm_computed_entries: list[ComputedStructureEntry] = df_wbm.query(    #! UndefinedVariableError: name 'n_elements' is not defined
#     "n_elements > 1"
# ).cse.map(ComputedStructureEntry.from_dict)
wbm_computed_entries: list[ComputedStructureEntry] = df_wbm.cse.map(ComputedStructureEntry.from_dict)

wbm_computed_entries = MaterialsProject2020Compatibility().process_entries(
    wbm_computed_entries, verbose=True, clean=True
)

n_skipped = len(df_wbm) - len(wbm_computed_entries)
assert n_skipped == 0
print(f"{n_skipped:,} ({n_skipped / len(df_wbm):.1%}) entries not processed")


# %% merge MP and WBM entries into a single PatchedPhaseDiagram
mp_wbm_ppd = PatchedPhaseDiagram(
    wbm_computed_entries + mp_computed_entries, verbose=True
)

# save MP+WBM PPD to disk (was not run)
with gzip.open(f"{module_dir}/{today}-ppd-mp-wbm.pkl.gz", "wb") as zip_file:
    pickle.dump(mp_wbm_ppd, zip_file)


# %% compute terminal reference entries across all MP (can be used to compute MP
# compatible formation energies quickly)
elemental_ref_entries = get_elemental_ref_entries(mp_computed_entries)

# save MP phase diagram elemental_ref_entries as JSON
mp_pd_elem_ref_entries_path = f"{MP_DIR}/{today}-mp-elemental-reference-entries.json.gz"
with gzip.open(mp_pd_elem_ref_entries_path, "wt") as file:
    json.dump(elemental_ref_entries, file, default=lambda x: x.as_dict())


df_mp = pd.read_csv(DATA_FILES.mp_energies, na_filter=False).set_index(Key.mat_id)


# %%
# df_mp["our_mp_e_form"] = [
#     get_e_form_per_atom(mp_computed_entries[mp_id]) for mp_id in df_mp.index
# ] #!  TypeError: list indices must be integers or slices, not str
df_mp["our_mp_e_form"] = [
    get_e_form_per_atom(mp_computed_entries[i]) for i, mp_id in enumerate(df_mp.index)
]

# make sure get_form_energy_per_atom() reproduces MP formation energies
ax = pymatviz.density_scatter(df_mp[Key.form_energy], df_mp["our_mp_e_form"])
ax.set(
    title="MP Formation Energy Comparison",
    xlabel="MP Formation Energy (eV/atom)",
    ylabel="Our Formation Energy (eV/atom)",
)
ax.figure.savefig(f"{ROOT}/tmp/{today}-our-vs-mp-formation-energies.webp", dpi=300)

# %%

[CompRhys]

I think looking at your scatter plot could you try divide all the energies on the y-axis again by the number of atoms per entry? I think that's the most likely simple bug we may be seeing and so getting confirmation would help in resolving your issue

[RyotaroOKabe]

Thanks for your reply. Maybe one possible reason for the error is the mismatch of df_mp and mp_computed_entries. I saw that len(df_mp), and len(mp_computed_entries) are (154718, 341401), respectively. I loaded the following files to get df_mp

(mp_computed_entries) 2024-03-05-mp-computed-structure-entries.json.gz
(df_mp) 2023-01-10-mp-energies.csv.gz

I am wondering if this mismatch of the data sizes might cause error.
Also, I realized my code has a problem about associating the data in mp_computed_entries and df_mp. In my previous code, I presented the following lines to get mp_computed_entries and called each ComputedEntry of mp_computed_entries for e_form calculation. However, it seems to be wrong, since the data size and the order of materials presented in mp_computed_entries and df_mp are different.

df = pd.read_json(data_path) 
mp_computed_entries = [ComputedEntry.from_dict(dct) for dct in tqdm(df.entry)]
...
df_mp["our_mp_e_form"] = [
    get_e_form_per_atom(mp_computed_entries[i]) for i, mp_id in enumerate(df_mp.index)
]

I then tried modifying the error by defining a dictionary of ComputedEntry, and call ComputedEntry by the Mp ID.

df = pd.read_json(data_path)
df = df.set_index('index')
df.index = df.index.map(lambda x: '-'.join(x.split('-')[:2]))   #!correct the mpid

mp_computed_entries_dict = {mpid: ComputedEntry.from_dict(dct) for (mpid, dct) in tqdm(zip(df.index, df.entry))}
df_mp["our_mp_e_form"] = [
    get_e_form_per_atom(mp_computed_entries_dict[mp_id]) for i, mp_id in enumerate(df_mp.index)
]

However, this still shows key error, like KeyError: 'mp-1179655'. I would be happy if you have suggestion about the debugging.

[CompRhys]

As the database evolves over time we can't guarentee that the methods originally used will work as the API/data evolves, a missing mpid is indicative of a change in dataset. We froze all our data sets on figshare if you try with those downloads we're still confident that things will work: https://matbench-discovery.materialsproject.org/contribute#--direct-download. If you wish to download "fresh" data please be aware that it may use pseudo potentials inconsistent with WBM and that we define the valid submissions as only those trained on the v2022.10.28 MP release or earlier releases.

We still are happy to upload models trained on other data but they will not be shown by default.

[RyotaroOKabe]

Thank you. Now I am still confused by calculation of Ehull with PatchedPhaseDiagram.

Below I show the code I want to run.

# load PatchedPhaseDiagram downloaded from figshare
ppd_file = "2023-02-07-ppd-mp.pkl.gz"
with gzip.open(ppd_file, "rb") as zip_file:
     ppd_mp: PatchedPhaseDiagram = pkl.load(zip_file)

file_name = "mpdata"    # load a pd.dataframe which have mpid, structure, ehull, formation energy
df_gt = pd.read_pickle(join(file_location, file_name + '.pkl'))  

ehulls_gt, ehulls_pred = [], []
for idx in range(len(df_gt)):
    row = df_gt.iloc[idx]
    mpid = row['mpid']
    ehull = row['ehull']
    structure = row['structure']
    comp = structure.composition
    f_energy = row['f_energy']
    cse = ComputedStructureEntry(structure=structure, composition=comp, energy=f_energy, entry_id=mpid)
    ehull_p = ppd_mp.get_e_above_hull(cse, allow_negative=True)
    ehulls_gt.append(ehull)
    ehulls_pred.append(ehull_p/1000)  # Is the predicted Ehull in [eV] or [meV]? 

plt.scatter(ehulls_gt, ehulls_pred)
plt.xlabel('True Ehull [eV]')
plt.ylabel('Pred Ehull [eV]')

Here are the issues I face with:
(1) ppd_mp
I wanted to load ppd_mp from "2023-02-07-ppd-mp.pkl.gz", which I downloaded directly from your figshare. I did not see any error when loading the ppd_mp with a line ppd_mp: PatchedPhaseDiagram = pkl.load(zip_file). However, once I try to use ppd_mp for Ehull calculation, I always saw the error below. Due to the error, I could not use the loaded ppd_mp, so then I tried generating the PatchedPhaseDiagram from scratch following build_phase_diagram.py. Could you tell me how I can load the ppm_mp correctly and use it for ehull calculation?

---

AttributeError Traceback (most recent call last)
File ~/anaconda3/envs/py311/lib/python3.11/site-packages/pymatgen/analysis/phase_diagram.py:770, in PhaseDiagram.get_decomp_and_e_above_hull(self, entry, allow_negative, check_stable, on_error)
769 try:
--> 770 decomp, hull_energy = self.get_decomp_and_hull_energy_per_atom(entry.composition)
771 except Exception as exc:

File ~/anaconda3/envs/py311/lib/python3.11/site-packages/pymatgen/analysis/phase_diagram.py:707, in PhaseDiagram.get_decomp_and_hull_energy_per_atom(self, comp)
706 decomp = self.get_decomposition(comp)
--> 707 return decomp, sum(e.energy_per_atom * n for e, n in decomp.items())

File ~/anaconda3/envs/py311/lib/python3.11/site-packages/pymatgen/analysis/phase_diagram.py:707, in (.0)
706 decomp = self.get_decomposition(comp)
--> 707 return decomp, sum(e.energy_per_atom * n for e, n in decomp.items())

File ~/anaconda3/envs/py311/lib/python3.11/site-packages/pymatgen/entries/init.py:86, in Entry.energy_per_atom(self)
85 """The energy per atom of the entry."""
---> 86 return self.energy / self.composition.num_atoms

File ~/anaconda3/envs/py311/lib/python3.11/site-packages/pymatgen/core/composition.py:459, in Composition.num_atoms(self)
456 """Total number of atoms in Composition. For negative amounts, sum
457 of absolute values.
458 """
--> 459 return self._n_atoms

AttributeError: 'Composition' object has no attribute '_n_atoms'

(2) Ehull calculation with ppd_mp.
As mentioned above I couldn't load ppd_mp from the saved data in your figshare, so I built from scratch. However, the predicted ehull values are so weird, even using correct structure and formation as the input of ComputedStructureEntry. Here I attach the correlation plot which i got from the code above. The output ehull values from ppd_mp.get_e_above_hull(cse, allow_negative=True) are too small if that is meV unit, but it is too high if it is eV unit. Please let me know if I am making any mistake.