From fd7957608a8ec94d8b6349fa9415245a82ab4e63 Mon Sep 17 00:00:00 2001
From: Anyang Peng <137014849+anyangml@users.noreply.github.com>
Date: Thu, 27 Aug 2026 15:05:43 +0800
Subject: [PATCH 1/2] fix: normalize roughness

---
 lambench/metrics/utils.py                     | 14 ++-
 .../tasks/calculator/diatomics/diatomics.py   | 78 +++++++++++----
 tests/tasks/calculator/test_diatomics.py      | 97 ++++++++++++-------
 3 files changed, 129 insertions(+), 60 deletions(-)

diff --git a/lambench/metrics/utils.py b/lambench/metrics/utils.py
index 5bc706e..58d00ec 100644
--- a/lambench/metrics/utils.py
+++ b/lambench/metrics/utils.py
@@ -168,28 +168,32 @@ def _diatomics_molecule_names() -> tuple[str, ...]:
         return tuple(entry["name"] for entry in json.load(fh))
 
 
+def _empty_diatomics_agg() -> dict[str, float | None]:
+    return {"avg_roughness": None}
+
+
 def aggregated_diatomics_results(results: dict[str, dict]) -> dict[str, float]:
     """
     Aggregate per-molecule diatomics results.
 
-    avg_roughness: arithmetic mean of curvature RMSE (eV/Ų), used on the leaderboard.
+    avg_roughness: mean slope MAE relative to a constant dummy, capped at 1.
     Requires a finite roughness for every molecule in diatomics.json; otherwise None.
     """
     if not results:
-        return {"avg_roughness": None}
+        return _empty_diatomics_agg()
 
     names = _diatomics_molecule_names()
     if not names:
-        return {"avg_roughness": None}
+        return _empty_diatomics_agg()
 
     roughness_values = []
     for name in names:
         mol_results = results.get(name)
         if mol_results is None:
-            return {"avg_roughness": None}
+            return _empty_diatomics_agg()
         roughness = mol_results.get("roughness")
         if roughness is None or not np.isfinite(roughness):
-            return {"avg_roughness": None}
+            return _empty_diatomics_agg()
         roughness_values.append(roughness)
 
     return {"avg_roughness": float(np.mean(roughness_values))}
diff --git a/lambench/tasks/calculator/diatomics/diatomics.py b/lambench/tasks/calculator/diatomics/diatomics.py
index c092e94..542030e 100644
--- a/lambench/tasks/calculator/diatomics/diatomics.py
+++ b/lambench/tasks/calculator/diatomics/diatomics.py
@@ -20,15 +20,15 @@
 
 Homonuclear diatomics dissociation curves (Applicability).
 
-Per-molecule metric, then arithmetic mean over molecules:
+Same comparison as stacking_fault, on the DFT bond lengths (no PCHIP):
+shift each scan so min(E)=0, then MAE of d(E / max(E)) / dr, divided
+by a constant-energy dummy (zeros after the shift) and capped at 1:
 
-    roughness – RMSE of d²(E_model - E_DFT)/dr²  (eV/Ų)
+    0 = perfect match, 1 = dummy or worse.
 
-Leaderboard (Applicability-Roughness ↓) uses avg_roughness.
+    roughness – slope MAE relative to dummy; leaderboard
 
-The last point of every scan is dropped: some PBE dissociation tails have a
-spurious endpoint jump, and trimming all scans the same way avoids per-molecule
-special cases.
+The last point of every scan is dropped (PBE tail artifacts).
 
 Reference data: lambench/tasks/calculator/diatomics/diatomics.json
     name, method, R, E, F, S^2
@@ -73,13 +73,58 @@ def _scan_arrays(entry: dict) -> tuple[np.ndarray, np.ndarray]:
     return bond_lengths, dft_energies
 
 
-def _compute_roughness(residuals: np.ndarray, dr: float) -> float | None:
-    """RMSE of d² residual / dr².  None if too few finite points."""
-    delta2 = np.diff(residuals, n=2)
-    valid = delta2[np.isfinite(delta2)]
-    if len(valid) == 0:
+def _shift_to_min(energies: np.ndarray) -> np.ndarray:
+    return energies - np.min(energies)
+
+
+def _normalized_slopes(bond_lengths: np.ndarray, y: np.ndarray) -> np.ndarray:
+    """d(E / max(E)) / dr on the native DFT scan (Å⁻¹)."""
+    peak = float(np.max(y))
+    slopes = np.diff(y) / np.diff(bond_lengths)
+    if peak > 0:
+        return slopes / peak
+    return np.zeros_like(slopes)
+
+
+def _mae(left: np.ndarray, right: np.ndarray) -> float:
+    return float(np.mean(np.abs(left - right)))
+
+
+def _ratio_capped_at_dummy(value: float, dummy: float) -> float | None:
+    if dummy <= 0 or not np.isfinite(dummy) or not np.isfinite(value):
+        return None
+    return float(min(value / dummy, 1.0))
+
+
+def _curve_metrics(
+    bond_lengths: np.ndarray, model_energies: np.ndarray, dft_energies: np.ndarray
+) -> dict[str, float] | None:
+    if not np.all(np.isfinite(model_energies)) or not np.all(np.isfinite(dft_energies)):
+        return None
+    if bond_lengths.size < 2:
+        return None
+    if bond_lengths.size != model_energies.size or bond_lengths.size != dft_energies.size:
+        return None
+    if np.any(np.diff(bond_lengths) == 0):
         return None
