Metal surface reconstuctions benchmark

docs/source/user_guide/benchmarks/surfaces.rst

@@ -145,3 +145,45 @@ Reference data:
   * S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. Chevrier, K. A. Persson, G. Ceder, "Python Materials Genomics (pymatgen): A Robust, Open-Source Python Library for Materials Analysis," Comput. Mater. Sci., 2013, 68, 314–319. https://doi.org/10.1016/j.commatsci.2012.10.028
 
 * Tran et al. relaxed the slabs using spin-polarized PBE calculations performed in VASP, with a cutoff energy of 400 eV.
+
+
+Metal Surface Reconstructions
+================================
+
+Summary
+-------
+
+Performance in the predicting surface energy and representing the surface reconstruction. 
+
+Metrics
+-------
+
+For each slab, a geometry optimization (F_max = 0.05 eV/Å) is performed, where the bottom most layer is fixed. After that three metrics are evaluated:
+
+* Surface energy error (the surface energy is calculated with regard to the bulk and gas phase references)
+
+* Displacement with regard to the reference configuration
+
+* Energetic ranking (using the surface energy)
+
+Computational cost
+------------------
+
+Very low: tests are likely to take half an hour.
+
+Data availability
+-----------------
+
+Input data:
+
+* Stable and metastable Cu(111)O surface reconstructions.
+
+  * Zhu, B., Huang, Y., Lv, J., Huang, W., Lian, Z., Ouyang, R.,  Yang, F. "Dynamic Evolution and Transformation of Copper Oxides on Cu (111)." The Journal of Physical Chemistry C, 2025, 129(29), 13497-13504. https://doi.org/10.1021/acs.jpcc.5c03416
+
+* Stable and metastable Pd/Au reconstructions with C, O and N atoms.
+
+  * Vinogradova, O. V., Reuter, K., Bukas, V. J. "Trends of Pd3Au (111) alloy surface segregation in oxygen, carbon, and nitrogen environments." The Journal of Physical Chemistry C, 2023, 127(45), 22060-22066. https://doi.org/10.1021/acs.jpcc.3c05640
+
+Reference data:
+
+* PBE conjugate gradient geomerty optimizations are performed via VASP, with a cutoff energy of 500 eV.

ml_peg/analysis/surfaces/metal_surface_reconstructions/analyse_metal_surface_reconstructions.py

