New benchmark: Elemental Transition Metals Vacancy Formation Energies

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

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

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

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)

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])