Add bond-length and well-depth metrics to existing diatomic references

lambench/metrics/utils.py

@@ -170,34 +170,88 @@ def _diatomics_molecule_names() -> tuple[str, ...]:
 
 
 def _empty_diatomics_agg() -> dict[str, float | None]:
-    return {"avg_roughness": None}
+    return {
+        "avg_roughness": None,
+        "bond_length_mae": None,
+        "well_depth_mae": None,
+        "bond_length_score": None,
+        "well_depth_score": None,
+        "score": None,
+    }
+
 
+def _diatomic_dummy_score(error: float, dummy_error: float) -> float:
+    # A zero-error dummy is already exact: only an exact prediction matches it.
+    if dummy_error == 0:
+        return 0.0 if error == 0 else 1.0
+    return float(min(error / dummy_error, 1.0))
 
-def aggregated_diatomics_results(results: dict[str, dict]) -> dict[str, float]:
+
+def aggregated_diatomics_results(results: dict) -> dict[str, float | None]:
     """
     Aggregate per-molecule diatomics results.
 
     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.
+    New results carry the tested molecule set and per-molecule dummy errors.
+    The three scores are averaged with equal weight only for complete results.
+    Legacy roughness-only records remain readable but have no combined score.
     """
     if not results:
         return _empty_diatomics_agg()
 
-    names = _diatomics_molecule_names()
+    names = results.get("reference_molecules", _diatomics_molecule_names())
+    curves = results.get("curves", results)
     if not names:
         return _empty_diatomics_agg()
 
     roughness_values = []
+    errors = {"bond_length": [], "well_depth": []}
+    scores = {"bond_length": [], "well_depth": []}
+    complete_metrics = True
     for name in names:
-        mol_results = results.get(name)
+        mol_results = curves.get(name)
         if mol_results is None:
             return _empty_diatomics_agg()
         roughness = mol_results.get("roughness")
-        if roughness is None or not np.isfinite(roughness):
+        if roughness is None or not np.isfinite(roughness) or not 0 <= roughness <= 1:
             return _empty_diatomics_agg()
         roughness_values.append(roughness)
-
-    return {"avg_roughness": float(np.mean(roughness_values))}
+        for metric, metric_errors in errors.items():
+            error = mol_results.get(f"{metric}_error")
+            dummy = mol_results.get(f"{metric}_dummy_error")
+            if metric == "bond_length" and error is None and dummy is None:
+                # Only references below the binding threshold skip bond position.
+                depth = mol_results.get("well_depth_dummy_error")
+                if depth is None or depth >= 0.05:
+                    complete_metrics = False
+                continue
+            if (
+                error is None
+                or dummy is None
+                or not np.isfinite(error)
+                or not np.isfinite(dummy)
+                or error < 0
+                or dummy < 0
+            ):
+                complete_metrics = False
+                continue
+            metric_errors.append(error)
+            scores[metric].append(_diatomic_dummy_score(error, dummy))
+
+    aggregated = _empty_diatomics_agg()
+    aggregated["avg_roughness"] = float(np.mean(roughness_values))
+    for metric, metric_errors in errors.items():
+        if metric_errors:
+            aggregated[f"{metric}_mae"] = float(np.mean(metric_errors))
+            aggregated[f"{metric}_score"] = float(np.mean(scores[metric]))
+    terms = [
+        aggregated["avg_roughness"],
+        aggregated["bond_length_score"],
+        aggregated["well_depth_score"],
+    ]
+    if complete_metrics and all(term is not None for term in terms):
+        aggregated["score"] = float(np.mean(terms))
+    return aggregated
 
 
 ####################################

lambench/metrics/vishelper/metrics_calculations.py

@@ -199,17 +199,17 @@ def _calculate_instability_error(self, cell: dict, lambda_0: float = 5e-4) -> fl
             else:
                 return np.clip(np.log10(slope / lambda_0), a_min=0, a_max=None)
 