-    return float(np.sqrt(np.mean((valid / dr**2) ** 2)))
+
+    y_dft = _shift_to_min(dft_energies)
+    y_model = _shift_to_min(model_energies)
+    y_dummy = np.zeros_like(y_dft)
+
+    roughness = _ratio_capped_at_dummy(
+        _mae(
+            _normalized_slopes(bond_lengths, y_model),
+            _normalized_slopes(bond_lengths, y_dft),
+        ),
+        _mae(
+            _normalized_slopes(bond_lengths, y_dummy),
+            _normalized_slopes(bond_lengths, y_dft),
+        ),
+    )
+    if roughness is None:
+        return None
+    return {"roughness": roughness}
 
 
 def _predict_energies(
@@ -111,7 +156,7 @@ def _predict_energies(
 
 
 def run_inference(model: ASEModel, test_data: Path | None = None) -> dict[str, dict]:
-    """Evaluate curvature roughness on homonuclear dimers."""
+    """Compare model and PBE dissociation curves with the stacking_fault metric."""
     label_path = _LABEL_FILE if test_data is None else test_data / "diatomics.json"
 
     with open(label_path) as fh:
@@ -124,13 +169,8 @@ def run_inference(model: ASEModel, test_data: Path | None = None) -> dict[str, d
         mol_name: str = entry["name"]
         bond_lengths, dft_energies = _scan_arrays(entry)
         element = _element_from_name(mol_name)
-        dr = float(np.mean(np.diff(bond_lengths)))
-
         model_energies = _predict_energies(calc, element, bond_lengths)
-
-        mol_result = {
-            "roughness": _compute_roughness(model_energies - dft_energies, dr),
-        }
+        mol_result = _curve_metrics(bond_lengths, model_energies, dft_energies)
         results[mol_name] = mol_result
         logging.info(f"{mol_name}: {mol_result}")
 
diff --git a/tests/tasks/calculator/test_diatomics.py b/tests/tasks/calculator/test_diatomics.py
index 3b8f520..fb97684 100644
--- a/tests/tasks/calculator/test_diatomics.py
+++ b/tests/tasks/calculator/test_diatomics.py
@@ -2,7 +2,7 @@
 import pytest
 
 from lambench.tasks.calculator.diatomics.diatomics import (
-    _compute_roughness,
+    _curve_metrics,
     _element_from_name,
     _scan_arrays,
 )
@@ -18,45 +18,76 @@ def _full_results(**overrides: dict | None) -> dict:
     return results
 
 
-def test_compute_roughness_flat_residuals():
-    residuals = np.zeros(10)
-    assert _compute_roughness(residuals, dr=0.1) == pytest.approx(0.0)
+def _well_scan():
+    r = np.linspace(0.8, 4.0, 25)
+    e = (r - 1.4) ** 2 - 1.0
+    return r, e
 
 
-def test_compute_roughness_oscillating():
-    residuals = np.array([0.0, 1.0, 0.0, 1.0, 0.0, 1.0], dtype=float)
-    roughness = _compute_roughness(residuals, dr=0.2)
-    assert roughness is not None
-    assert roughness > 0
+def test_scan_arrays_drops_last_point():
+    r, e = _scan_arrays(
+        {"name": "SiSi", "R": [1.0, 1.2, 1.4, 1.6], "E": [0.0, -1.0, -0.5, 10.0]}
+    )
+    np.testing.assert_array_equal(r, [1.0, 1.2, 1.4])
+    np.testing.assert_array_equal(e, [0.0, -1.0, -0.5])
 
 
-def test_compute_roughness_dr_scaling():
-    residuals = np.array([0.0, 0.1, 0.0, 0.1, 0.0], dtype=float)
-    r1 = _compute_roughness(residuals, dr=0.2)
-    r2 = _compute_roughness(residuals, dr=0.1)
-    assert r2 == pytest.approx(4 * r1, rel=1e-6)
+def test_perfect_match_is_zero():
+    r, e = _well_scan()
+    metrics = _curve_metrics(r, e, e)
+    assert metrics is not None
+    assert metrics["roughness"] == pytest.approx(0.0)
 
 
-def test_compute_roughness_too_few_valid_points():
-    residuals = np.array([np.nan, 0.5, np.nan])
-    assert _compute_roughness(residuals, dr=0.1) is None
+def test_constant_offset_is_removed_by_min_shift():
+    r, e = _well_scan()
+    metrics = _curve_metrics(r, e + 1.5, e)
+    assert metrics is not None
+    assert metrics["roughness"] == pytest.approx(0.0, abs=1e-10)
 
 
-def test_compute_roughness_all_nan():
-    assert _compute_roughness(np.array([np.nan, np.nan, np.nan]), dr=0.1) is None
+def test_constant_dummy_scores_one():
+    r, e = _well_scan()
+    dummy = _curve_metrics(r, np.full_like(e, e[-1]), e)
+    match = _curve_metrics(r, e, e)
+    assert dummy is not None and match is not None
+    assert dummy["roughness"] == pytest.approx(1.0)
+    assert dummy["roughness"] > match["roughness"]
 
 
-def test_compute_roughness_matching_raw_reference_is_zero():
-    residuals = np.zeros(20)
-    assert _compute_roughness(residuals, dr=0.2) == pytest.approx(0.0)
+def test_worse_than_dummy_is_capped_at_one():
+    r, e = _well_scan()
+    y = e - np.min(e)
+    inverted = _curve_metrics(r, np.max(y) - y, e)
+    assert inverted is not None
+    assert inverted["roughness"] == pytest.approx(1.0)
 
 
-def test_scan_arrays_drops_last_point():
-    r, e = _scan_arrays(
-        {"name": "SiSi", "R": [1.0, 1.2, 1.4, 1.6], "E": [0.0, -1.0, -0.5, 10.0]}
-    )
-    np.testing.assert_array_equal(r, [1.0, 1.2, 1.4])
-    np.testing.assert_array_equal(e, [0.0, -1.0, -0.5])
+def test_zero_bond_step_returns_none():
+    r = np.array([1.0, 1.0, 1.4])
+    e = np.array([0.0, -1.0, -0.5])
+    assert _curve_metrics(r, e, e) is None
+
+
+def test_shape_match_is_scale_invariant():
+    r, e = _well_scan()
+    scaled = _curve_metrics(r, 10 * e, e)
+    assert scaled is not None
+    assert scaled["roughness"] == pytest.approx(0.0, abs=1e-10)
+
+
+def test_oscillation_increases_roughness():
+    r, e = _well_scan()
+    match = _curve_metrics(r, e, e)
+    wiggly = _curve_metrics(r, e + 0.2 * np.sin(25 * (r - r[0])), e)
+    assert match is not None and wiggly is not None
+    assert wiggly["roughness"] > match["roughness"]
+
+
+def test_nonfinite_returns_none():
+    r = np.array([1.0, 1.2, 1.4, 1.6])
+    e = np.array([0.0, np.nan, -0.5, -0.4])
+    assert _curve_metrics(r, e, e) is None
 
 
 def test_element_from_name():
@@ -76,9 +107,9 @@ def test_aggregated_means():
     results = {name: {"roughness": 0.02} for name in names}
     results["HH"] = {"roughness": 0.01}
     results["NN"] = {"roughness": 0.03}
-    expected = (0.01 + 0.03 + 0.02 * (len(names) - 2)) / len(names)
+    n = len(names)
     agg = aggregated_diatomics_results(results)
-    assert agg["avg_roughness"] == pytest.approx(expected)
+    assert agg["avg_roughness"] == pytest.approx((0.01 + 0.03 + 0.02 * (n - 2)) / n)
 
 
 def test_aggregated_empty_results():
@@ -88,12 +119,6 @@ def test_aggregated_empty_results():
 
 def test_aggregated_incomplete_coverage_is_none():
     assert aggregated_diatomics_results(_full_results(HH=None))["avg_roughness"] is None
-    assert (
-        aggregated_diatomics_results(_full_results(HH={"roughness": np.nan}))[
-            "avg_roughness"
-        ]
-        is None
-    )
     incomplete = _full_results()
     del incomplete["HH"]
     assert aggregated_diatomics_results(incomplete)["avg_roughness"] is None

From 5b5e21d1457acd697dea973a535b5df727793a22 Mon Sep 17 00:00:00 2001
From: Anyang Peng <137014849+anyangml@users.noreply.github.com>
Date: Thu, 27 Aug 2026 15:07:46 +0800
Subject: [PATCH 2/2] chore:lint

---
 lambench/tasks/calculator/diatomics/diatomics.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/lambench/tasks/calculator/diatomics/diatomics.py b/lambench/tasks/calculator/diatomics/diatomics.py
index 542030e..dda5419 100644
--- a/lambench/tasks/calculator/diatomics/diatomics.py
+++ b/lambench/tasks/calculator/diatomics/diatomics.py
@@ -103,7 +103,10 @@ def _curve_metrics(
         return None
     if bond_lengths.size < 2:
         return None
-    if bond_lengths.size != model_energies.size or bond_lengths.size != dft_energies.size:
+    if (
+        bond_lengths.size != model_energies.size
+        or bond_lengths.size != dft_energies.size
+    ):
         return None
     if np.any(np.diff(bond_lengths) == 0):
         return None