@@ -0,0 +1,311 @@
+
+from __future__ import annotations
+
+from pathlib import Path
+
+from ase import units
+from ase.io import read, write
+import numpy as np
+import pytest
+
+from ml_peg.analysis.utils.decorators import build_table, plot_parity
+from ml_peg.analysis.utils.utils import build_d3_name_map, 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)
+print(MODELS)
+#D3_MODEL_NAMES = build_d3_name_map(MODELS)
+CALC_PATH = CALCS_ROOT / "surfaces" / "metal_surfaces" / "outputs"
+OUT_PATH = APP_ROOT / "data" / "surfaces" / "metal_surfaces"
+
+METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml")
+DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config(
+    METRICS_CONFIG_PATH
+)
+
+# Unit conversion
+EV_TO_KJ_PER_MOL = units.mol / units.kJ
+
+
+def get_system_names() -> list[str]:
+    """
+    Get list of metal surface 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
+        print(model_dir)
+        if model_dir.exists():
+            print('yes')
+            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"])
+
+    return system_names
+
+
+
+@pytest.fixture
+@plot_parity(
+    filename=OUT_PATH / "slab_energies.json",
+    title="Metal Slab Energies",
+    x_label="Predicted Surface Free Energy / meV/Ų",
+    y_label="Reference EneSurface Free Energy / meV/Ų",
+    hoverdata={
+        "System": get_system_names(),
+    },
+)
+def slab_energies() -> dict[str, list]:
+    """
+    Get energies for all slabs systems.
+
+    Returns
+    -------
+    dict[str, list]
+        Dictionary of reference and predicted lattice energies.
+    """
+    results = {"ref": []} | {mlip: [] for mlip in MODELS}
+    ref_stored = False
+
+    ref_mu = {}
+
+    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
+
+        model_mu = {}
+        for xyz_file in xyz_files:
+            name = xyz_file.name
+            if name.startswith('bulk') or name.startswith('gas_phase'):
+                structs = read(xyz_file, index=":")[0]
+                model_mu[structs.get_chemical_symbols()[0]] = structs.get_potential_energy()/len(structs)
+                if not ref_stored:
+                    ref_mu[structs.get_chemical_symbols()[0]] = structs.info['DFT_energy']/len(structs)
+
+
+
+        for xyz_file in xyz_files:
+            name = xyz_file.name
+            if not (name.startswith('bulk') or name.startswith('gas_phase')):
+                structs = read(xyz_file, index=":")[0]
+                system = structs.info["system"]
+                symbols = structs.get_chemical_symbols()
+                cell = structs.cell
+                A = np.linalg.norm(np.cross(cell[0], cell[1]))
+
+                results[model_name].append((structs.get_potential_energy()-np.sum([model_mu[s] for s in symbols]))*1000/A)
+
+
+
+                # 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.info["DFT_energy"]-np.sum([ref_mu[s] for s in symbols]))*1000/A)
+
+        ref_stored = True
+
+    return results
+
+
+
+
+
+
+
+@pytest.fixture
+def slab_positions() -> dict[str, list]:
+    """
+    Get energies for all slabs systems.
+
+    Returns
+    -------
+    dict[str, list]
+        Dictionary of reference and predicted lattice 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:
+            name = xyz_file.name
+            if not (name.startswith('bulk') or name.startswith('gas_phase')):
+                structs = read(xyz_file, index=":")[0]
+                z_min = np.min(structs.positions[:,2])
+                moving = structs.positions[:,2]>z_min+0.1
+                results[model_name].append(structs.positions[moving])
+
+
+
+                # Store reference energies (only once)
+                if not ref_stored:
+                    results["ref"].append(structs.arrays["DFT_positions"][moving])
+
+        ref_stored = True
+
+    return results
+
+
+@pytest.fixture
+def ranking_error(slab_energies) -> dict[str, float]:
+    """
+    Get ranking error across all triplets.
+
+    Parameters
+    ----------
+    relative_energies
+        Dictionary of reference and predicted relative energies.
+
+    Returns
+    -------
+    dict[str, float]
+        Dictionary of predicted ranking errors for all models.
+    """
+    print(slab_energies.keys())
+    results = {}
+    ref_min = []
+    ref_max = []
+    for i in range(len(slab_energies["ref"]) // 3):
+        ref_energies = slab_energies["ref"][3 * i : 3 * i + 3]
+        ref_min.append(np.argmin(ref_energies))
+        ref_max.append(np.argmax(ref_energies))
+
+    for model_name in MODELS:
+        if slab_energies[model_name]:
+            pred_min = []
+            pred_max = []
+            for i in range(len(slab_energies[model_name]) // 3):
+                pred_energies = slab_energies[model_name][3 * i : 3 * i + 3]
+                pred_min.append(np.argmin(pred_energies))
+                pred_max.append(np.argmax(pred_energies))
+
+            results[model_name] = (
+                1
+                - 0.5 * np.mean(np.array(ref_min) == np.array(pred_min))
+                - 0.5 * np.mean(np.array(ref_max) == np.array(pred_max))
+            )
+        else:
+            results[model_name] = None
+
+    return results
+
+
+
+
+@pytest.fixture
+def metal_surfaces_errors(slab_energies) -> dict[str, float]:
+    """
+    Get mean absolute error for lattice energies.
+
+    Parameters
+    ----------
+    lattice_energies
+        Dictionary of reference and predicted lattice energies.
+
+    Returns
+    -------
+    dict[str, float]
+        Dictionary of predicted lattice energy errors for all models.
+    """
+    results = {}
+    for model_name in MODELS:
+        if slab_energies[model_name]:
+            results[model_name] = mae(
+                slab_energies["ref"], slab_energies[model_name]
+            )
+        else:
+            results[model_name] = None
+    return results
+
+
+
+@pytest.fixture
+def metal_position_errors(slab_positions) -> dict[str, float]:
+    """
+    Get mean absolute error for lattice energies.
+
+    Parameters
+    ----------
+    lattice_energies
+        Dictionary of reference and predicted lattice energies.
+
+    Returns
+    -------
+    dict[str, float]
+        Dictionary of predicted lattice energy errors for all models.
+    """
+    results = {}
+    for model_name in MODELS:
+        if slab_positions[model_name]:
+            results[model_name] = (
+                np.mean(np.linalg.norm(np.concat(slab_positions["ref"])-np.concat(slab_positions[model_name]),axis=1))
+            )
+        else:
+            results[model_name] = None
+    return results
+
+@pytest.fixture
+@build_table(
+    filename=OUT_PATH / "metal_surfaces_metrics_table.json",
+    metric_tooltips=DEFAULT_TOOLTIPS,
+    thresholds=DEFAULT_THRESHOLDS,
+    mlip_name_map=None,
+)
+def metrics(metal_surfaces_errors: dict[str, float], metal_position_errors: dict[str, float], ranking_error:  dict[str, float]) -> dict[str, dict]:
+    """
+    Get all metal surface metrics.
+
+    Parameters
+    ----------
+    metal_surfaces_errors
+        Mean absolute errors for all systems.
+
+    Returns
+    -------
+    dict[str, dict]
+        Metric names and values for all models.
+    """
+    return {
+        "MAE": metal_surfaces_errors,
+        "Displacement": metal_position_errors,
+        "Ranking Error": ranking_error,
+    }
+
+def test_metal_surfaces(metrics: dict[str, dict]) -> None:
+    """
+    Run metal surface test.
+
+    Parameters
+    ----------
+    metrics
+        All metal surface metrics.
+    """
+    return

ml_peg/analysis/surfaces/metal_surface_reconstructions/metrics.yml

@@ -0,0 +1,19 @@
+metrics:
+    MAE:
+        good: 0.0
+        bad: 20.0
+        unit: meV/Ų
+        tooltip: "Mean Absolute Error for all systems"
+        level_of_theory: PBE
+    Displacement:
+        good: 0.0
+        bad: 0.1
+        unit: Å
+        tooltip: "Mean Atom displacement"
+        level_of_theory: PBE
+    Ranking Error:
+        good: 0.0
+        bad: 1.0
+        unit: null
+        tooltip: "Error in ranking stability across neighboring triplets"
+        level_of_theory: PBE