-    def calculate_diatomics_roughness_results(self) -> dict[str, float]:
+    def calculate_diatomics_results(self) -> dict[str, float]:
         """
-        Leaderboard scores for homonuclear diatomics: avg_roughness (lower is better).
+        Leaderboard scores for homonuclear diatomics: combined score (lower is better).
         Models with missing results are omitted from this dict; the final ranking
-        still includes them with a null roughness column.
+        still includes them with a null diatomics column.
         """
         raw = self.fetcher.fetch_diatomics_results()
         return {
-            model: metrics["avg_roughness"]
+            model: metrics["score"]
             for model, metrics in raw.items()
-            if metrics is not None and metrics.get("avg_roughness") is not None
+            if metrics is not None and metrics.get("score") is not None
         }
 
     def calculate_efficiency_results(self) -> dict[str, float]:
@@ -238,7 +238,7 @@ def summarize_final_rankings(self):
         )
         stability_results = self.calculate_stability_results()
         efficiency_results = self.calculate_efficiency_results()
-        roughness_results = self.calculate_diatomics_roughness_results()
+        diatomics_results = self.calculate_diatomics_results()
         if not generalizability_ood or not generalizability_downstream:
             logging.warning(
                 "Missing data for generalizability metrics (ood or downstream)"
@@ -272,8 +272,8 @@ def summarize_final_rankings(self):
             "Applicability-Efficiency ↑": [
                 efficiency_results[model] for model in shared_models
             ],
-            "Applicability-Roughness ↓": [
-                roughness_results.get(model) for model in shared_models
+            "Applicability-Diatomics ↓": [
+                diatomics_results.get(model) for model in shared_models
             ],
         }
 
@@ -293,7 +293,7 @@ def summarize_final_rankings(self):
                 "Generalizability-PC Error ↓",
                 "Applicability-Instability ↓",
                 "Applicability-Efficiency ↑",
-                "Applicability-Roughness ↓",
+                "Applicability-Diatomics ↓",
             ],
             ascending=[True, True, True, False, True],
         )

lambench/tasks/calculator/diatomics/diatomics.py

@@ -26,7 +26,17 @@
 
     0 = perfect match, 1 = dummy or worse.
 
-    roughness – slope MAE relative to dummy; leaderboard
+    roughness – slope MAE relative to dummy
+
+Two additional metrics use the same trimmed scan: equilibrium bond-length
+error (angstrom) and well-depth error (eV), with depth defined relative to
+the last retained separation rather than an infinite-distance atomic limit.
+Bond length is fitted to up to five points around the sampled minimum and
+skipped when reference depth is below 0.05 eV. Its dummy guesses the scan
+midpoint; the well-depth dummy predicts constant energy (zero well depth).
+Each error is dummy-normalized and capped per molecule, then averaged.
+The leaderboard uses the equal-weight mean of these two scores and roughness,
+only when every reference molecule has a valid result.
 
 The last point of every scan is dropped (PBE tail artifacts).
 
@@ -50,6 +60,8 @@
 
 _LABEL_FILE = Path(__file__).parent / "diatomics.json"
 _DROP_LAST_POINTS = 1
+_MIN_BINDING_EV = 0.05
+_FIT_POINTS = 5
 
 
 def _element_from_name(mol_name: str) -> str:
@@ -96,9 +108,29 @@ def _ratio_capped_at_dummy(value: float, dummy: float) -> float | None:
     return float(min(value / dummy, 1.0))
 
 
