diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 2a6ebb8..09c6f39 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -25,15 +25,27 @@ jobs: uses: actions/setup-python@v7 with: python-version: "3.12" + - name: Cache DFTB+ parameters + uses: actions/cache@v6 + with: + path: ~/.local/share/thermoscreening/slakos/3ob-3-1 + key: dftbparams-3ob-3-1-v1 - name: Add conda to system path run: | # $CONDA is an environment variable pointing to the root of the miniconda directory echo $CONDA/bin >> $GITHUB_PATH - name: Install dependencies run: | - conda install -y -c conda-forge python=3.12 gsl dftbplus + conda install -y -c conda-forge python=3.12 gsl 'dftbplus=25.1=nompi_*' python -m pip install --upgrade pip python -m pip install '.[test,lint]' + - name: Verify DFTB+ workflow + env: + THERMOSCREENING_REQUIRE_DFTB: "1" + run: | + thermo setup-dftb --parameter-set 3ob + thermo doctor --engine dftb+ + python -m pytest tests/integration/test_dftb_workflow.py -q - name: Lint with Pylint run: python -m pylint ThermoScreening - name: Test with pytest diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py index 0b7db23..153174a 100644 --- a/ThermoScreening/calculator/dftbplus.py +++ b/ThermoScreening/calculator/dftbplus.py @@ -582,9 +582,6 @@ def read(self): Hamiltonian_Filling="Fermi {", Hamiltonian_Filling_empty="Temperature [Kelvin] = 300", - # Convergence helper - Hamiltonian_Mixer="DIIS{}", - # Are guessed by ase Hamiltonian_MaxAngularMomentum_="", @@ -632,9 +629,6 @@ def read(self): Hamiltonian_Filling="Fermi {", Hamiltonian_Filling_empty="Temperature [Kelvin] = 300", - # Convergence helper - Hamiltonian_Mixer="DIIS{}", - # Are guessed by ase Hamiltonian_MaxAngularMomentum_="", diff --git a/docs/benchmarks/anthraquinone_workflow.rst b/docs/benchmarks/anthraquinone_workflow.rst new file mode 100644 index 0000000..beab1ec --- /dev/null +++ b/docs/benchmarks/anthraquinone_workflow.rst @@ -0,0 +1,17 @@ +Anthraquinone workflow regression +================================= + +The test suite contains a compact three-molecule anthraquinone fixture covering +the parent compound and two hydroxy substitution positions. It validates: + +- deterministic SMILES embedding and molecular formulas; +- construction of the oxidized, singly reduced and doubly reduced states; +- separate reference calibration of both reduction steps; +- potential-inversion classification; and +- reproducible input-set provenance. + +The fixture tests workflow invariants with controlled state energies. It is not +an absolute-potential accuracy benchmark and does not copy results from an +external dataset. Quantitative comparisons require the same structures, +electronic-structure method, solvent model, reference electrode and standard +states. diff --git a/docs/benchmarks/entropy_accuracy.md b/docs/benchmarks/entropy_accuracy.md deleted file mode 100644 index 9e03ccd..0000000 --- a/docs/benchmarks/entropy_accuracy.md +++ /dev/null @@ -1,45 +0,0 @@ -# Gas-phase entropy accuracy - -Standard molar entropies S°(298.15 K, 1 bar) from ThermoScreening compared with -experiment, for both engines. These are *method* accuracies (GFN2-xTB and -DFTB3/3ob), not tool errors — the tool's thermochemistry itself matches ASE -`IdealGasThermo` to < 0.0003 cal/mol/K (a non-circular cross-check). - -| Molecule | σ | xTB (GFN2) | DFTB+ 3ob | Experiment | -|----------|---|-----------|-----------|------------| -| H₂O | 2 | 45.08 | 45.07 | 45.1 | -| CH₄ | 12| 44.44 | 44.47 | 44.5 | -| NH₃ | 3 | 45.96 | 46.00 | 46.0 | -| N₂ | 2 | 45.79 | 45.76 | 45.8 | -| CO₂ | 2 | 51.23 | 51.53 | 51.1 | -| CO | 1 | 47.23 | 47.14 | 47.2 | -| CH₃OH | 1 | 56.63 | 57.00 | 57.3 | -| C₂H₆ | 6 | 54.38 | 54.54 | 54.8 | - -(cal/mol/K) - -**Deviation from experiment** - -| Engine | MAD | RMSD | max\|dev\| | mean signed | -|------------|------|------|----------------|-------------| -| xTB (GFN2) | 0.17 | 0.28 | 0.67 (CH₃OH) | −0.13 | -| DFTB+ 3ob | 0.14 | 0.21 | 0.43 (CO₂) | −0.04 | - -Both engines reproduce experimental gas-phase entropies to well under -1 cal/mol/K. - -## Notes - -- **Symmetry numbers** are auto-detected correctly for all eight molecules - (H₂O = 2, CH₄ = 12, C₂H₆ = 6, …). -- **Standard-state convention.** The values above use `pressure=100000` (1 bar, - the S° convention). The tool's default `pressure=101325` (1 atm) raises the - entropy by `R·ln(101325/100000) = 0.026 cal/mol/K` — pass `pressure=100000` - to compare with tabulated S°. -- **Residual error** on CH₃OH and C₂H₆ comes from treating low-frequency - torsions as harmonic oscillators; `quasi_rrho=True` softens this. -- These numbers require the xtb/DFTB+ engines to regenerate optimized - geometries and frequencies, so they live as skippable integration tests - (`tests/calculator/test_xtb.py::test_xtb_thermo_entropy_vs_experiment`, - `tests/calculator/test_dftbplus.py::TestDftbplus::test_dftbplus_entropy_vs_experiment`) - rather than in the CI-safe suite. diff --git a/docs/benchmarks/entropy_accuracy.rst b/docs/benchmarks/entropy_accuracy.rst new file mode 100644 index 0000000..4c33e3c --- /dev/null +++ b/docs/benchmarks/entropy_accuracy.rst @@ -0,0 +1,94 @@ +Gas-phase entropy accuracy +========================== + +Standard molar entropies S°(298.15 K, 1 bar) from ThermoScreening compared with +experiment, for both engines. These are method accuracies (GFN2-xTB and +DFTB3/3ob), not tool errors. The thermochemistry implementation itself matches +ASE ``IdealGasThermo`` to within 0.0003 cal/(mol K). + +.. list-table:: + :header-rows: 1 + + * - Molecule + - Symmetry number + - xTB (GFN2) + - DFTB+ 3ob + - Experiment + * - H₂O + - 2 + - 45.08 + - 45.07 + - 45.1 + * - CH₄ + - 12 + - 44.44 + - 44.47 + - 44.5 + * - NH₃ + - 3 + - 45.96 + - 46.00 + - 46.0 + * - N₂ + - 2 + - 45.79 + - 45.76 + - 45.8 + * - CO₂ + - 2 + - 51.23 + - 51.53 + - 51.1 + * - CO + - 1 + - 47.23 + - 47.14 + - 47.2 + * - CH₃OH + - 1 + - 56.63 + - 57.00 + - 57.3 + * - C₂H₆ + - 6 + - 54.38 + - 54.54 + - 54.8 + +Values are in cal/(mol K). + +Deviation from experiment +------------------------- + +.. list-table:: + :header-rows: 1 + + * - Engine + - MAD + - RMSD + - Maximum absolute deviation + - Mean signed deviation + * - xTB (GFN2) + - 0.17 + - 0.28 + - 0.67 (CH₃OH) + - -0.13 + * - DFTB+ 3ob + - 0.14 + - 0.21 + - 0.43 (CO₂) + - -0.04 + +Both engines reproduce experimental gas-phase entropies to well under +1 cal/(mol K). + +Notes +----- + +- Symmetry numbers are detected correctly for all eight molecules. +- The values use ``pressure=100000`` (1 bar). The default pressure of 101325 Pa + raises the entropy by 0.026 cal/(mol K). +- Residual error for methanol and ethane comes from treating low-frequency + torsions as harmonic oscillators; ``quasi_rrho=True`` softens this. +- The Linux release gate runs a mandatory DFTB+ calculation. Broader xTB and + solvent checks remain optional integration tests. diff --git a/docs/index.rst b/docs/index.rst index 2397d7d..0ae2ead 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -107,6 +107,8 @@ Features usage configuration api + benchmarks/entropy_accuracy + benchmarks/anthraquinone_workflow Indices ------- diff --git a/pytest.ini b/pytest.ini index 5ee6477..2511dd3 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,2 +1,4 @@ [pytest] testpaths = tests +markers = + integration: requires an external calculation backend and parameter data diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index 6a66982..56fb655 100644 --- a/tests/calculator/test_dftbplus.py +++ b/tests/calculator/test_dftbplus.py @@ -485,6 +485,8 @@ def test_resolve_parameter_set_selects_hamiltonian_and_constants(): # 3ob is DFTB3 (third order), mio is DFTB2 (no third order) assert params_3ob["Hamiltonian_ThirdOrderFull"] == "Yes" assert "Hamiltonian_ThirdOrderFull" not in params_mio + assert "Hamiltonian_Mixer" not in params_3ob + assert "Hamiltonian_Mixer" not in params_mio assert spin_3ob is dftbplus_module.SPIN_CONSTANTS_3OB assert spin_mio is dftbplus_module.SPIN_CONSTANTS_MIO diff --git a/tests/data/regression/anthraquinones.csv b/tests/data/regression/anthraquinones.csv new file mode 100644 index 0000000..8a823c5 --- /dev/null +++ b/tests/data/regression/anthraquinones.csv @@ -0,0 +1,4 @@ +name,smiles,charge +anthraquinone,O=C1c2ccccc2C(=O)c2ccccc12,0 +1-hydroxyanthraquinone,O=C1c2ccccc2C(=O)c2cccc(O)c12,0 +2-hydroxyanthraquinone,O=C1c2ccccc2C(=O)c2ccc(O)cc12,0 diff --git a/tests/data/regression/water.xyz b/tests/data/regression/water.xyz new file mode 100644 index 0000000..ff7472c --- /dev/null +++ b/tests/data/regression/water.xyz @@ -0,0 +1,5 @@ +3 +DFTB+ release-gate geometry +O -0.00566600 0.40130200 0.00000000 +H 0.79261600 -0.19855100 0.00000000 +H -0.78695000 -0.20275100 0.00000000 diff --git a/tests/integration/test_dftb_workflow.py b/tests/integration/test_dftb_workflow.py new file mode 100644 index 0000000..158dfc6 --- /dev/null +++ b/tests/integration/test_dftb_workflow.py @@ -0,0 +1,77 @@ +"""Mandatory DFTB+ release-gate workflow.""" + +import json +import math +import os +from pathlib import Path +import subprocess +import sys + +import pytest + + +@pytest.mark.integration +def test_dftb_screen_end_to_end(tmp_path): + if os.getenv("THERMOSCREENING_REQUIRE_DFTB") != "1": + pytest.skip("set THERMOSCREENING_REQUIRE_DFTB=1 for the release gate") + + structure = ( + Path(__file__).resolve().parents[1] / "data" / "regression" / "water.xyz" + ) + manifest = tmp_path / "molecules.csv" + manifest.write_text( + f"name,path,charge\nwater,{structure},0\n", + encoding="utf-8", + ) + out = tmp_path / "results" + completed = subprocess.run( + [ + sys.executable, + "-m", + "ThermoScreening", + "screen", + str(manifest), + "--out", + str(out), + "--directory", + str(tmp_path / "calculations"), + "--engine", + "dftb+", + ], + capture_output=True, + text=True, + timeout=300, + check=False, + ) + + if completed.returncode: + calculation_dir = tmp_path / "calculations" / "water" + diagnostics = [completed.stdout, completed.stderr] + for filename in ("dftb_in.hsd", "second_derivative.out", "detailed.out"): + path = calculation_dir / filename + if path.is_file(): + diagnostics.append( + f"\n--- {filename} ---\n{path.read_text(encoding='utf-8')}" + ) + pytest.fail("".join(diagnostics)) + records = json.loads(out.with_suffix(".json").read_text(encoding="utf-8")) + assert len(records) == 1 + record = records[0] + assert record["name"] == "water" + assert record["formula"] == "H2O" + assert record["status"] == "ok" + for field in ( + "Eelec_hartree", + "G_total_hartree", + "S_cal_per_mol_K", + "Cv_cal_per_mol_K", + ): + assert math.isfinite(float(record[field])) + + metadata = json.loads( + out.with_name(f"{out.name}-run.json").read_text(encoding="utf-8") + ) + assert metadata["workflow"] == "thermochemistry_screen" + assert metadata["settings"]["engine"] == "dftb+" + assert metadata["input_set_fingerprint"] + assert metadata["jobs"][0]["structure_sha256"] diff --git a/tests/thermo/test_anthraquinone_regression.py b/tests/thermo/test_anthraquinone_regression.py new file mode 100644 index 0000000..b5147d1 --- /dev/null +++ b/tests/thermo/test_anthraquinone_regression.py @@ -0,0 +1,113 @@ +"""Scientific invariants for the anthraquinone redox workflow.""" + +import csv +import json +from pathlib import Path + +import ase.io +import pytest + +from ThermoScreening.thermo import screening +from ThermoScreening.thermo._units import HARTREE_TO_EV + + +_DATA = ( + Path(__file__).resolve().parents[1] + / "data" + / "regression" + / "anthraquinones.csv" +) +_ENERGIES = { + "anthraquinone": { + "oxidized": -100.0, + "reduced_once": -100.10, + "reduced_twice": -100.18, + }, + "1-hydroxyanthraquinone": { + "oxidized": -110.0, + "reduced_once": -110.115, + "reduced_twice": -110.205, + }, + "2-hydroxyanthraquinone": { + "oxidized": -120.0, + "reduced_once": -120.08, + "reduced_twice": -120.20, + }, +} + + +def _state_screen(source, **kwargs): + with open(source, newline="", encoding="utf-8") as handle: + rows = list(csv.DictReader(handle)) + records = [] + for row in rows: + molecule, state = row["name"].rsplit("--", 1) + records.append( + { + "name": row["name"], + "status": "ok", + "G_total_hartree": _ENERGIES[molecule][state], + } + ) + screening._write_results(records, kwargs["out"]) + return records + + +def _run_regression(monkeypatch, tmp_path, name): + monkeypatch.setattr(screening, "screen", _state_screen) + out = tmp_path / name + results = screening.redox_screen( + _DATA, + out=out, + directory=tmp_path / f"{name}-work", + reference="anthraquinone", + reference_e1=-0.75, + reference_e2=-1.40, + potential_scale="AQ reference", + max_conformers=3, + ) + metadata = json.loads( + out.with_name(f"{out.name}-run.json").read_text(encoding="utf-8") + ) + return results, metadata + + +def test_anthraquinone_workflow_regression(monkeypatch, tmp_path): + results, metadata = _run_regression(monkeypatch, tmp_path, "first") + by_name = {result["name"]: result for result in results} + + assert list(by_name) == [ + "anthraquinone", + "1-hydroxyanthraquinone", + "2-hydroxyanthraquinone", + ] + assert all(result["status"] == "ok" for result in results) + assert by_name["anthraquinone"]["E1_V"] == pytest.approx(-0.75) + assert by_name["anthraquinone"]["E2_V"] == pytest.approx(-1.40) + assert by_name["1-hydroxyanthraquinone"]["E1_V"] == pytest.approx( + -0.75 + 0.015 * HARTREE_TO_EV + ) + assert by_name["1-hydroxyanthraquinone"]["E2_V"] == pytest.approx( + -1.40 + 0.010 * HARTREE_TO_EV + ) + assert by_name["2-hydroxyanthraquinone"]["E1_V"] == pytest.approx( + -0.75 - 0.020 * HARTREE_TO_EV + ) + assert by_name["2-hydroxyanthraquinone"]["E2_V"] == pytest.approx( + -1.40 + 0.040 * HARTREE_TO_EV + ) + assert by_name["2-hydroxyanthraquinone"]["potential_inversion"] is True + + formulas = { + candidate["name"]: ase.io.read(candidate["path"]).get_chemical_formula() + for candidate in metadata["candidates"] + } + assert formulas == { + "anthraquinone": "C14H8O2", + "1-hydroxyanthraquinone": "C14H8O3", + "2-hydroxyanthraquinone": "C14H8O3", + } + assert metadata["input_set_fingerprint"] + + _results, repeated = _run_regression(monkeypatch, tmp_path, "repeated") + assert repeated["input_set_fingerprint"] == metadata["input_set_fingerprint"]