diff --git a/.github/workflows/scientific-report.yml b/.github/workflows/scientific-report.yml new file mode 100644 index 00000000..772fa06d --- /dev/null +++ b/.github/workflows/scientific-report.yml @@ -0,0 +1,122 @@ +name: Scientific and performance report + +on: + workflow_dispatch: + inputs: + eri_filter: + description: Substring filter for ERI cases + required: true + default: h2-sto3g + type: string + mp2_sizes: + description: Comma-separated MP2 AO dimensions + required: true + default: '12,24' + type: string + mp2_memory: + description: Configured MP2 workspace budget + required: true + default: 32 MiB + type: string + include_heavy_eri: + description: Include the heavy ERI cases + required: true + default: false + type: boolean + +permissions: + contents: read + +jobs: + report: + name: Scientific validation and Criterion reports + runs-on: ubuntu-latest + env: + RUSTIQ_BENCH_FILTER: ${{ inputs.eri_filter }} + RUSTIQ_BENCH_HEAVY: ${{ inputs.include_heavy }} + RUSTIQ_MP2_SIZES: ${{ inputs.mp2_sizes }} + RUSTIQ_MP2_MEMORY: ${{ inputs.mp2_memory }} + + steps: + - name: Check out selected commit + uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@stable + + - name: Set up uv + uses: astral-sh/setup-uv@c18668ad3cf93ea998bef934396af7bb5c839dc7 # v10.2.0 + + - name: Install locked Python environment + run: uv sync --locked + + - name: Remove stale Criterion measurements + run: rm -rf target/criterion + + - name: Run canonical PySCF reference suite + id: reference + continue-on-error: true + run: uv run --locked pytest tools/reference + + - name: Generate scientific report + id: scientific + if: always() + continue-on-error: true + run: uv run --locked python tools/reference/report.py --output-dir target/scientific-report + + - name: Run selected ERI Criterion benchmarks + id: eri + continue-on-error: true + run: cargo bench --features bench-support --bench eri_timings + + - name: Run selected MP2 Criterion benchmarks + id: mp2 + continue-on-error: true + run: cargo bench --features bench-support --bench mp2_timings + + - name: Generate performance report + id: performance + if: always() + continue-on-error: true + run: uv run --locked python tools/benchmarks/report.py --criterion-dir target/criterion --output-dir target/performance-report + + - name: Upload generated reports + if: always() + uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2 + with: + name: rustiq-scientific-performance-reports + if-no-files-found: warn + path: | + target/scientific-report/ + target/performance-report/ + + - name: Write job summary + if: always() + env: + REFERENCE: ${{ steps.reference.outcome }} + SCIENTIFIC: ${{ steps.scientific.outcome }} + ERI: ${{ steps.eri.outcome }} + MP2: ${{ steps.mp2.outcome }} + PERFORMANCE: ${{ steps.performance.outcome }} + run: | + { + echo '# RustiQ scientific and performance reports' + echo + echo '| Stage | Result |' + echo '|---|---|' + echo "| Canonical PySCF suite | $REFERENCE |" + echo "| Scientific report | $SCIENTIFIC |" + echo "| ERI Criterion | $ERI |" + echo "| MP2 Criterion | $MP2 |" + echo "| Performance report | $PERFORMANCE |" + echo + echo 'Scientific validation and performance failures are reported independently. Download the artifact for detailed diagnostics.' + } >> "$GITHUB_STEP_SUMMARY" + if [[ "$REFERENCE" != success || "$SCIENTIFIC" != success ]]; then + echo 'Scientific validation failed; inspect the PySCF test and validation report steps.' >&2 + exit 1 + fi + if [[ "$ERI" != success || "$MP2" != success || "$PERFORMANCE" != success ]]; then + echo 'Performance reporting failed; inspect the Criterion and performance report steps.' >&2 + exit 1 + fi diff --git a/benches/eri_timings.rs b/benches/eri_timings.rs index abb72a0e..a6a0d63e 100644 --- a/benches/eri_timings.rs +++ b/benches/eri_timings.rs @@ -59,6 +59,12 @@ fn eri_benchmarks(criterion: &mut Criterion) { let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR")); let basis_store = basis_store(); let mut group = criterion.benchmark_group("eri"); + let mut metadata = Vec::new(); + let run_metadata = run_metadata(); + let path = PathBuf::from("target/criterion/eri-metadata.json"); + if path.exists() { + std::fs::remove_file(&path).expect("stale ERI metadata should be removable"); + } for case in selected_cases() { // Ten is Criterion's minimum sample size; use it for heavy systems to // keep the requested measurement target from multiplying into hours. @@ -82,6 +88,19 @@ fn eri_benchmarks(criterion: &mut Criterion) { let result = input .run_once() .expect("benchmark input should produce valid Coulomb self-integrals"); + metadata.push(serde_json::json!({ + "benchmark_id": format!("eri/{}/{}", case.name, case.basis), + "case": case.name, + "molecule": case.geometry, + "basis": case.basis, + "basis_functions": result.basis_functions, + "compact_integrals": result.compact_integrals, + "pair_count": result.pair_count, + "heavy": case.heavy, + "criterion_sample_size": if case.heavy { 10 } else { 30 }, + "criterion_warm_up_seconds": if case.heavy { 1 } else { 2 }, + "criterion_measurement_seconds": if case.heavy { 2 } else { 5 }, + })); eprintln!( "ERI {}: {} basis functions, {} pairs, {} compact integrals, {:?} pair expansion, {:?} Schwarz bounds, {:?} fill, {} cached Coulomb values", case.name, @@ -112,6 +131,57 @@ fn eri_benchmarks(criterion: &mut Criterion) { } group.finish(); + if let Some(parent) = path.parent() { + std::fs::create_dir_all(parent).expect("Criterion metadata directory should be writable"); + } + std::fs::write( + path, + serde_json::to_vec_pretty(&serde_json::json!({"run": run_metadata, "cases": metadata})) + .expect("metadata should serialize"), + ) + .expect("Criterion metadata should be writable"); +} + +fn run_metadata() -> serde_json::Value { + let output = std::process::Command::new("git") + .args(["rev-parse", "HEAD"]) + .output() + .ok(); + let commit = output + .filter(|result| result.status.success()) + .map(|result| String::from_utf8_lossy(&result.stdout).trim().to_owned()); + let compiler = std::process::Command::new("rustc") + .arg("--version") + .output() + .ok(); + let rustc = compiler + .filter(|result| result.status.success()) + .map(|result| String::from_utf8_lossy(&result.stdout).trim().to_owned()); + let clean = std::process::Command::new("git") + .args(["status", "--porcelain"]) + .output() + .ok() + .filter(|result| result.status.success()) + .map(|result| result.stdout.is_empty()); + let cpu = std::fs::read_to_string("/proc/cpuinfo") + .ok() + .and_then(|info| { + info.lines() + .find(|line| line.starts_with("model name\t:")) + .map(str::to_owned) + }) + .and_then(|line| { + line.split_once(':') + .map(|(_, value)| value.trim().to_owned()) + }); + let timestamp = std::time::SystemTime::now() + .duration_since(std::time::UNIX_EPOCH) + .ok() + .map(|duration| duration.as_secs_f64()); + serde_json::json!({"git_commit": commit, "working_tree_clean": clean, "rustc_version": rustc, + "operating_system": std::env::consts::OS, "architecture": std::env::consts::ARCH, + "cpu_model": cpu, "benchmark_started_unix_seconds": timestamp, + "benchmark_profile": "cargo bench default bench profile"}) } fn selected_cases() -> Vec<&'static BenchCase> { diff --git a/benches/mp2_timings.rs b/benches/mp2_timings.rs index 0d362bd4..d4428667 100644 --- a/benches/mp2_timings.rs +++ b/benches/mp2_timings.rs @@ -7,6 +7,10 @@ use rustiq_core::bench_support::Mp2BenchInput; use rustiq_core::config::MemoryLimit; fn mp2_benchmarks(criterion: &mut Criterion) { + let path = std::path::PathBuf::from("target/criterion/mp2-metadata.json"); + if path.exists() { + std::fs::remove_file(&path).expect("stale MP2 metadata should be removable"); + } let requested = std::env::var("RUSTIQ_MP2_MEMORY").unwrap_or_else(|_| "auto".into()); let budget = if requested.eq_ignore_ascii_case("auto") { MemoryLimit::Auto.resolve() @@ -18,6 +22,8 @@ fn mp2_benchmarks(criterion: &mut Criterion) { .split(',') .map(|size| size.parse().expect("integer AO dimension")) .collect(); + let mut metadata = Vec::new(); + let run_metadata = run_metadata(); let mut group = criterion.benchmark_group("mp2"); for n in sizes { @@ -63,9 +69,73 @@ fn mp2_benchmarks(criterion: &mut Criterion) { budget.display().iec(), rayon::current_num_threads() ); + metadata.push(serde_json::json!({ + "ao_dimension": n, + "occupied_orbitals": occupied, + "workspace_budget_bytes": budget.as_u64(), + "block_size": memory.block_size, + "rayon_threads": rayon::current_num_threads(), + "dense_energy": dense_energy, + "blocked_energy": blocked_energy, + "energy_absolute_difference": (dense_energy - blocked_energy).abs(), + "criterion_sample_size": if n >= 48 { 10 } else { 30 }, + "criterion_warm_up_seconds": if n >= 48 { 1 } else { 2 }, + "criterion_measurement_seconds": if n >= 48 { 2 } else { 5 }, + })); } group.finish(); + if let Some(parent) = path.parent() { + std::fs::create_dir_all(parent).expect("Criterion metadata directory should be writable"); + } + std::fs::write( + path, + serde_json::to_vec_pretty(&serde_json::json!({"run": run_metadata, "cases": metadata})) + .expect("metadata should serialize"), + ) + .expect("Criterion metadata should be writable"); +} + +fn run_metadata() -> serde_json::Value { + let output = std::process::Command::new("git") + .args(["rev-parse", "HEAD"]) + .output() + .ok(); + let commit = output + .filter(|result| result.status.success()) + .map(|result| String::from_utf8_lossy(&result.stdout).trim().to_owned()); + let compiler = std::process::Command::new("rustc") + .arg("--version") + .output() + .ok(); + let rustc = compiler + .filter(|result| result.status.success()) + .map(|result| String::from_utf8_lossy(&result.stdout).trim().to_owned()); + let clean = std::process::Command::new("git") + .args(["status", "--porcelain"]) + .output() + .ok() + .filter(|result| result.status.success()) + .map(|result| result.stdout.is_empty()); + let cpu = std::fs::read_to_string("/proc/cpuinfo") + .ok() + .and_then(|info| { + info.lines() + .find(|line| line.starts_with("model name\t:")) + .map(str::to_owned) + }) + .and_then(|line| { + line.split_once(':') + .map(|(_, value)| value.trim().to_owned()) + }); + let timestamp = std::time::SystemTime::now() + .duration_since(std::time::UNIX_EPOCH) + .ok() + .map(|duration| duration.as_secs_f64()); + serde_json::json!({"git_commit": commit, "working_tree_clean": clean, "rustc_version": rustc, + "operating_system": std::env::consts::OS, "architecture": std::env::consts::ARCH, + "cpu_model": cpu, "benchmark_started_unix_seconds": timestamp, + "benchmark_profile": "cargo bench default bench profile"}) } criterion_group!(benches, mp2_benchmarks); diff --git a/docs/src/development/contributing.md b/docs/src/development/contributing.md index b5b2a04e..b2cbe0d0 100644 --- a/docs/src/development/contributing.md +++ b/docs/src/development/contributing.md @@ -81,7 +81,7 @@ CLI TOML/Nickel adapter or another frontend -> iterative SCF and final canonicalization -> HfOutcome::Converged or HfOutcome::Unconverged -> converged HfSolution::mp2() - blocked AO-to-MO transformation and correlation energy + blocked AO-to-MO integral transformation and MP2 correlation-energy evaluation ``` `prepare()` constructs reusable input data; it does not compute integrals or @@ -143,9 +143,10 @@ cargo bench --bench mp2_timings --features bench-support cargo bench --bench eri_timings --features bench-support ``` -The MP2 benchmark compares dense and memory-blocked kernels on deterministic -synthetic inputs. It validates that both paths produce matching energies before -timing them; input and ERI construction remain outside the measured iterations. +The MP2 benchmark compares dense and memory-blocked AO-to-MO transformations +and their subsequent correlation-energy evaluation on deterministic synthetic +inputs. It validates that both paths produce matching energies before timing +them; input and ERI construction remain outside the measured iterations. `RUSTIQ_MP2_MEMORY`, `RUSTIQ_MP2_SIZES` (comma-separated AO dimensions), and `RAYON_NUM_THREADS` select the memory budget, AO dimensions, and thread count. The ERI benchmark times repeated integral construction for the checked-in @@ -182,3 +183,53 @@ cargo bench --bench mp2_timings --features bench-support -- --baseline before See the [MP2 benchmark source](https://github.com/mveril/RustiQ/blob/main/benches/mp2_timings.rs) and [ERI benchmark source](https://github.com/mveril/RustiQ/blob/main/benches/eri_timings.rs). Benchmarks provide performance measurements, not scientific validation. + +### Generate reproducible reports + +From the repository root, the scientific report reuses the canonical PySCF +reference cases and their absolute tolerances: + +```sh +uv run --locked python tools/reference/report.py --output-dir target/scientific-report +``` + +It writes full precision results to `validation.json` and readable tables to +`validation.md`. In addition to the cases run by `pytest tools/reference`, it +includes the frozen-core MP2 variants exercised by those tests at 32 KiB and +512 MiB workspace budgets. Total HF energies and MP2 correlation energies are +reported separately. A case fails if either implementation does not converge, +an energy comparison exceeds its existing tolerance, or the UHF spin expectation +value misses its tolerance. Failed cases remain in the report, and the command +returns a nonzero status. The report records software versions, commit and +platform metadata so results can be compared across revisions. + +Generate Criterion measurements and their report with: + +```sh +cargo bench --features bench-support --bench eri_timings +cargo bench --features bench-support --bench mp2_timings +uv run --locked python tools/benchmarks/report.py +``` + +The report is written to `target/performance-report/benchmarks.json` and +`benchmarks.md`. To produce both reports on GitHub, open **Actions → Scientific +and performance report → Run workflow** on the selected branch or commit. The +workflow uploads `rustiq-scientific-performance-reports`; download that artifact +from the completed workflow run to inspect the Markdown and JSON files. + +Criterion stores structured estimates in `target/criterion`; the benchmark +executables also write the case dimensions and configuration used by the report +generator there. The default runs include the small ERI cases and MP2 AO sizes +12, 24, and 48. To select a smaller local run or opt into the expensive ERI +cases, set `RUSTIQ_BENCH_FILTER`, `RUSTIQ_MP2_SIZES`, `RUSTIQ_MP2_MEMORY`, and +`RUSTIQ_BENCH_HEAVY` as described above. The report includes Criterion means +and confidence intervals. Its configured MP2 workspace is not measured peak +RSS, and synthetic MP2 cases do not model complete molecular calculations. +Results are machine-dependent; the on-demand GitHub Actions workflow starts +with one small ERI case and AO sizes 12 and 24. Heavy ERI cases require explicit +opt-in in the workflow. + +These reports provide numerical regression checks for selected cases. They do +not establish general chemical accuracy; see the [scientific scope](../scientific-scope.md) +for the limits of HF, MP2, and numerical agreement. Equal numeric convergence +thresholds do not mean that RustiQ and PySCF use identical convergence tests. diff --git a/tools/benchmarks/report.py b/tools/benchmarks/report.py new file mode 100644 index 00000000..89f87477 --- /dev/null +++ b/tools/benchmarks/report.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 +"""Summarize structured Criterion results and metadata emitted by RustiQ benches.""" +from __future__ import annotations + +import argparse +import json +import math +import re +import subprocess +import tomllib +from datetime import UTC, datetime +from pathlib import Path +from typing import Any + +ROOT = Path(__file__).resolve().parents[2] + + +class ReportError(ValueError): + """Criterion data is absent or cannot support the requested report.""" + + +def read_json(path: Path) -> Any: + try: + return json.loads(path.read_text(encoding="utf-8")) + except (OSError, json.JSONDecodeError) as error: + raise ReportError(f"Invalid or missing JSON: {path}: {error}") from error + + +def finite(value: Any, label: str) -> float: + try: + result = float(value) + except (TypeError, ValueError) as error: + raise ReportError(f"{label} must be numeric") from error + if not math.isfinite(result): + raise ReportError(f"{label} must be finite") + return result + + +def criterion_measurement( + root: Path, benchmark_id: str, not_before: float +) -> dict[str, Any]: + # Criterion 0.8 emits structured benchmark and estimate JSON under new/. + encoded = benchmark_id.replace("/", "_") + for path in root.rglob("new/benchmark.json"): + data = read_json(path) + if not isinstance(data, dict): + raise ReportError(f"Criterion benchmark metadata must be a JSON object: {path}") + if data.get("full_id") == benchmark_id or path.parent.parent.name == encoded: + estimates_path = path.parent / "estimates.json" + try: + if estimates_path.stat().st_mtime < not_before: + raise ReportError( + f"Criterion results for {benchmark_id} predate this benchmark run" + ) + except OSError as error: + raise ReportError(f"Missing Criterion estimates: {estimates_path}") from error + estimates = read_json(estimates_path) + if not isinstance(estimates, dict): + raise ReportError(f"Criterion estimates must be a JSON object: {path.parent / 'estimates.json'}") + mean = estimates.get("mean", {}) + if not isinstance(mean, dict): + raise ReportError(f"Criterion mean estimate is invalid for {benchmark_id}") + interval = mean.get("confidence_interval", {}) + if not isinstance(interval, dict): + raise ReportError(f"Criterion confidence interval is invalid for {benchmark_id}") + estimate = finite(mean.get("point_estimate"), f"{benchmark_id} mean") + lower = finite(interval.get("lower_bound"), f"{benchmark_id} CI lower") + upper = finite(interval.get("upper_bound"), f"{benchmark_id} CI upper") + confidence = finite(interval.get("confidence_level"), f"{benchmark_id} CI confidence") + if estimate <= 0 or lower < 0 or lower > upper or not lower <= estimate <= upper or not 0 < confidence < 1: + raise ReportError(f"Invalid confidence interval for {benchmark_id}") + return {"benchmark_id": benchmark_id, "mean_seconds": estimate / 1e9, + "confidence_interval_seconds": [lower / 1e9, upper / 1e9], + "confidence_level": confidence, "throughput": data.get("throughput")} + raise ReportError(f"Missing Criterion results for {benchmark_id} under {root}") + + +def machine_metadata(eri_run: dict[str, Any], mp2_run: dict[str, Any]) -> dict[str, Any]: + fields = ("git_commit", "working_tree_clean", "rustc_version", "operating_system", "architecture", "cpu_model", "benchmark_profile") + for key in fields: + if eri_run.get(key) != mp2_run.get(key): + raise ReportError(f"ERI and MP2 benchmark provenance differs for {key}") + if not eri_run.get("git_commit") or not eri_run.get("rustc_version"): + raise ReportError("Benchmark execution did not record a Git commit and Rust compiler version") + try: + current_commit = subprocess.run(["git", "rev-parse", "HEAD"], cwd=ROOT, text=True, + capture_output=True, check=True).stdout.strip() + current_status = subprocess.run(["git", "status", "--porcelain"], cwd=ROOT, text=True, + capture_output=True, check=True).stdout + except (OSError, subprocess.CalledProcessError) as error: + raise ReportError(f"Could not verify benchmark source revision: {error}") from error + if eri_run["git_commit"] != current_commit: + raise ReportError("Criterion measurements were made at a different Git commit") + if eri_run.get("working_tree_clean") != (current_status == ""): + raise ReportError("Benchmark working-tree state differs from the report-generation source tree") + try: + started = min(finite(eri_run["benchmark_started_unix_seconds"], "ERI timestamp"), + finite(mp2_run["benchmark_started_unix_seconds"], "MP2 timestamp")) + except KeyError as error: + raise ReportError(f"Missing benchmark execution timestamp: {error}") from error + try: + cargo = tomllib.loads((ROOT / "Cargo.toml").read_text(encoding="utf-8")) + version = cargo["package"]["version"] + except (OSError, KeyError, tomllib.TOMLDecodeError): + version = None + return {"rustiq_version": version, **{key: eri_run[key] for key in fields}, + "benchmark_started_at_utc": datetime.fromtimestamp(started, UTC).isoformat(), + "eri_benchmark_started_at_utc": datetime.fromtimestamp(finite(eri_run["benchmark_started_unix_seconds"], "ERI timestamp"), UTC).isoformat(), + "mp2_benchmark_started_at_utc": datetime.fromtimestamp(finite(mp2_run["benchmark_started_unix_seconds"], "MP2 timestamp"), UTC).isoformat(), + "report_generated_at_utc": datetime.now(UTC).isoformat(), + "criterion_version": _criterion_version(), "results_machine_dependent": True} + + +def _criterion_version() -> str | None: + try: + lock = (ROOT / "Cargo.lock").read_text(encoding="utf-8") + except OSError: + return None + match = re.search(r'(?ms)^\[\[package\]\]\nname = "criterion"\nversion = "([^\"]+)"', lock) + return match.group(1) if match else None + + +def generate(criterion_dir: Path, output_dir: Path) -> dict[str, Any]: + eri_file = read_json(criterion_dir / "eri-metadata.json") + mp2_file = read_json(criterion_dir / "mp2-metadata.json") + if not isinstance(eri_file, dict) or not isinstance(mp2_file, dict): + raise ReportError("ERI and MP2 metadata must be JSON objects") + eri_meta, mp2_meta = eri_file.get("cases"), mp2_file.get("cases") + eri_run, mp2_run = eri_file.get("run"), mp2_file.get("run") + if not isinstance(eri_meta, list) or not eri_meta or not isinstance(mp2_meta, list) or not mp2_meta: + raise ReportError("ERI and MP2 benchmark metadata must both contain executed cases") + if not isinstance(eri_run, dict) or not isinstance(mp2_run, dict): + raise ReportError("Missing benchmark execution provenance") + eri = [] + for item in eri_meta: + if not isinstance(item, dict) or not isinstance(item.get("benchmark_id"), str): + raise ReportError("Invalid ERI benchmark metadata entry") + measurement = criterion_measurement( + criterion_dir, + item["benchmark_id"], + finite(eri_run.get("benchmark_started_unix_seconds"), "ERI timestamp"), + ) + throughput = measurement.pop("throughput") + eps = None + if isinstance(throughput, dict) and len(throughput) == 1: + eps = finite(next(iter(throughput.values())), "ERI throughput") / measurement["mean_seconds"] + eri.append({**item, **measurement, "throughput_elements_per_second": eps}) + mp2 = [] + for item in mp2_meta: + if not isinstance(item, dict) or not isinstance(item.get("ao_dimension"), int): + raise ReportError("Invalid MP2 benchmark metadata entry") + n = item["ao_dimension"] + mp2_started = finite(mp2_run.get("benchmark_started_unix_seconds"), "MP2 timestamp") + dense = criterion_measurement(criterion_dir, f"mp2/dense/{n}", mp2_started) + blocked = criterion_measurement(criterion_dir, f"mp2/blocked/{n}", mp2_started) + if finite(item.get("energy_absolute_difference"), f"MP2 N={n} consistency") >= 1e-12: + raise ReportError(f"Dense and blocked MP2 energy consistency check failed for N={n}") + if finite(item.get("workspace_budget_bytes"), "workspace budget") <= 0 or item.get("block_size", 0) <= 0: + raise ReportError(f"Invalid MP2 configuration metadata for N={n}") + mp2.append({**item, "dense": dense, "blocked": blocked, + "dense_over_blocked_ratio": dense["mean_seconds"] / blocked["mean_seconds"]}) + report = {"schema_version": 1, "metadata": machine_metadata(eri_run, mp2_run), + "criterion_configuration": { + "eri": [{key: item[key] for key in ("benchmark_id", "criterion_sample_size", "criterion_warm_up_seconds", "criterion_measurement_seconds")} for item in eri_meta], + "mp2": [{key: item[key] for key in ("ao_dimension", "criterion_sample_size", "criterion_warm_up_seconds", "criterion_measurement_seconds")} for item in mp2_meta]}, + "rayon_threads": sorted({item["rayon_threads"] for item in mp2}), + "eri": eri, "mp2": mp2, + "memory_note": "Configured MP2 workspace budget is not measured peak process memory. Peak RSS was not independently measured.", + "mp2_scope_note": "Synthetic MP2 cases time AO-to-MO transformation and correlation-energy evaluation; they do not represent complete molecular calculations."} + output_dir.mkdir(parents=True, exist_ok=True) + (output_dir / "benchmarks.json").write_text(json.dumps(report, indent=2, sort_keys=True, allow_nan=False) + "\n", encoding="utf-8") + (output_dir / "benchmarks.md").write_text(render_markdown(report), encoding="utf-8") + return report + + +def render_markdown(report: dict[str, Any]) -> str: + fmt = lambda value: f"{value:.6g}" + lines = ["# RustiQ Performance Benchmarks", "", "Results are machine-dependent and come from Criterion measurements. No comparison against other quantum chemistry packages is made.", "", "## Reproducibility information", "", "```json", json.dumps(report["metadata"], indent=2, sort_keys=True), "```", "", "Sample sizes, warm-up and measurement targets are recorded per case in `benchmarks.json`.", "", "## ERI computation", "", "| Case | Basis | AO functions | Compact integrals | Mean (s) | Criterion CI (s) | Compact elements/s |", "|---|---|---:|---:|---:|---:|---:|"] + for row in report["eri"]: + lo, hi = row["confidence_interval_seconds"] + level = round(row["confidence_level"] * 100) + throughput = fmt(row["throughput_elements_per_second"]) if row["throughput_elements_per_second"] is not None else "—" + lines.append(f"| {row['case']} | {row['basis']} | {row['basis_functions']} | {row['compact_integrals']} | {fmt(row['mean_seconds'])} | {level}% [{fmt(lo)}, {fmt(hi)}] | {throughput} |") + lines += ["", "Throughput counts compact storage elements per second; Schwarz screening may skip evaluations or leave stored values zero.", "", "## MP2 AO-to-MO transformation and correlation-energy evaluation", "", "Synthetic AO-dimension cases time the transformation and correlation-energy evaluation, not full molecular calculations.", "", "| AO dim | Occupied | Dense mean (s) | Dense CI (s) | Blocked mean (s) | Blocked CI (s) | Dense/blocked | Workspace budget | Block size | Rayon threads |", "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|"] + for row in report["mp2"]: + dlo, dhi = row["dense"]["confidence_interval_seconds"] + blo, bhi = row["blocked"]["confidence_interval_seconds"] + dlevel = round(row["dense"]["confidence_level"] * 100) + blevel = round(row["blocked"]["confidence_level"] * 100) + lines.append(f"| {row['ao_dimension']} | {row['occupied_orbitals']} | {fmt(row['dense']['mean_seconds'])} | {dlevel}% [{fmt(dlo)}, {fmt(dhi)}] | {fmt(row['blocked']['mean_seconds'])} | {blevel}% [{fmt(blo)}, {fmt(bhi)}] | {fmt(row['dense_over_blocked_ratio'])} | {row['workspace_budget_bytes']} B | {row['block_size']} | {row['rayon_threads']} |") + lines += ["", "## Memory interpretation", "", report["memory_note"], "", "Configured workspace is a per-calculation budget for the blocked MP2 transformation; it is not measured peak process memory.", ""] + return "\n".join(lines) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--criterion-dir", type=Path, default=Path("target/criterion")) + parser.add_argument("--output-dir", type=Path, default=Path("target/performance-report")) + args = parser.parse_args() + try: + report = generate(args.criterion_dir, args.output_dir) + except ReportError as error: + parser.error(str(error)) + print(f"Performance report written for {len(report['eri'])} ERI and {len(report['mp2'])} MP2 cases") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/benchmarks/test_bench_report.py b/tools/benchmarks/test_bench_report.py new file mode 100644 index 00000000..e8f16c53 --- /dev/null +++ b/tools/benchmarks/test_bench_report.py @@ -0,0 +1,76 @@ +from __future__ import annotations + +import importlib.util +import json +from pathlib import Path + +import pytest + +SPEC = importlib.util.spec_from_file_location("performance_report", Path(__file__).with_name("report.py")) +assert SPEC and SPEC.loader +report = importlib.util.module_from_spec(SPEC) +SPEC.loader.exec_module(report) + + +def write_criterion(root: Path, benchmark_id: str, mean: float = 2e6) -> None: + directory = root / "group" / benchmark_id.replace("/", "_") / "new" + directory.mkdir(parents=True) + (directory / "benchmark.json").write_text(json.dumps({"full_id": benchmark_id, + "throughput": {"Elements": 16}})) + (directory / "estimates.json").write_text(json.dumps({"mean": { + "point_estimate": mean, "confidence_interval": {"confidence_level": 0.95, + "lower_bound": mean * 0.9, "upper_bound": mean * 1.1}}})) + + +def test_reads_criterion_estimates_and_confidence_interval(tmp_path: Path) -> None: + write_criterion(tmp_path, "eri/h2-sto3g/sto-3g") + result = report.criterion_measurement(tmp_path, "eri/h2-sto3g/sto-3g", 0) + assert result["mean_seconds"] == 0.002 + assert result["confidence_interval_seconds"] == [0.0018, 0.0022] + + +def test_invalid_and_missing_criterion_results_fail(tmp_path: Path) -> None: + with pytest.raises(report.ReportError, match="Missing Criterion"): + report.criterion_measurement(tmp_path, "mp2/dense/12", 0) + write_criterion(tmp_path, "mp2/dense/12", float("nan")) + with pytest.raises(report.ReportError, match="finite"): + report.criterion_measurement(tmp_path, "mp2/dense/12", 0) + + +def test_stale_criterion_results_are_rejected(tmp_path: Path) -> None: + write_criterion(tmp_path, "eri/h2-sto3g/sto-3g") + with pytest.raises(report.ReportError, match="predate this benchmark run"): + report.criterion_measurement(tmp_path, "eri/h2-sto3g/sto-3g", 4e9) + + +def test_invalid_json_is_diagnostic(tmp_path: Path) -> None: + path = tmp_path / "invalid.json" + path.write_text("not json") + with pytest.raises(report.ReportError, match="Invalid or missing JSON"): + report.read_json(path) + + +def test_performance_markdown_separates_budget_from_peak_memory() -> None: + row = {"case": "h2", "basis": "sto-3g", "basis_functions": 2, + "compact_integrals": 6, "mean_seconds": 0.01, + "confidence_interval_seconds": [0.009, 0.011], + "confidence_level": 0.95, + "throughput_elements_per_second": 600.0} + mp = {"ao_dimension": 12, "occupied_orbitals": 4, "dense": {"mean_seconds": 0.02, + "confidence_interval_seconds": [0.019, 0.021], "confidence_level": 0.95}, "blocked": {"mean_seconds": 0.01, + "confidence_interval_seconds": [0.009, 0.011], "confidence_level": 0.95}, "dense_over_blocked_ratio": 2.0, + "workspace_budget_bytes": 32768, "block_size": 2, "rayon_threads": 4} + md = report.render_markdown({"metadata": {}, "eri": [row], "mp2": [mp], + "memory_note": "Configured budget; peak RSS was not measured."}) + assert "Dense/blocked" in md + assert "peak RSS was not measured" in md + assert "Compact elements/s" in md + + +def test_execution_provenance_must_match() -> None: + run = {"git_commit": "abc", "working_tree_clean": True, "rustc_version": "rustc test", "operating_system": "linux", + "architecture": "x86_64", "cpu_model": "cpu", "benchmark_profile": "bench", + "benchmark_started_unix_seconds": 1.0} + other = {**run, "git_commit": "different"} + with pytest.raises(report.ReportError, match="provenance differs"): + report.machine_metadata(run, other) diff --git a/tools/reference/README.md b/tools/reference/README.md index c1041ec1..6849e442 100644 --- a/tools/reference/README.md +++ b/tools/reference/README.md @@ -110,6 +110,22 @@ refresh a fixture while keeping the source explicit, point hashes used by both uv and Nix; update it intentionally with `uv lock` whenever the Python dependency declarations change. +Generate a full precision scientific validation report from the same cases and +comparisons: + +```sh +uv run --locked python tools/reference/report.py --output-dir target/scientific-report +``` + +The command writes `validation.json` and `validation.md`. It includes the +canonical cases plus the two frozen-core MP2 cases at each memory budget used +by pytest. A failed or non-finite calculation remains visible in the report and +causes a nonzero exit status. The JSON contains machine-readable values; +Markdown formats energies for reading. The comparison only validates selected +numerical cases against PySCF. It does not establish general chemical accuracy +or make HF/MP2 reliable for every system. Matching thresholds also do not imply +identical convergence criteria. + RustiQ is invoked with `run --format json`; pytest reads schema version 1 JSON directly, checks convergence, methods, and energies, and provides all comparison reporting. This keeps reference validation independent of terminal report wording diff --git a/tools/reference/compare_pyscf.py b/tools/reference/compare_pyscf.py old mode 100755 new mode 100644 index 0be7fbe4..94334690 --- a/tools/reference/compare_pyscf.py +++ b/tools/reference/compare_pyscf.py @@ -9,6 +9,7 @@ import json import os +import re import shutil import subprocess import sys @@ -36,6 +37,19 @@ class ReferenceCase: mp2: bool = False mp2_tolerance: float | None = None frozen_orbitals: int = 0 + frozen_core_reference: bool = False + + +class PyscfMp2Error(RuntimeError): + """MP2 failed after converged HF; retain the completed HF diagnostics.""" + + def __init__(self, case_name: str, hf_result: tuple[float, None, float | None], cause: Exception): + super().__init__(f"PySCF MP2 failed for {case_name}: {cause}") + self.hf_result = hf_result + + +class PyscfScfConvergenceError(RuntimeError): + """PySCF reached its cycle limit without converging HF.""" CASES = [ @@ -166,6 +180,7 @@ class ReferenceCase: spin_tolerance=2e-5, mp2=True, mp2_tolerance=5e-9, + frozen_core_reference=True, ), ReferenceCase( name="h2-sto-3g-rhf-mp2", @@ -213,9 +228,30 @@ class ReferenceCase: # The residual correlation-energy difference is about 5.1e-9 Hartree # with the checked-in BSE basis and a 1e-12 SCF threshold. mp2_tolerance=1e-8, + frozen_core_reference=True, ), ] +FROZEN_CORE_MEMORY_BUDGETS = ("32 KiB", "512 MiB") +RUSTIQ_SCHEMA_VERSION = 1 + + +def validate_rustiq_result(result: object, case_name: str) -> dict[str, object]: + """Check the stable JSON envelope consumed by reference validation.""" + if not isinstance(result, dict): + raise TypeError(f"RustiQ JSON output for {case_name} must be an object") + version = result.get("schema_version") + if type(version) is not int or version != RUSTIQ_SCHEMA_VERSION: + raise ValueError( + f"Unsupported or missing RustiQ schema_version for {case_name}: {version!r}" + ) + calculation = result.get("calculation") + if not isinstance(calculation, dict): + raise TypeError(f"RustiQ JSON output for {case_name} is missing calculation object") + if not isinstance(calculation.get("hf"), dict): + raise TypeError(f"RustiQ JSON output for {case_name} is missing calculation.hf object") + return result + def run_command(args: list[str], *, env: dict[str, str] | None = None) -> str: completed = subprocess.run( @@ -278,11 +314,12 @@ def rustiq_result(case: ReferenceCase, env: dict[str, str]) -> dict[str, object] env=env, ) try: - return json.loads(stdout) + parsed = json.loads(stdout) except json.JSONDecodeError as error: raise RuntimeError( f"Could not parse RustiQ JSON output for {case.name}: {error}\nstdout:\n{stdout}" ) from error + return validate_rustiq_result(parsed, case.name) def pyscf_result(case: ReferenceCase) -> tuple[float, float | None, float | None]: @@ -307,11 +344,16 @@ def pyscf_result(case: ReferenceCase) -> tuple[float, float | None, float | None mf.max_cycle = case.max_cycle energy = mf.kernel() if not mf.converged: - raise RuntimeError(f"PySCF did not converge for {case.name}.") + raise PyscfScfConvergenceError(f"PySCF did not converge for {case.name}.") s_squared = float(mf.spin_square()[0]) if case.method == "uhf" else None mp2_correlation_energy = None if case.mp2: - mp2_correlation_energy, _ = mp.MP2(mf, frozen=case.frozen_orbitals).kernel() + try: + mp2_correlation_energy, _ = mp.MP2(mf, frozen=case.frozen_orbitals).kernel() + except Exception as error: + raise PyscfMp2Error( + case.name, (float(energy), None, s_squared), error + ) from error return ( float(energy), (float(mp2_correlation_energy) if mp2_correlation_energy is not None else None), @@ -319,6 +361,37 @@ def pyscf_result(case: ReferenceCase) -> tuple[float, float | None, float | None ) +def frozen_core_case( + case: ReferenceCase, budget: str, directory: Path +) -> ReferenceCase: + """Construct the frozen-core variants shared by pytest and reporting.""" + source = case.runfile.read_text(encoding="utf-8") + source = re.sub( + r"^geometry = .*$", + f"geometry = {json.dumps(str(case.xyz))}", + source, + flags=re.MULTILINE, + ) + marker = "frozen_orbitals = 0" + if source.count(marker) != 1: + raise ValueError(f"Expected one MP2 frozen-orbital setting in {case.runfile}") + source = source.replace( + marker, f'frozen_orbitals = 1\nmemory_limit = "{budget}"' + ) + directory.mkdir(parents=True, exist_ok=True) + runfile = directory / f"{case.name}-frozen-mp2-{budget.replace(' ', '-')}.toml" + runfile.write_text(source, encoding="utf-8") + import tomllib + + generated = tomllib.loads(source) + mp2 = generated.get("method", {}).get("mp2", {}) + if mp2.get("frozen_orbitals") != 1 or mp2.get("memory_limit") != budget: + raise ValueError("Generated frozen-core runfile did not apply its requested settings") + from dataclasses import replace + + return replace(case, runfile=runfile, frozen_orbitals=1) + + def main(pytest_args: list[str] | None = None) -> int: """Run the canonical pytest reference suite, forwarding any arguments.""" import pytest diff --git a/tools/reference/report.py b/tools/reference/report.py new file mode 100644 index 00000000..f7187a8d --- /dev/null +++ b/tools/reference/report.py @@ -0,0 +1,273 @@ +#!/usr/bin/env python3 +"""Generate a full precision RustiQ/PySCF scientific validation report.""" +from __future__ import annotations + +import argparse +import json +import math +import platform +import shutil +import subprocess +import tomllib +from dataclasses import replace +from datetime import UTC, datetime +from pathlib import Path +from typing import Any + +import compare_pyscf + +REPO_ROOT = compare_pyscf.REPO_ROOT + + +def frozen_runfile(case: compare_pyscf.ReferenceCase, budget: str, directory: Path) -> compare_pyscf.ReferenceCase: + """Use the same frozen-core runfile construction as the pytest suite.""" + return compare_pyscf.frozen_core_case(case, budget, directory) + + +def number(value: Any) -> float | None: + if isinstance(value, bool): + return None + try: + result = float(value) + except (TypeError, ValueError): + return None + return result if math.isfinite(result) else None + + +def difference(left: Any, right: Any) -> float | None: + a, b = number(left), number(right) + return abs(a - b) if a is not None and b is not None else None + + +def passed(error: float | None, tolerance: float | None, *required: Any) -> bool: + return all(number(value) is not None for value in required) and error is not None and tolerance is not None and error <= tolerance + + +def validate_case(case: compare_pyscf.ReferenceCase, rustiq: dict[str, Any] | None, pyscf: tuple[float, float | None, float | None] | None, errors: list[str], pyscf_hf_converged: bool | None = None) -> dict[str, Any]: + if rustiq is not None: + try: + rustiq = compare_pyscf.validate_rustiq_result(rustiq, case.name) + except (TypeError, ValueError) as error: + errors.append(str(error)) + rustiq = None + calc: dict[str, Any] = {} + hf: dict[str, Any] = {} + mp2: dict[str, Any] | None = None + spin: dict[str, Any] | None = None + try: + calc = rustiq["calculation"] if rustiq is not None else {} + hf = calc["hf"] + if not isinstance(hf, dict): + hf = {} + except (KeyError, TypeError): + errors.append("RustiQ result is missing calculation.hf") + rhf_energy = number(hf.get("total_energy")) + py_hf = number(pyscf[0]) if pyscf else None + error_hf = difference(rhf_energy, py_hf) + converged = hf.get("converged") is True + hf_method = hf.get("method") + orthogonalization = hf.get("orthogonalization") + ao_dimension = orthogonalization.get("ao_basis_dimension") if isinstance(orthogonalization, dict) else None + method_matches = hf_method == case.method.upper() + dimension_matches = ao_dimension == case.ao_dimension + if not method_matches: + errors.append("RustiQ HF method does not match the reference case") + if not dimension_matches: + errors.append("RustiQ AO basis dimension does not match the reference case") + py_converged = pyscf_hf_converged is True + hf_ok = (converged and py_converged and method_matches and dimension_matches + and passed(error_hf, case.tolerance, rhf_energy, py_hf)) + if not converged: + errors.append("RustiQ HF did not converge") + if pyscf_hf_converged is not True: + errors.append("PySCF HF did not converge" if pyscf_hf_converged is False + else "PySCF HF convergence status is unknown") + if py_hf is None: + errors.append("PySCF HF calculation did not return a finite energy; convergence could not be established") + if error_hf is None or (error_hf is not None and error_hf > case.tolerance): + errors.append("HF energy comparison failed") + result: dict[str, Any] = { + "case_id": case.name, "molecule": molecular_formula(case.xyz), + "geometry_path": str(case.xyz.relative_to(REPO_ROOT)), "basis": case.basis, + "charge": case.charge, "spin_n_alpha_minus_n_beta": case.spin, + "hf_method": hf_method, "expected_hf_method": case.method.upper(), + "ao_basis_dimension": ao_dimension, "expected_ao_basis_dimension": case.ao_dimension, + "rustiq_hf_total_energy_hartree": rhf_energy, + "pyscf_hf_total_energy_hartree": py_hf, "hf_absolute_difference_hartree": error_hf, + "hf_tolerance_hartree": case.tolerance, "rustiq_hf_converged": converged, + "pyscf_hf_converged": pyscf_hf_converged, "hf_pass": hf_ok, + "status": "FAIL", "errors": errors, + } + if case.mp2: + try: + mp2 = calc["mp2"] + except (KeyError, TypeError): + mp2 = None + py_mp2 = number(pyscf[1]) if pyscf else None + rust_mp2 = number(mp2.get("correlation_energy")) if isinstance(mp2, dict) else None + mp2_error = difference(rust_mp2, py_mp2) + tolerance = case.mp2_tolerance + expected_mp2_method = f"{case.method.upper()}-MP2" + mp2_method = mp2.get("method") if isinstance(mp2, dict) else None + mp2_ok = mp2_method == expected_mp2_method and passed(mp2_error, tolerance, rust_mp2, py_mp2) + if mp2_method != expected_mp2_method: + errors.append("RustiQ MP2 method does not match the reference case") + if rust_mp2 is None or py_mp2 is None: + errors.append("MP2 correlation energy is missing or non-finite") + if not mp2_ok: + errors.append("MP2 correlation energy comparison failed") + result.update({"mp2_reference_frozen_orbitals": case.frozen_orbitals, + "mp2_frozen_orbitals_convention": "spatial orbitals" if case.method == "rhf" else "orbitals in each spin channel", + "rustiq_mp2_method": mp2_method, + "rustiq_mp2_correlation_energy_hartree": rust_mp2, + "pyscf_mp2_correlation_energy_hartree": py_mp2, + "mp2_absolute_difference_hartree": mp2_error, + "mp2_tolerance_hartree": tolerance, + "pyscf_mp2_result_available": py_mp2 is not None, "mp2_pass": mp2_ok}) + if case.method == "uhf": + try: + spin = hf["spin"] + except (KeyError, TypeError): + spin = None + rust_spin = number(spin.get("s_squared")) if isinstance(spin, dict) else None + py_spin = number(pyscf[2]) if pyscf else None + spin_error = difference(rust_spin, py_spin) + spin_ok = passed(spin_error, case.spin_tolerance, rust_spin, py_spin) + if not spin_ok: + errors.append("Spin expectation comparison failed or is missing") + result.update({"rustiq_s_squared": rust_spin, "pyscf_s_squared": py_spin, + "s_squared_absolute_difference": spin_error, + "spin_tolerance": case.spin_tolerance, "spin_pass": spin_ok}) + elif isinstance(hf, dict) and hf.get("spin") is not None: + errors.append("RustiQ reported an unexpected open-shell spin diagnostic for RHF") + if not case.mp2: + rustiq_mp2 = calc.get("mp2") if isinstance(calc, dict) else None + pyscf_mp2 = pyscf[1] if pyscf is not None else None + if rustiq_mp2 is not None or pyscf_mp2 is not None: + errors.append("An unexpected MP2 result was returned for an HF-only case") + result["status"] = "PASS" if ( + not errors and hf_ok + and (not case.mp2 or result.get("mp2_pass") is True) + and (case.method != "uhf" or result.get("spin_pass") is True) + ) else "FAIL" + return result + + +def metadata() -> dict[str, Any]: + def command(*args: str) -> str | None: + try: + return subprocess.run(args, cwd=REPO_ROOT, text=True, capture_output=True, check=True).stdout.strip() + except (OSError, subprocess.CalledProcessError): + return None + try: + from pyscf import __version__ as pyscf_version + except ImportError: + pyscf_version = None + cargo = tomllib.loads((REPO_ROOT / "Cargo.toml").read_text(encoding="utf-8")) + return {"rustiq_version": cargo["package"]["version"], + "git_commit": command("git", "rev-parse", "HEAD"), + "working_tree_clean": command("git", "status", "--porcelain") == "", + "pyscf_version": pyscf_version, "python_version": platform.python_version(), + "rustc_version": command("rustc", "--version"), "operating_system": platform.platform(), + "cpu_architecture": platform.machine(), "generated_at_utc": datetime.now(UTC).isoformat()} + + +def molecular_formula(path: Path) -> str: + counts: dict[str, int] = {} + for line in path.read_text(encoding="utf-8").splitlines()[2:]: + fields = line.split() + if fields: + symbol = fields[0] + if symbol.isdecimal(): + # The shared H2 geometry uses XYZ's atomic-number spelling for hydrogen. + symbol = {"1": "H"}.get(symbol, symbol) + counts[symbol] = counts.get(symbol, 0) + 1 + if counts == {"H": 1, "O": 1}: + return "OH" + return "".join(symbol + (str(counts[symbol]) if counts[symbol] > 1 else "") for symbol in sorted(counts)) + + +def render_markdown(report: dict[str, Any]) -> str: + rows = report["cases"] + fmt = lambda n: "—" if n is None else f"{n:.10e}" + lines = ["# RustiQ Scientific Validation", "", "## Experimental status and limitations", "", + "Agreement with PySCF validates only the selected numerical cases; it does not establish general chemical accuracy. The reference suite is not exhaustive. HF omits electron correlation, and MP2 can be unreliable for strong correlation or small energy denominators. Matching convergence thresholds do not imply identical convergence criteria.", "", + "## Reproducibility information", "", "```json", json.dumps(report["metadata"], indent=2, sort_keys=True), "```", "", + "## Reference cases", "", "Spin is Nα−Nβ=2M_S. The AO dimension and method shown are read from RustiQ JSON and checked against the case definition.", "", + "| Case | Molecule | Geometry | Basis | Charge | Spin (Nα−Nβ) | AO dimension | Method |", "|---|---|---|---|---:|---:|---:|---|"] + for r in rows: + lines.append(f"| {r['case_id']} | {r['molecule']} | `{r['geometry_path']}` | {r['basis']} | {r['charge']} | {r['spin_n_alpha_minus_n_beta']} | {r['ao_basis_dimension'] if r['ao_basis_dimension'] is not None else '—'} | {r['hf_method'] or '—'} |") + lines += ["", "## Hartree-Fock validation", "", "Energies below are total molecular HF energies in Hartree, including nuclear repulsion.", "", + "| Case | Method | RustiQ total (Ha) | PySCF total (Ha) | Absolute error (Ha) | Tolerance (Ha) | Status |", "|---|---:|---:|---:|---:|---:|---|"] + for r in rows: + lines.append(f"| {r['case_id']} | {r['hf_method']} | {fmt(r['rustiq_hf_total_energy_hartree'])} | {fmt(r['pyscf_hf_total_energy_hartree'])} | {fmt(r['hf_absolute_difference_hartree'])} | {fmt(r['hf_tolerance_hartree'])} | {'PASS' if r['hf_pass'] else 'FAIL'} |") + lines += ["", "## MP2 validation", "", "MP2 values are correlation energies only; they exclude the HF energy and nuclear repulsion.", "", + "| Case | HF method | Frozen orbitals | RustiQ correlation (Ha) | PySCF correlation (Ha) | Absolute error (Ha) | Tolerance (Ha) | Status |", "|---|---|---:|---:|---:|---:|---:|---|"] + for r in rows: + if "mp2_pass" in r: + lines.append(f"| {r['case_id']} | {r['hf_method'] or '—'} | {r.get('mp2_reference_frozen_orbitals', 0)} | {fmt(r['rustiq_mp2_correlation_energy_hartree'])} | {fmt(r['pyscf_mp2_correlation_energy_hartree'])} | {fmt(r['mp2_absolute_difference_hartree'])} | {fmt(r['mp2_tolerance_hartree'])} | {'PASS' if r['mp2_pass'] else 'FAIL'} |") + lines += ["", "## Open-shell spin diagnostics", "", "| Case | RustiQ | PySCF | Absolute difference | Tolerance | Status |", "|---|---:|---:|---:|---:|---|"] + for r in rows: + if "spin_pass" in r: + lines.append(f"| {r['case_id']} | {fmt(r['rustiq_s_squared'])} | {fmt(r['pyscf_s_squared'])} | {fmt(r['s_squared_absolute_difference'])} | {fmt(r['spin_tolerance'])} | {'PASS' if r['spin_pass'] else 'FAIL'} |") + failed = [r for r in rows if r["status"] != "PASS"] + lines += ["", "## Summary", "", f"{len(rows) - len(failed)}/{len(rows)} cases passed; {len(failed)} failed."] + if failed: + lines += ["", "Failures:"] + for r in failed: + lines.append(f"- **{r['case_id']}**: {'; '.join(r['errors']) or 'validation failed'}") + lines += ["", "A frozen count applies to spatial orbitals for RHF and to orbitals in each spin channel for UHF; the exact convention is included in the JSON.", "", "## Scientific limitations", "", "This suite checks numerical agreement for the listed geometries, basis sets, charge and spin states, and methods. It does not assess basis-set completeness, broad chemical accuracy, or the suitability of HF/MP2 for arbitrary systems. The reported is the HF reference spin diagnostic.", ""] + return "\n".join(lines) + + +def generate(output_dir: Path) -> dict[str, Any]: + env = compare_pyscf.prepare_rustiq_env() + cases = list(compare_pyscf.CASES) + work = output_dir / ".inputs" + for original in (case for case in compare_pyscf.CASES if case.frozen_core_reference): + for budget in compare_pyscf.FROZEN_CORE_MEMORY_BUDGETS: + variant = frozen_runfile(original, budget, work) + cases.append(replace(variant, name=f"{original.name}-frozen1-{budget.replace(' ', '-')}")) + results = [] + try: + for case in cases: + errors: list[str] = [] + rustiq = pyscf = None + try: + rustiq = compare_pyscf.rustiq_result(case, env) + except Exception as error: # noqa: BLE001 - preserve every calculation failure in the report + errors.append(f"RustiQ calculation failed: {error}") + try: + pyscf = compare_pyscf.pyscf_result(case) + pyscf_hf_converged = True + except compare_pyscf.PyscfMp2Error as error: + pyscf = error.hf_result + pyscf_hf_converged = True + errors.append(str(error)) + except compare_pyscf.PyscfScfConvergenceError as error: + pyscf_hf_converged = False + errors.append(f"PySCF HF did not converge: {error}") + except Exception as error: # noqa: BLE001 - preserve every calculation failure in the report + errors.append(f"PySCF calculation failed: {error}") + results.append(validate_case(case, rustiq, pyscf, errors, pyscf_hf_converged)) + finally: + shutil.rmtree(Path(env["RUSTIQ_DATA_HOME"]), ignore_errors=True) + shutil.rmtree(work, ignore_errors=True) + report = {"schema_version": 1, "metadata": metadata(), "tolerances": {c.name: {"hf_hartree": c.tolerance, "mp2_hartree": c.mp2_tolerance, "spin_s_squared": c.spin_tolerance if c.method == "uhf" else None} for c in cases}, "reference_case_ids": [c.name for c in cases], "cases": results, "passed": all(r["status"] == "PASS" for r in results)} + output_dir.mkdir(parents=True, exist_ok=True) + (output_dir / "validation.json").write_text(json.dumps(report, indent=2, sort_keys=True, allow_nan=False) + "\n", encoding="utf-8") + (output_dir / "validation.md").write_text(render_markdown(report), encoding="utf-8") + return report + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output-dir", type=Path, default=Path("target/scientific-report")) + args = parser.parse_args() + report = generate(args.output_dir) + print(f"Scientific validation: {sum(r['status'] == 'PASS' for r in report['cases'])}/{len(report['cases'])} passed") + return 0 if report["passed"] else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/reference/test_compare_pyscf.py b/tools/reference/test_compare_pyscf.py index 11293601..0181b4c7 100644 --- a/tools/reference/test_compare_pyscf.py +++ b/tools/reference/test_compare_pyscf.py @@ -1,11 +1,9 @@ from __future__ import annotations import json -import re import shutil import tomllib from collections.abc import Iterator -from dataclasses import replace from pathlib import Path from typing import Any @@ -94,8 +92,30 @@ def malformed_output(*args: Any, **kwargs: Any) -> str: compare_pyscf.rustiq_result(compare_pyscf.CASES[0], {}) -@pytest.mark.parametrize("name", ["h2o-6-31g-rhf-mp2", "oh-sto-3g-uhf-mp2"]) -@pytest.mark.parametrize("budget", ["32 KiB", "512 MiB"]) +@pytest.mark.parametrize( + "payload", + [ + {"calculation": {"hf": {}}}, + {"schema_version": 2, "calculation": {"hf": {}}}, + {"schema_version": 1}, + {"schema_version": 1, "calculation": []}, + ], +) +def test_rustiq_result_rejects_incompatible_json_contract( + monkeypatch: pytest.MonkeyPatch, payload: dict[str, Any] +) -> None: + monkeypatch.setattr( + compare_pyscf, "run_command", lambda *args, **kwargs: json.dumps(payload) + ) + + with pytest.raises((TypeError, ValueError), match="RustiQ JSON output|schema_version"): + compare_pyscf.rustiq_result(compare_pyscf.CASES[0], {}) + + +@pytest.mark.parametrize( + "name", [case.name for case in compare_pyscf.CASES if case.frozen_core_reference] +) +@pytest.mark.parametrize("budget", compare_pyscf.FROZEN_CORE_MEMORY_BUDGETS) def test_frozen_core_mp2_matches_pyscf( name: str, budget: str, @@ -103,17 +123,5 @@ def test_frozen_core_mp2_matches_pyscf( rustiq_env: dict[str, str], ) -> None: original = next(case for case in compare_pyscf.CASES if case.name == name) - source = original.runfile.read_text() - source = re.sub( - r"^geometry = .*$", - f"geometry = {json.dumps(str(original.xyz))}", - source, - flags=re.MULTILINE, - ) - source = source.replace( - "frozen_orbitals = 0", f'frozen_orbitals = 1\nmemory_limit = "{budget}"' - ) - runfile = tmp_path / "frozen-mp2.toml" - runfile.write_text(source) - case = replace(original, runfile=runfile, frozen_orbitals=1) + case = compare_pyscf.frozen_core_case(original, budget, tmp_path) test_rustiq_matches_pyscf(case, rustiq_env) diff --git a/tools/reference/test_report.py b/tools/reference/test_report.py new file mode 100644 index 00000000..030a0073 --- /dev/null +++ b/tools/reference/test_report.py @@ -0,0 +1,156 @@ +from __future__ import annotations + +import json +import math +from pathlib import Path + +import compare_pyscf +import pytest +import report + + +def test_hf_mp2_and_spin_differences_and_markdown() -> None: + case = next(c for c in compare_pyscf.CASES if c.name == "oh-sto-3g-uhf-mp2") + rustiq = {"schema_version": 1, "calculation": {"hf": {"converged": True, "total_energy": -75.0, + "method": "UHF", "orthogonalization": {"ao_basis_dimension": case.ao_dimension}, + "spin": {"s_squared": 0.752}}, "mp2": {"method": "UHF-MP2", + "correlation_energy": -0.2}}} + result = report.validate_case(case, rustiq, (-75.0 + 1e-9, -0.2 + 1e-10, 0.752 + 1e-6), [], True) + assert math.isclose(result["hf_absolute_difference_hartree"], 1e-9, abs_tol=1e-12) + assert math.isclose(result["mp2_absolute_difference_hartree"], 1e-10, abs_tol=1e-12) + assert math.isclose(result["s_squared_absolute_difference"], 1e-6, abs_tol=1e-12) + assert result["status"] == "PASS" + md = report.render_markdown({"cases": [result], "metadata": {}}) + assert "total molecular HF energies" in md + assert "correlation energies only" in md + + +def test_missing_and_non_finite_values_fail() -> None: + case = compare_pyscf.CASES[0] + result = report.validate_case(case, {"schema_version": 1, "calculation": {"hf": {"converged": True, + "total_energy": float("nan"), "method": "RHF", + "orthogonalization": {"ao_basis_dimension": case.ao_dimension}}}}, (float("inf"), None, None), []) + assert result["rustiq_hf_total_energy_hartree"] is None + assert result["pyscf_hf_total_energy_hartree"] is None + assert result["status"] == "FAIL" + missing = report.validate_case(case, None, None, ["calculation failed"]) + assert missing["status"] == "FAIL" + assert missing["errors"] + + +def test_failed_tolerance_is_reported() -> None: + case = compare_pyscf.CASES[0] + row = report.validate_case(case, {"schema_version": 1, "calculation": {"hf": {"converged": True, + "total_energy": -1.0, "method": "RHF", + "orthogonalization": {"ao_basis_dimension": case.ao_dimension}}}}, (-0.9, None, None), []) + assert not row["hf_pass"] + assert row["status"] == "FAIL" + + +def valid_rustiq(case: compare_pyscf.ReferenceCase, *, energy: float = -75.0) -> dict: + hf = {"converged": True, "total_energy": energy, "method": case.method.upper(), + "orthogonalization": {"ao_basis_dimension": case.ao_dimension}} + calculation = {"hf": hf} + if case.method == "uhf": + hf["spin"] = {"s_squared": 0.75} + if case.mp2: + calculation["mp2"] = {"method": f"{case.method.upper()}-MP2", "correlation_energy": -0.2} + return {"schema_version": 1, "calculation": calculation} + + +@pytest.mark.parametrize("convergence", [False, None]) +def test_pyscf_nonconvergence_or_unknown_fails_even_when_energies_match(convergence) -> None: + case = compare_pyscf.CASES[0] + result = report.validate_case(case, valid_rustiq(case), (-75.0, None, None), [], convergence) + assert result["hf_absolute_difference_hartree"] == 0.0 + assert not result["hf_pass"] + assert result["status"] == "FAIL" + assert result["pyscf_hf_converged"] is convergence + + +@pytest.mark.parametrize("version", [None, 2, True]) +def test_missing_or_unsupported_schema_version_fails(version) -> None: + case = compare_pyscf.CASES[0] + rustiq = valid_rustiq(case) + if version is None: + rustiq.pop("schema_version") + else: + rustiq["schema_version"] = version + errors: list[str] = [] + result = report.validate_case(case, rustiq, (-75.0, None, None), errors, True) + assert result["status"] == "FAIL" + assert any("schema_version" in error for error in errors) + + +@pytest.mark.parametrize("payload", [{"schema_version": 1}, {"schema_version": 1, "calculation": []}]) +def test_malformed_json_structure_or_missing_calculation_section_fails(payload) -> None: + case = compare_pyscf.CASES[0] + result = report.validate_case(case, payload, (-75.0, None, None), [], True) + assert result["status"] == "FAIL" + assert result["errors"] + + +@pytest.mark.parametrize( + ("case_name", "rustiq_change", "pyscf", "converged", "diagnostic"), + [ + ("h2-sto-3g-rhf", lambda data: data["calculation"]["hf"].update(total_energy=float("inf")), (-75.0, None, None), True, "HF energy"), + ("h2-sto-3g-rhf-mp2", lambda data: data["calculation"]["mp2"].update(correlation_energy=float("nan")), (-75.0, -0.2, None), True, "MP2 correlation"), + ("oh-sto-3g-uhf", lambda data: data["calculation"]["hf"].update(spin={"s_squared": float("nan")}), (-75.0, None, 0.75), True, "Spin expectation"), + ], +) +def test_non_finite_scientific_values_fail(case_name, rustiq_change, pyscf, converged, diagnostic) -> None: + case = next(item for item in compare_pyscf.CASES if item.name == case_name) + rustiq = valid_rustiq(case) + rustiq_change(rustiq) + errors: list[str] = [] + result = report.validate_case(case, rustiq, pyscf, errors, converged) + assert result["status"] == "FAIL" + assert any(diagnostic in error for error in errors) + + +def test_recorded_execution_error_forces_failure() -> None: + case = compare_pyscf.CASES[0] + errors = ["RustiQ calculation failed: synthetic diagnostic"] + result = report.validate_case(case, valid_rustiq(case), (-75.0, None, None), errors, True) + assert result["hf_pass"] + assert result["status"] == "FAIL" + + +def test_combined_uhf_mp2_spin_validation_and_failed_case_serialization() -> None: + case = next(item for item in compare_pyscf.CASES if item.name == "oh-sto-3g-uhf-mp2") + result = report.validate_case(case, valid_rustiq(case), (-75.0, -0.2, 0.75), [], True) + assert result["hf_pass"] and result["mp2_pass"] and result["spin_pass"] + assert result["status"] == "PASS" + failed = report.validate_case(case, valid_rustiq(case), (-75.0, -0.2, 0.75), ["recorded execution error"], True) + payload = json.loads(json.dumps({"cases": [failed]}, allow_nan=False)) + assert payload["cases"][0]["status"] == "FAIL" + markdown = report.render_markdown({**payload, "metadata": {}}) + assert "recorded execution error" in markdown + assert "0/1 cases passed" in markdown + + +def test_json_serialization_preserves_full_precision_and_rejects_nan(tmp_path: Path) -> None: + payload = {"energy": -1.2345678901234567} + path = tmp_path / "report.json" + path.write_text(json.dumps(payload, allow_nan=False), encoding="utf-8") + assert json.loads(path.read_text())["energy"] == payload["energy"] + + +def test_markdown_scientific_notation_and_failure_detail() -> None: + item = {"case_id": "failed", "molecule": "H2", "geometry_path": "samples/h2/molecule.xyz", + "basis": "sto-3g", "charge": 0, "spin_n_alpha_minus_n_beta": 0, + "ao_basis_dimension": 2, "hf_method": "RHF", "rustiq_hf_total_energy_hartree": None, + "pyscf_hf_total_energy_hartree": None, "hf_absolute_difference_hartree": None, + "hf_tolerance_hartree": 1e-10, "hf_pass": False, "status": "FAIL", + "errors": ["HF energy comparison failed"]} + md = report.render_markdown({"cases": [item], "metadata": {}}) + assert "0/1 cases passed" in md + assert "HF energy comparison failed" in md + assert "—" in md + + +def test_molecular_formula_uses_geometry_atoms() -> None: + case = next(c for c in compare_pyscf.CASES if c.name == "h2o-sto-3g-rhf") + assert report.molecular_formula(case.xyz) == "H2O" + h2_case = compare_pyscf.CASES[0] + assert report.molecular_formula(h2_case.xyz) == "H2"