NEB surface reaction

docs/source/user_guide/benchmarks/nebs.rst

@@ -46,3 +46,42 @@ Reference data:
 
 * Manually taken from https://doi.org/10.1149/1.1633511.
 * Meta-GGA (Perdew-Wang) exchange correlation functional
+
+
+Surface reaction
+============
+
+Summary
+-------
+
+Performance in running NEB for three surface reactions (adsorption/desorption, dissociation, transfer) from OC20NEB dataset.
+
+Metrics
+-------
+
+1. Activation barrier error
+
+Initial and final geometries are from OC20NEB dataset and relaxed with each model with several layers of slab fixed. And interpolation is generated with 10 images including initial and final images. Following NEB setting from original paper, initially it runs without climbing image until fmax=0.45 eV/A or max steps 200 and converts to climbing image mode with fmax=0.05 eV/A or max steps 300. Barrier is measured between the highest energy point and initial image.
+
+The benchmark includes following reactions:
+
+- desorption_ood_87_9841_0_111-1
+- dissociation_ood_268_6292_46_211-5
+- transfer_id_601_1482_1_211-5
+
+Computational cost
+------------------
+
+Slow: tests are likely to take more than 2 hours to run on single GPU.
+
+Data availability
+-----------------
+
+Input structure:
+
+* OC20NEB dataset : https://dl.fbaipublicfiles.com/opencatalystproject/data/oc20neb/oc20neb_dft_trajectories_04_23_24.tar.gz
+
+Reference data:
+
+* Manually taken from https://doi.org/10.1021/acscatal.4c04272
+* GGA RPBE exchange correlation functional

ml_peg/analysis/nebs/surface_reaction/analyse_surface_reaction.py

