"""Cross-validate H2O / RHF / 6-31G* against PySCF and (optionally) ORCA.

The "Hello, validation" of vibe-qc cross-checking. Same molecule, same
basis, three calculators, same number to ~µHa.

Run:
    .venv/bin/python examples/ase_compare/compare-h2o-hf.py

ORCA is optional — provide its path via the ``ORCA_COMMAND`` (or
``ASE_ORCA_COMMAND``) environment variable, or have ``orca`` on
``$PATH``. If unavailable, the ORCA row is reported as "unavailable"
and the comparison continues with vibe-qc + PySCF.

Produces:
    output-compare-h2o-hf.csv   — full comparison table
    (stdout)                    — pretty-printed table

Expected outcome: all three rows agree on E and |F|max to ~µHa
precision (the basis is identical and the SCF threshold is tight
enough). If they disagree, that's a real bug worth investigating —
either in vibe-qc or in the calling convention.
"""

from __future__ import annotations

from pathlib import Path

from ase.build import molecule

from vibeqc.ase import VibeQC
from vibeqc.benchmark import (
    PySCFCalculator,
    compare_calculators,
    make_orca_calculator,
    print_calculator_availability,
)

HERE = Path(__file__).resolve().parent


def main() -> None:
    print("=" * 72)
    print(" Cross-validation:  H2O / RHF / 6-31G*")
    print("=" * 72)
    print()
    print("Calculator availability:")
    print_calculator_availability()
    print()

    atoms = molecule("H2O")  # ASE's stock H2O geometry (Angstrom)

    calculators: list = [
        ("vibe-qc/RHF/6-31G*", VibeQC(basis="6-31g*")),
        ("PySCF/RHF/6-31G*",   PySCFCalculator(basis="6-31g*", method="RHF")),
    ]

    orca = make_orca_calculator(
        orcasimpleinput="HF 6-31G* EnGrad",
        label="orca-h2o-hf",
    )
    if orca is not None:
        calculators.append(("ORCA/RHF/6-31G*", orca))
    else:
        print("(ORCA not found — set ORCA_COMMAND or add orca to $PATH "
              "to include it in the comparison)")
        print()

    results = compare_calculators(
        atoms,
        calculators,
        properties=("energy", "forces"),
    )
    results.print_table()
    csv_path = HERE / "output-compare-h2o-hf.csv"
    results.to_csv(csv_path)
    print(f"\nFull table written to {csv_path.relative_to(HERE.parent.parent)}")

    # Tight tolerance — these are the same SCF on the same basis. Any
    # disagreement past 1e-5 eV (~0.4 µHa) is a real problem.
    results.assert_agreement(
        reference="vibe-qc/RHF/6-31G*",
        tol={"energy": 1e-5, "forces": 1e-4},
    )
    print("✓ all calculators agree to within 1e-5 eV / 1e-4 eV·Å⁻¹ — "
          "identical SCF on identical basis, as expected.")


if __name__ == "__main__":
    main()