diff --git a/lambench/metrics/utils.py b/lambench/metrics/utils.py
index 6c9fc74..acb4940 100644
--- a/lambench/metrics/utils.py
+++ b/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
 
 
 ####################################
diff --git a/lambench/metrics/vishelper/metrics_calculations.py b/lambench/metrics/vishelper/metrics_calculations.py
index 0fe5d19..1346fc7 100644
--- a/lambench/metrics/vishelper/metrics_calculations.py
+++ b/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],
         )
diff --git a/lambench/tasks/calculator/diatomics/diatomics.py b/lambench/tasks/calculator/diatomics/diatomics.py
index dda5419..8cdbe16 100644
--- a/lambench/tasks/calculator/diatomics/diatomics.py
+++ b/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,
+    }
diff --git a/tests/metrics/test_metrics_calculator.py b/tests/metrics/test_metrics_calculator.py
index 9344e1e..01e8793 100644
--- a/tests/metrics/test_metrics_calculator.py
+++ b/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)
diff --git a/tests/tasks/calculator/test_diatomics.py b/tests/tasks/calculator/test_diatomics.py
index 9a17c29..d63af22 100644
--- a/tests/tasks/calculator/test_diatomics.py
+++ b/tests/tasks/calculator/test_diatomics.py
@@ -1,7 +1,12 @@
+import json
+from pathlib import Path
+from types import SimpleNamespace
+
 import numpy as np
 import pytest
 
 from lambench.metrics.utils import (
+    _diatomic_dummy_score,
     _diatomics_molecule_names,
     aggregated_diatomics_results,
 )
@@ -9,6 +14,7 @@
     _curve_metrics,
     _element_from_name,
     _scan_arrays,
+    run_inference,
 )
 
 
@@ -122,3 +128,138 @@ def test_aggregated_incomplete_coverage_is_none():
     incomplete = _full_results()
     del incomplete["HH"]
     assert aggregated_diatomics_results(incomplete)["avg_roughness"] is None
+
+
+def test_bond_position_and_depth_errors():
+    r, e = _well_scan()
+    exact = _curve_metrics(r, e, e)
+    assert exact["bond_length_error"] == pytest.approx(0.0)
+    assert exact["well_depth_error"] == pytest.approx(0.0)
+    shifted = _curve_metrics(r, (r - 1.6) ** 2 - 1.0, e)
+    assert shifted["bond_length_error"] == pytest.approx(0.2)
+    scaled = _curve_metrics(r, e * 2, e)
+    assert scaled["well_depth_error"] == pytest.approx(e[-1] - e.min())
+    assert scaled["roughness"] == pytest.approx(0.0, abs=1e-10)
+
+
+def test_constant_energy_dummy_scores_one_for_all_terms():
+    r, e = _well_scan()
+    record = _curve_metrics(r, np.zeros_like(e), e)
+    agg = aggregated_diatomics_results(
+        {"reference_molecules": ["HH"], "curves": {"HH": record}}
+    )
+    assert agg["avg_roughness"] == pytest.approx(1.0)
+    assert agg["bond_length_score"] == pytest.approx(1.0)
+    assert agg["well_depth_score"] == pytest.approx(1.0)
+    assert agg["score"] == pytest.approx(1.0)
+
+
+def test_weak_reference_still_penalizes_false_binding():
+    r, e = _well_scan()
+    weak = e * 0.001
+    result = _curve_metrics(r, e, weak)
+    assert result["bond_length_error"] is None
+    assert result["well_depth_error"] > 1.0
+
+
+def test_per_molecule_normalization_precedes_averaging():
+    curves = {
+        "HH": {
+            "roughness": 0.3,
+            "bond_length_error": 0.1,
+            "bond_length_dummy_error": 0.2,
+            "well_depth_error": 100.0,
+            "well_depth_dummy_error": 1.0,
+        },
+        "LiLi": {
+            "roughness": 0.1,
+            "bond_length_error": 0.0,
+            "bond_length_dummy_error": 0.1,
+            "well_depth_error": 0.0,
+            "well_depth_dummy_error": 1.0,
+        },
+    }
+    result = aggregated_diatomics_results(
+        {"reference_molecules": list(curves), "curves": curves}
+    )
+    assert result["bond_length_score"] == pytest.approx(0.25)
+    assert result["well_depth_score"] == pytest.approx(0.5)
+    assert result["score"] == pytest.approx((0.2 + 0.25 + 0.5) / 3)
+
+
+def test_inference_uses_trimmed_custom_reference_and_roundtrips(tmp_path, monkeypatch):
+    r, e = _well_scan()
+    reference = [{"name": "HH", "R": r.tolist() + [5.0], "E": e.tolist() + [999.0]}]
+    (tmp_path / "diatomics.json").write_text(json.dumps(reference))
+
+    def predict(calc, element, distances):
+        assert element == "H"
+        np.testing.assert_array_equal(distances, r)
+        return e.copy()
+
+    monkeypatch.setattr(
+        "lambench.tasks.calculator.diatomics.diatomics._predict_energies", predict
+    )
+    result = json.loads(json.dumps(run_inference(SimpleNamespace(calc=None), tmp_path)))
+    assert result["reference_molecules"] == ["HH"]
+    assert aggregated_diatomics_results(result)["score"] == pytest.approx(0.0)
+    del result["curves"]["HH"]
+    assert aggregated_diatomics_results(result)["score"] is None
+
+
+def test_builtin_references_self_score_zero():
+    path = (
+        Path(__file__).resolve().parents[3]
+        / "lambench/tasks/calculator/diatomics/diatomics.json"
+    )
+    reference = json.loads(path.read_text())
+    assert len(reference) == 14
+    curves = {}
+    for entry in reference:
+        r, e = _scan_arrays(entry)
+        curves[entry["name"]] = _curve_metrics(r, e, e)
+    result = aggregated_diatomics_results(
+        {"reference_molecules": list(curves), "curves": curves}
+    )
+    assert result["score"] == pytest.approx(0.0)
+    dummy_curves = {}
+    for entry in reference:
+        r, e = _scan_arrays(entry)
+        dummy_curves[entry["name"]] = _curve_metrics(r, np.zeros_like(e), e)
+    dummy = aggregated_diatomics_results(
+        {"reference_molecules": list(dummy_curves), "curves": dummy_curves}
+    )
+    assert dummy["score"] == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize(("error", "expected"), [(0.0, 0.0), (0.1, 1.0)])
+def test_zero_dummy_error_is_well_defined(error, expected):
+    assert _diatomic_dummy_score(error, 0.0) == expected
+
+
+def test_legacy_records_do_not_get_a_partial_combined_score():
+    result = aggregated_diatomics_results(_full_results())
+    assert result["avg_roughness"] == pytest.approx(0.02)
+    assert result["score"] is None
+
+
+def test_missing_bound_metric_does_not_get_a_score():
+    r, e = _well_scan()
+    result = _curve_metrics(r, e, e)
+    result["bond_length_error"] = None
+    agg = aggregated_diatomics_results(
+        {"reference_molecules": ["HH"], "curves": {"HH": result}}
+    )
+    assert agg["score"] is None
+
+
+@pytest.mark.parametrize("invalid", ["nan", "reversed", "matrix"])
+def test_invalid_distance_grid_returns_none(invalid):
+    r, e = _well_scan()
+    if invalid == "nan":
+        r[2] = np.nan
+    elif invalid == "reversed":
+        r = r[::-1]
+    else:
+        r = r[:, None]
+    assert _curve_metrics(r, e, e) is None