@@ -0,0 +1,151 @@
+"""Analyse Li diffusion benchmark."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+from ase.io import read, write
+import numpy as np
+import pytest
+
+from ml_peg.analysis.utils.decorators import build_table, plot_scatter
+from ml_peg.analysis.utils.utils import load_metrics_config
+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)
+DATA_PATH = CALCS_ROOT / "nebs" / "surface_reaction" / "data"
+CALC_PATH = CALCS_ROOT / "nebs" / "surface_reaction" / "outputs"
+OUT_PATH = APP_ROOT / "data" / "nebs" / "surface_reaction"
+
+METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml")
+DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config(
+    METRICS_CONFIG_PATH
+)
+
+REF_VALUES = {
+    "desorption_ood_87_9841_0_111-1": 2.061,
+    "dissociation_ood_268_6292_46_211-5": 1.505,
+    "transfer_id_601_1482_1_211-5": 0.868,
+}
+REACTIONS = [
+    "desorption_ood_87_9841_0_111-1",
+    "dissociation_ood_268_6292_46_211-5",
+    "transfer_id_601_1482_1_211-5",
+]
+
+
+def plot_nebs(model: str, reaction: str) -> None:
+    """
+    Plot NEB paths and save all structure files.
+
+    Parameters
+    ----------
+    model
+        Name of MLIP.
+    reaction
+        Reaction id for NEB.
+    """
+
+    @plot_scatter(
+        filename=OUT_PATH / f"figure_{model}_neb_{reaction}.json",
+        title=f"NEB path {reaction}",
+        x_label="Image",
+        y_label="Energy / eV",
+        show_line=True,
+    )
+    def plot_neb() -> dict[str, tuple[list[float], list[float]]]:
+        """
+        Plot a NEB and save the structure file.
+
+        Returns
+        -------
+        dict[str, tuple[list[float], list[float]]]
+            Dictionary of tuples of image/energy for each model.
+        """
+        results = {}
+        structs = read(
+            #            CALC_PATH / f"surface_reaction_{reaction}-{model}.xyz",
+            CALC_PATH / f"{reaction}_{model}.xyz",
+            index=":",
+        )
+        results[model] = [
+            list(range(len(structs))),
+            [struct.info["mlip_energy"] for struct in structs],
+        ]
+        structs_dir = OUT_PATH / model
+        structs_dir.mkdir(parents=True, exist_ok=True)
+        write(structs_dir / f"{model}-{reaction}.xyz", structs)
+
+        return results
+
+    plot_neb()
+
+
+@pytest.fixture
+def barrier_error() -> dict[str, dict[str, float]]:
+    """
+    Get error in energy barrier for all reactions.
+
+    Returns
+    -------
+    dict[str, float]
+        Dictionary of predicted barrier errors for all models.
+    """
+    #    OUT_PATH.mkdir(parents=True, exist_ok=True)
+    results = {}
+    for model_name in MODELS:
+        reaction_dict = {}
+        for reaction in REACTIONS:
+            plot_nebs(model_name, reaction)
+            structs = read(CALC_PATH / f"{reaction}_{model_name}.xyz", ":")
+            energies = [struct.info["mlip_energy"] for struct in structs]
+            pred_forward_barrier = np.max(energies) - energies[0]
+            reaction_dict[reaction] = np.abs(
+                pred_forward_barrier - REF_VALUES[reaction]
+            )
+        results[model_name] = reaction_dict
+
+    return results
+
+
+@pytest.fixture
+@build_table(
+    filename=OUT_PATH / "surface_reaction_metrics_table.json",
+    metric_tooltips=DEFAULT_TOOLTIPS,
+    thresholds=DEFAULT_THRESHOLDS,
+)
+def metrics(barrier_error: dict[str, dict[str, float]]) -> dict[str, dict]:
+    """
+    Get all surface reactions metrics.
+
+    Parameters
+    ----------
+    barrier_error
+        Activation barriers for all reactions and all models.
+
+    Returns
+    -------
+    dict[str, dict]
+        Metric names and values for all models.
+    """
+    metrics_dict = {}
+    for reaction in REACTIONS:
+        metrics_dict[f"{reaction} barrier error"] = {
+            model: barrier_error[model][reaction] for model in MODELS
+        }
+    return metrics_dict
+
+
+def test_surface_reaction(metrics: dict[str, dict]) -> None:
+    """
+    Run surface reaction test.
+
+    Parameters
+    ----------
+    metrics
+        All surface reaction metrics.
+    """
+    return

ml_peg/analysis/nebs/surface_reaction/metrics.yml

@@ -0,0 +1,19 @@
+metrics:
+  desorption_ood_87_9841_0_111-1 barrier error:
+    good: 0.1
+    bad: 0.5
+    unit: eV
+    tooltip: "desorption_ood_87_9841_0_111-1 barrier error"
+    level_of_theory: RPBE
+  dissociation_ood_268_6292_46_211-5 barrier error:
+    good: 0.1
+    bad: 0.5
+    unit: eV
+    tooltip: "dissociation_ood_268_6292_46_211-5 barrier error"
+    level_of_theory: RPBE
+  transfer_id_601_1482_1_211-5 barrier error:
+    good: 0.1
+    bad: 0.5
+    unit: eV
+    tooltip: "transfer_id_601_1482_1_211-5 barrier error"
+    level_of_theory: RPBE

ml_peg/app/nebs/surface_reaction/app_surface_reaction.py

@@ -0,0 +1,106 @@
+"""Run Li diffusion 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_cell,
+    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 = "Surface reaction"
+DOCS_URL = (
+    "https://ddmms.github.io/ml-peg/user_guide/benchmarks/nebs.html#surface-reactionn"
+)
+DATA_PATH = APP_ROOT / "data" / "nebs" / "surface_reaction"
+
+REACTIONS = [
+    "desorption_ood_87_9841_0_111-1",
+    "dissociation_ood_268_6292_46_211-5",
+    "transfer_id_601_1482_1_211-5",
+]
+
+
+class SurfaceReactionApp(BaseApp):
+    """Surface reaction benchmark app layout and callbacks."""
+
+    def register_callbacks(self) -> None:
+        """Register callbacks to app."""
+        scatter_plots = {
+            model: {
+                f"{reaction} barrier error": read_plot(
+                    DATA_PATH / f"figure_{model}_neb_{reaction}.json",
+                    id=f"{BENCHMARK_NAME}-{model}-figure-{reaction}",
+                )
+                for reaction in REACTIONS
+            }
+            for model in MODELS
+        }
+
+        # Assets dir will be parent directory
+        assets_dir = "assets/nebs/surface_reaction"
+        structs = {
+            model: {
+                f"{rxn} barrier error": f"{assets_dir}/{model}/{model}-{rxn}.xyz"
+                for rxn in REACTIONS
+            }
+            for model in MODELS
+        }
+
+        plot_from_table_cell(
+            table_id=self.table_id,
+            plot_id=f"{BENCHMARK_NAME}-figure-placeholder",
+            cell_to_plot=scatter_plots,
+        )
+
+        for model in MODELS:
+            for reaction in REACTIONS:
+                struct_from_scatter(
+                    scatter_id=f"{BENCHMARK_NAME}-{model}-figure-{reaction}",
+                    struct_id=f"{BENCHMARK_NAME}-struct-placeholder",
+                    structs=structs[model][f"{reaction} barrier error"],
+                    mode="traj",
+                )
+
+
+def get_app() -> SurfaceReactionApp:
+    """
+    Get Li diffusion benchmark app layout and callback registration.
+
+    Returns
+    -------
+    SurfaceReactionApp
+        Benchmark layout and callback registration.
+    """
+    return SurfaceReactionApp(
+        name=BENCHMARK_NAME,
+        description=("Performance in predicting energy barriers for Surface reaction."),
+        docs_url=DOCS_URL,
+        table_path=DATA_PATH / "surface_reaction_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)
+
+    # Construct layout and register callbacks
+    surface_reaction_app = get_app()
+    full_app.layout = surface_reaction_app.layout
+    surface_reaction_app.register_callbacks()
+
+    # Run app
+    full_app.run(port=8051, debug=True)