+def _equilibrium_distance(bond_lengths: np.ndarray, energies: np.ndarray) -> float:
+    """Fit the local well, falling back to the sampled minimum if not convex."""
+    if np.ptp(energies) == 0:
+        return float((bond_lengths[0] + bond_lengths[-1]) / 2)
+    minimum = int(np.argmin(energies))
+    sampled = float(bond_lengths[minimum])
+    if len(bond_lengths) < 3 or minimum in (0, len(bond_lengths) - 1):
+        return sampled
+    start = min(max(0, minimum - _FIT_POINTS // 2), max(0, len(energies) - _FIT_POINTS))
+    r = bond_lengths[start : start + _FIT_POINTS]
+    e = energies[start : start + _FIT_POINTS]
+    a, b, _ = np.polyfit(r, e - energies[minimum], 2)
+    if a <= 0:
+        return sampled
+    fitted = -b / (2 * a)
+    return float(fitted) if r[0] <= fitted <= r[-1] else sampled
+
+
 def _curve_metrics(
     bond_lengths: np.ndarray, model_energies: np.ndarray, dft_energies: np.ndarray
-) -> dict[str, float] | None:
+) -> dict[str, float | None] | None:
+    if any(values.ndim != 1 for values in (bond_lengths, model_energies, dft_energies)):
+        return None
     if not np.all(np.isfinite(model_energies)) or not np.all(np.isfinite(dft_energies)):
         return None
     if bond_lengths.size < 2:
@@ -108,7 +140,7 @@ def _curve_metrics(
         or bond_lengths.size != dft_energies.size
     ):
         return None
-    if np.any(np.diff(bond_lengths) == 0):
+    if not np.isfinite(bond_lengths).all() or np.any(np.diff(bond_lengths) <= 0):
         return None
 
     y_dft = _shift_to_min(dft_energies)
@@ -127,7 +159,22 @@ def _curve_metrics(
     )
     if roughness is None:
         return None
-    return {"roughness": roughness}
+    ref_depth = float(y_dft[-1])
+    model_depth = float(y_model[-1])
+    result: dict[str, float | None] = {
+        "roughness": roughness,
+        "bond_length_error": None,
+        "bond_length_dummy_error": None,
+        "well_depth_error": abs(model_depth - ref_depth),
+        "well_depth_dummy_error": ref_depth,
+    }
+    if ref_depth >= _MIN_BINDING_EV:
+        ref_distance = _equilibrium_distance(bond_lengths, dft_energies)
+        model_distance = _equilibrium_distance(bond_lengths, model_energies)
+        midpoint = float((bond_lengths[0] + bond_lengths[-1]) / 2)
+        result["bond_length_error"] = abs(model_distance - ref_distance)
+        result["bond_length_dummy_error"] = abs(midpoint - ref_distance)
+    return result
 
 
 def _predict_energies(
@@ -158,14 +205,14 @@ def _predict_energies(
     return np.array(energies)
 
 
-def run_inference(model: ASEModel, test_data: Path | None = None) -> dict[str, dict]:
-    """Compare model and PBE dissociation curves with the stacking_fault metric."""
+def run_inference(model: ASEModel, test_data: Path | None = None) -> dict:
+    """Score trimmed scans and save the exact reference molecule set."""
     label_path = _LABEL_FILE if test_data is None else test_data / "diatomics.json"
 
     with open(label_path) as fh:
         reference_data: list[dict] = json.load(fh)
 
-    results: dict[str, dict] = {}
+    results: dict[str, dict | None] = {}
     calc = model.calc
 
     for entry in reference_data:
@@ -177,4 +224,7 @@ def run_inference(model: ASEModel, test_data: Path | None = None) -> dict[str, d
         results[mol_name] = mol_result
         logging.info(f"{mol_name}: {mol_result}")
 
-    return results
+    return {
+        "reference_molecules": [entry["name"] for entry in reference_data],
+        "curves": results,
+    }

tests/metrics/test_metrics_calculator.py

@@ -92,7 +92,7 @@ def test_summarize_final_rankings(metrics_calculator):
     metrics_calculator.calculate_stability_results = MagicMock(
         return_value={"model1": 0.2, "model2": 0.5}
     )
-    metrics_calculator.calculate_diatomics_roughness_results = MagicMock(
+    metrics_calculator.calculate_diatomics_results = MagicMock(
         return_value={"model1": 0.05, "model2": 0.03}
     )
     _, result = metrics_calculator.summarize_final_rankings()
@@ -101,14 +101,14 @@ def test_summarize_final_rankings(metrics_calculator):
     assert result.iloc[1]["Model"] == "model1"
 
 
-def test_calculate_diatomics_roughness_results(metrics_calculator, mock_raw_results):
+def test_calculate_diatomics_results(metrics_calculator, mock_raw_results):
     mock_raw_results.fetch_diatomics_results.return_value = {
-        "model1": {"avg_roughness": 0.05},
-        "model2": {"avg_roughness": 0.03},
+        "model1": {"score": 0.05},
+        "model2": {"score": 0.03},
         "model3": None,
-        "model4": {"avg_roughness": None},
+        "model4": {"score": None},
     }
-    result = metrics_calculator.calculate_diatomics_roughness_results()
+    result = metrics_calculator.calculate_diatomics_results()
     assert set(result.keys()) == {"model1", "model2"}
     np.testing.assert_almost_equal(result["model1"], 0.05)
     np.testing.assert_almost_equal(result["model2"], 0.03)