From 89f3f9aeb91feebf2354d4092c52fba5e942389c Mon Sep 17 00:00:00 2001 From: alxdr-peuch Date: Mon, 23 Feb 2026 14:23:11 +0000 Subject: [PATCH] Add new benchmark for Elemental Transition Metals Vacancy Formation Energies --- .../user_guide/benchmarks/bulk_crystal.rst | 40 +++++ .../analyse_elemental_tm_vacancies.py | 169 ++++++++++++++++++ .../elemental_tm_vacancies/metrics.yml | 7 + .../app_elemental_tm_vacancies.py | 87 +++++++++ .../calc_elemental_tm_vacancies.py | 77 ++++++++ 5 files changed, 380 insertions(+) create mode 100644 ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/analyse_elemental_tm_vacancies.py create mode 100644 ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/metrics.yml create mode 100644 ml_peg/app/bulk_crystal/elemental_tm_vacancies/app_elemental_tm_vacancies.py create mode 100644 ml_peg/calcs/bulk_crystal/elemental_tm_vacancies/calc_elemental_tm_vacancies.py diff --git a/docs/source/user_guide/benchmarks/bulk_crystal.rst b/docs/source/user_guide/benchmarks/bulk_crystal.rst index 575bbcc11..55c825cd4 100644 --- a/docs/source/user_guide/benchmarks/bulk_crystal.rst +++ b/docs/source/user_guide/benchmarks/bulk_crystal.rst @@ -115,3 +115,43 @@ Reference data: * Same as input data * PBE + + +Elemental TM Vacancy Formation Energies +======================================= + +Summary +------- + +Performance in predicting vacancy formation energies for 42 fcc or hcp elemental transition metal structures. + +Metrics +------- + +1. MAE + +Mean absolute error (MAE) between predicted and reference (PBE) vacancy formation energies values. + +For each elemental transition metal structure, the vacancy formation enthalpy is determined by evaluating the total energy of the cell containing the monovacancy and evaluating the total energy of the undefected cell and comparing these. Note: the undefected cell energy is scaled by a fraction in order to conserve the appropriate number of atoms. + +From the reference paper providing the structures, the elemental transition metals are a subset. Note: + * For magnetic structures, the reference spin-polarized calculations only is used. + * Unstable element-structure combinations are not included. + +Computational cost +------------------ + +Low: tests are likely to take a couple of minutes to run on CPU. + + +Data availability +----------------- + +Input structures: + +* T. Angsten et al. Elemental vacancy diffusion database from high-throughput first-principles calculations for fcc and hcp structures. New J. Phys. 2024 16 015018 + +Reference data: + +* Same as input data +* PBE diff --git a/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/analyse_elemental_tm_vacancies.py b/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/analyse_elemental_tm_vacancies.py new file mode 100644 index 000000000..6916ba37b --- /dev/null +++ b/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/analyse_elemental_tm_vacancies.py @@ -0,0 +1,169 @@ +"""Analyse Elemental TM Vacancy Formation Energies benchmark.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read, write +import pytest + +from ml_peg.analysis.utils.decorators import build_table, plot_parity +from ml_peg.analysis.utils.utils import load_metrics_config, mae +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +MODELS = get_model_names(current_models) +CALC_PATH = CALCS_ROOT / "bulk_crystal" / "elemental_tm_vacancies" / "outputs" +OUT_PATH = APP_ROOT / "data" / "bulk_crystal" / "elemental_tm_vacancies" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( + METRICS_CONFIG_PATH +) + + +def get_system_names() -> list[str]: + """ + Get list of all system names. + + Returns + ------- + list[str] + List of system names from structure files. + """ + system_names = [] + for model_name in MODELS: + model_dir = CALC_PATH / model_name + if model_dir.exists(): + xyz_files = sorted(model_dir.glob("*.xyz")) + if xyz_files: + for xyz_file in xyz_files: + atoms = read(xyz_file) + system_names.append(atoms.info["system"]) + break + return system_names + + +@pytest.fixture +@plot_parity( + filename=OUT_PATH / "figure_vacancy_formation_energies.json", + title="Elemental TM Vacancy Formation Energies", + x_label="Predicted Vacancy Formation Energy / eV", + y_label="Reference Vacancy Formation Energy / eV", + hoverdata={ + "System": get_system_names(), + }, +) +def vacancy_formation_energies() -> dict[str, list]: + """ + Get vacancy formation energies for all systems. + + Returns + ------- + dict[str, list] + Dictionary of reference and predicted vacancy formation energies. + """ + results = {"ref": []} | {mlip: [] for mlip in MODELS} + ref_stored = False + + for model_name in MODELS: + model_dir = CALC_PATH / model_name + + if not model_dir.exists(): + continue + + xyz_files = sorted(model_dir.glob("*.xyz")) + if not xyz_files: + continue + + for xyz_file in xyz_files: + structs = read(xyz_file, index=":") + + bulk_energy = structs[0].get_potential_energy() + num_atoms = len(structs[0]) + system = structs[0].info["system"] + vacancy_energy = structs[1].get_potential_energy() + + vacancy_formation_energy = ( + vacancy_energy - ((num_atoms - 1) / num_atoms) * bulk_energy + ) + results[model_name].append(vacancy_formation_energy) + + # Copy individual structure files to app data directory + structs_dir = OUT_PATH / model_name + structs_dir.mkdir(parents=True, exist_ok=True) + write(structs_dir / f"{system}.xyz", structs) + + # Store reference energies (only once) + if not ref_stored: + results["ref"].append(structs[0].info["ref"]) + + ref_stored = True + + return results + + +@pytest.fixture +def vacancy_formation_errors(vacancy_formation_energies) -> dict[str, float]: + """ + Get mean absolute error for vacancy formation energies. + + Parameters + ---------- + vacancy_formation_energies + Dictionary of reference and predicted vacancy formation energies. + + Returns + ------- + dict[str, float] + Dictionary of predicted vacancy formation energy errors for all models. + """ + results = {} + for model_name in MODELS: + if vacancy_formation_energies[model_name]: + results[model_name] = mae( + vacancy_formation_energies["ref"], + vacancy_formation_energies[model_name], + ) + else: + results[model_name] = None + return results + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "vacancy_formation_energies_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, +) +def metrics(vacancy_formation_errors: dict[str, float]) -> dict[str, dict]: + """ + Get all Vacancy Formation Energies metrics. + + Parameters + ---------- + vacancy_formation_errors + Mean absolute errors for all systems. + + Returns + ------- + dict[str, dict] + Metric names and values for all models. + """ + return { + "MAE": vacancy_formation_errors, + } + + +def test_vacancy_formation_energies(metrics: dict[str, dict]) -> None: + """ + Run test to assess a model on elemental TM vac. form. energies. + + Parameters + ---------- + metrics + All elemental TM vacancy formation energies metrics. + """ + return diff --git a/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/metrics.yml b/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/metrics.yml new file mode 100644 index 000000000..fe2f1432d --- /dev/null +++ b/ml_peg/analysis/bulk_crystal/elemental_tm_vacancies/metrics.yml @@ -0,0 +1,7 @@ +metrics: + MAE: + good: 0.05 + bad: 0.5 + unit: eV + tooltip: "Mean Absolute Error for all systems" + level_of_theory: PBE diff --git a/ml_peg/app/bulk_crystal/elemental_tm_vacancies/app_elemental_tm_vacancies.py b/ml_peg/app/bulk_crystal/elemental_tm_vacancies/app_elemental_tm_vacancies.py new file mode 100644 index 000000000..7f0a84a1b --- /dev/null +++ b/ml_peg/app/bulk_crystal/elemental_tm_vacancies/app_elemental_tm_vacancies.py @@ -0,0 +1,87 @@ +"""Run Elemental TM Vacancy app.""" + +from __future__ import annotations + +from dash import Dash +from dash.html import Div + +from ml_peg.app import APP_ROOT +from ml_peg.app.base_app import BaseApp +from ml_peg.app.utils.build_callbacks import ( + plot_from_table_column, + struct_from_scatter, +) +from ml_peg.app.utils.load import read_plot +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +# Get all models +MODELS = get_model_names(current_models) +BENCHMARK_NAME = "Elemental Transition Metal Vacancy Formation Energies" +DOCS_URL = "https://ddmms.github.io/ml-peg/user_guide/benchmarks/bulk_crystal.html#elemental-tm-vacancy-formation-energies" +DATA_PATH = APP_ROOT / "data" / "bulk_crystal" / "elemental_tm_vacancies" + + +class ElementalTMVacanciesApp(BaseApp): + """Elemental TM Vacancies benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + scatter = read_plot( + DATA_PATH / "figure_vacancy_formation_energies.json", + id=f"{BENCHMARK_NAME}-figure", + ) + + # Assets dir will be parent directory - individual files for each system + structs_dir = DATA_PATH / MODELS[0] + structs = [ + f"assets/bulk_crystal/elemental_tm_vacancies/{MODELS[0]}/{struct_file.stem}.xyz" + for struct_file in sorted(structs_dir.glob("*.xyz")) + ] + + plot_from_table_column( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + column_to_plot={"MAE": scatter}, + ) + + struct_from_scatter( + scatter_id=f"{BENCHMARK_NAME}-figure", + struct_id=f"{BENCHMARK_NAME}-struct-placeholder", + structs=structs, + mode="traj", + ) + + +def get_app() -> ElementalTMVacanciesApp: + """ + Get Elemental TM Vacancies benchmark app layout and callback registration. + + Returns + ------- + ElementalTMVacanciesApp + Benchmark layout and callback registration. + """ + return ElementalTMVacanciesApp( + name=BENCHMARK_NAME, + description="Vacancy formation energies for 42 elemental TM structures.", + docs_url=DOCS_URL, + table_path=DATA_PATH / "vacancy_formation_energies_metrics_table.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + Div(id=f"{BENCHMARK_NAME}-struct-placeholder"), + ], + ) + + +if __name__ == "__main__": + # Create Dash app + full_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent) + + # Construct layout and register callbacks + elemental_tm_vacancies_app = get_app() + full_app.layout = elemental_tm_vacancies_app.layout + elemental_tm_vacancies_app.register_callbacks() + + # Run app + full_app.run(port=8055, debug=True) diff --git a/ml_peg/calcs/bulk_crystal/elemental_tm_vacancies/calc_elemental_tm_vacancies.py b/ml_peg/calcs/bulk_crystal/elemental_tm_vacancies/calc_elemental_tm_vacancies.py new file mode 100644 index 000000000..dfae0f168 --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/elemental_tm_vacancies/calc_elemental_tm_vacancies.py @@ -0,0 +1,77 @@ +"""Run calculations for Elemental Transition Metal Vacancy Formation tests.""" + +from __future__ import annotations + +from copy import copy +from pathlib import Path +from typing import Any + +from ase.io import read, write +import numpy as np +import pytest + +from ml_peg.calcs.utils.utils import download_s3_data +from ml_peg.models.get_models import load_models +from ml_peg.models.models import current_models + +MODELS = load_models(current_models) + +DATA_PATH = Path(__file__).parent / "data" +OUT_PATH = Path(__file__).parent / "outputs" + + +@pytest.mark.parametrize("mlip", MODELS.items()) +def test_vacancy_formation_energy(mlip: tuple[str, Any]) -> None: + """ + Run vacancy formation energy test. + + Parameters + ---------- + mlip + Name of model use and model to get calculator. + """ + model_name, model = mlip + calc = model.get_calculator() + + # Download dataset + elemental_tm_vacancies_dir = ( + download_s3_data( + key="inputs/bulk_crystal/elemental_tm_vacancies/Elemental_TM_Vacancy.zip", + filename="Elemental_TM_Vacancy.zip", + ) + / "Elemental_TM_Vacancy" + ) + + with open(elemental_tm_vacancies_dir / "list") as f: + systems = f.read().splitlines() + + for system in systems: + bulk_path = elemental_tm_vacancies_dir / system / "CONTCAR_bulk" + vacancy_path = elemental_tm_vacancies_dir / system / "CONTCAR_vacancy" + ref_path = elemental_tm_vacancies_dir / system / "vacancy_formation_energy_PBE" + + bulk = read(bulk_path, index=0, format="vasp") + # Set default charge and spin + bulk.info.setdefault("charge", 0) + bulk.info.setdefault("spin", 1) + bulk.calc = calc + bulk.get_potential_energy() + + vacancy = read(vacancy_path, index=0, format="vasp") + # Set default charge and spin + vacancy.info.setdefault("charge", 0) + vacancy.info.setdefault("spin", 1) + vacancy.calc = copy(calc) + vacancy.get_potential_energy() + + ref = np.loadtxt(ref_path) + + bulk.info["ref"] = ref + bulk.info["system"] = system + vacancy.info["ref"] = ref + vacancy.info["system"] = system + + # Write output structures + write_dir = OUT_PATH / model_name + write_dir.mkdir(parents=True, exist_ok=True) + write(write_dir / f"{system}.xyz", [bulk, vacancy])