Close validation gaps: broaden cross-checks, benchmark, and document

ThermoScreening/thermo/api.py

@@ -567,6 +567,17 @@ def dftbplus_thermo(
     -------
     Thermo
         The thermo calculation object.
+
+    Notes
+    -----
+    Reproducibility: identical inputs on the same machine give bit-identical
+    results (DFTB+ is thread-invariant). Across machines the optimiser stopping
+    point can vary by ~3e-5 A, giving ~2e-9 Hartree / ~3e-4 cal/mol/K scatter --
+    physically negligible (~1e-6 of the method error).
+
+    Solvation accuracy: with 3ob/mio the GBSA parameters are GFN-xTB-fit, so
+    solvation free energies are only qualitative and for small neutral solutes
+    can be positive/non-monotonic in the solvent dielectric.
     """
 
     # Resolve the spin (electron-count guess when not given) so the calculation
@@ -662,6 +673,13 @@ def xtb_thermo(
     -------
     Thermo
         The thermo calculation object.
+
+    Notes
+    -----
+    Reproducibility: with ``OMP_NUM_THREADS=1`` identical inputs give
+    bit-identical energies on the same machine; the entropy is thread-invariant
+    regardless. The optimiser stopping point is geometry/tolerance sensitive, so
+    tighten ``fmax`` (e.g. 1e-3) for more environment-stable geometries.
     """
 
     # Resolve the spin the same way as the DFTB+ path so the calculation and the

ThermoScreening/thermo/screening.py

@@ -177,7 +177,10 @@ def screen(
     solvent : str, optional
         Solvent name for GBSA/ALPB implicit solvation applied to every molecule
         (e.g. ``"water"``). Defaults to gas phase. The solvent parameter file
-        must be installed (``thermo setup-dftb --solvent <name>``).
+        must be installed (``thermo setup-dftb --solvent <name>``). With the
+        DFTB+ engine the GBSA parameters are GFN-xTB-fit, so solvation free
+        energies are only qualitative (and can be non-monotonic in the dielectric
+        for small neutral solutes).
     quasi_rrho : bool
         If True, use Grimme's quasi-RRHO treatment for the vibrational entropy
         (recommended for flexible molecules with low-frequency modes). Default

ThermoScreening/thermo/thermo.py

@@ -33,7 +33,7 @@ class Thermo:
     temperature : float
         The temperature of the system in Kelvin.
     pressure : float
-        The pressure of the system in bar.
+        The pressure of the system in Pa.
     system : System
         The system information of the system.
     engine : str
@@ -71,7 +71,7 @@ def __init__(
         temperature : float
             The temperature of the system in Kelvin.
         pressure : float
-            The pressure of the system in bar.
+            The pressure of the system in Pa.
         system : System
             The system information of the system.
         engine : str
@@ -882,7 +882,12 @@ def total_gibbs_free_energy(self, unit: str) -> float:
 
     def total_heat_capacity(self, unit: str) -> float:
         """
-        Returns the total heat capacity of the system based on the unit.
+        Returns the constant-volume (ideal-gas) heat capacity Cv of the system.
+
+        Cv = translational (3/2 R) + rotational (n_rot/2 R) + vibrational +
+        electronic (0). For the constant-pressure heat capacity of an ideal gas
+        add R: Cp = Cv + R. (Cross-checked against native xTB, whose reported
+        heat capacity is Cp: tool Cv == xtb Cp - R.)
 
         Parameters
         ----------
@@ -897,7 +902,7 @@ def total_heat_capacity(self, unit: str) -> float:
         Returns
         -------
         float
-            The total heat capacity of the system.
+            The total constant-volume heat capacity Cv of the system.
         """
 
         if unit == "cal/(mol*K)":

docs/benchmarks/entropy_accuracy.md

@@ -0,0 +1,45 @@
+# 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.

tests/calculator/test_dftbplus.py

@@ -9,8 +9,13 @@
 import ThermoScreening.calculator.dftbplus as dftbplus_module
 from ThermoScreening.calculator import Geoopt, Hessian, Modes
 from ThermoScreening.calculator.dftbplus import _slako_dir, dftb_3ob_parameters
+from ThermoScreening.cli.dftb_setup import gbsa_param_path
 from ThermoScreening.thermo.api import dftbplus_thermo
 
+_solvent_params_ready = all(
+    gbsa_param_path(name).exists() for name in ("dmso", "thf", "water")
+)
+
 # --------------------------------------------------------------------------- #
 
 
@@ -582,3 +587,64 @@ def test_dftbplus_thermo(self, test_with_data_dir):
         assert np.allclose(thermo.total_EeGtot(),
                            -4.0595598803365744,
                            atol=1e-5)
+
+    @pytest.mark.parametrize("name,exp_S", [("H2O", 45.1), ("CH4", 44.5), ("N2", 45.8)])
+    def test_dftbplus_entropy_vs_experiment(self, name, exp_S, tmp_path):
+        # DFTB+ 3ob gas-phase entropy within ~2 cal/mol/K of experiment, and the
+        # rotational symmetry number is auto-detected correctly
+        from ase.build import molecule
+
+        thermo = dftbplus_thermo(
+            molecule(name), pressure=100000.0,
+            directory=str(tmp_path / name), **dftb_3ob_parameters,
+        )
+        assert thermo.total_entropy("cal/(mol*K)") == pytest.approx(exp_S, abs=2.0)
+        assert thermo._system.rotational_symmetry_number == {"H2O": 2, "CH4": 12, "N2": 2}[name]
+
+    def test_dftbplus_thermo_is_reproducible(self, tmp_path):
+        # identical inputs on the same machine give bit-identical DFTB+ results
+        def oh():
+            return Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]])
+
+        a = dftbplus_thermo(oh(), directory=str(tmp_path / "a"), **dftb_3ob_parameters)
+        b = dftbplus_thermo(oh(), directory=str(tmp_path / "b"), **dftb_3ob_parameters)
+
+        assert abs(a.electronic_energy() - b.electronic_energy()) < 1e-12
+        assert abs(a.total_entropy("cal/(mol*K)") - b.total_entropy("cal/(mol*K)")) < 1e-9
+        # loose, version-robust regression anchors
+        assert a.electronic_energy() == pytest.approx(-3.579083, abs=1e-5)
+        assert a.total_entropy("cal/(mol*K)") == pytest.approx(42.6520, abs=1e-3)
+
+    @pytest.mark.skipif(
+        not _solvent_params_ready,
+        reason="GBSA params for dmso/thf/water not installed "
+        "(run: thermo setup-dftb --solvent <name>).",
+    )
+    def test_dftbplus_solvation_across_solvents(self, tmp_path):
+        # several solvents each apply a GBSA block; dmso/thf stabilise the solute
+        # and solvation enters the Hessian (frequencies shift). Magnitudes and the
+        # solvent ordering are conformer-sensitive, so only robust facts are asserted.
+        from ase.build import molecule
+
+        def run(name, solvent=None):
+            return dftbplus_thermo(
+                molecule("CH3OH"), directory=str(tmp_path / name),
+                solvent=solvent, **dftb_3ob_parameters,
+            )
+
+        gas = run("gas")
+        assert "GeneralizedBorn" not in (tmp_path / "gas" / "dftb_in.hsd").read_text()
+        e_gas = gas.electronic_energy()
+
+        solvated = {s: run(s, solvent=s) for s in ("dmso", "thf", "water")}
+        for s in solvated:
+            assert "GeneralizedBorn" in (tmp_path / s / "dftb_in.hsd").read_text()
+
+        # mid/polar solvents stabilise the solute (sign robust; magnitude is not)
+        assert solvated["dmso"].electronic_energy() - e_gas < -1e-3
+        assert solvated["thf"].electronic_energy() - e_gas < -1e-3
+
+        # solvation enters the Hessian: the lowest mode shifts by several cm^-1
+        gas_low = np.sort(gas._system.real_vibrational_frequencies)[0]
+        dmso_low = np.sort(solvated["dmso"]._system.real_vibrational_frequencies)[0]
+        assert abs(dmso_low - gas_low) > 5.0

tests/calculator/test_xtb.py

@@ -44,17 +44,88 @@ def test_optimise_and_frequencies_with_emt(monkeypatch, tmp_path):
     assert frequencies[-1] > 0
 
 
+def test_xtb_thermo_has_no_solvation_parameter():
+    # regression guard: xTB-side implicit solvation is intentionally not wired up
+    import inspect
+    from ThermoScreening.thermo.api import xtb_thermo
+
+    assert "solvent" not in inspect.signature(xtb_thermo).parameters
+
+
 tblite_available = importlib.util.find_spec("tblite") is not None
+xtb_skip = pytest.mark.skipif(not tblite_available, reason="tblite (GFN-xTB) is not installed.")
 
 
-@pytest.mark.skipif(not tblite_available, reason="tblite (GFN-xTB) is not installed.")
-def test_xtb_thermo_runs_real_gfn2(tmp_path):
-    # real GFN2-xTB end-to-end: gas-phase water entropy is close to experiment
+@xtb_skip
+@pytest.mark.parametrize("name,exp_S", [("H2O", 45.1), ("CH4", 44.5), ("N2", 45.8)])
+def test_xtb_thermo_entropy_vs_experiment(name, exp_S, tmp_path):
+    # real GFN2-xTB gas-phase entropy is within ~2 cal/mol/K of experiment
     from ase.build import molecule
     from ThermoScreening.thermo.api import xtb_thermo
 
-    thermo = xtb_thermo(molecule("H2O"), directory=str(tmp_path / "w"))
-    entropy = thermo.total_entropy("cal/(mol*K)")
+    thermo = xtb_thermo(molecule(name), pressure=100000.0, directory=str(tmp_path / name))
+    assert thermo.total_entropy("cal/(mol*K)") == pytest.approx(exp_S, abs=2.0)
+
+
+# n-octane geometry (embedded so the test needs only tblite, not a 3D builder)
+_OCTANE = [
+    ("C", (4.1592, -0.4455, -0.0648)), ("C", (2.7620, 0.1265, 0.1227)),
+    ("C", (1.7204, -0.9840, 0.2660)), ("C", (0.3124, -0.4653, 0.5771)),
+    ("C", (-0.3124, 0.3231, -0.5771)), ("C", (-1.7204, 0.8418, -0.2660)),
+    ("C", (-2.7620, -0.2686, -0.1227)), ("C", (-4.1592, 0.3033, 0.0648)),
+    ("H", (4.8889, 0.3636, -0.1690)), ("H", (4.4517, -1.0575, 0.7944)),
+    ("H", (4.2106, -1.0677, -0.9640)), ("H", (2.5272, 0.7621, -0.7373)),
+    ("H", (2.7520, 0.7646, 1.0138)), ("H", (1.6959, -1.5906, -0.6475)),
+    ("H", (2.0238, -1.6535, 1.0806)), ("H", (0.3407, 0.1568, 1.4797)),
+    ("H", (-0.3149, -1.3331, 0.8095)), ("H", (-0.3407, -0.2989, -1.4797)),
+    ("H", (0.3149, 1.1909, -0.8095)), ("H", (-1.6959, 1.4484, 0.6475)),
+    ("H", (-2.0238, 1.5113, -1.0806)), ("H", (-2.5272, -0.9043, 0.7374)),
+    ("H", (-2.7520, -0.9068, -1.0137)), ("H", (-4.8889, -0.5058, 0.1691)),
+    ("H", (-4.4517, 0.9153, -0.7944)), ("H", (-4.2106, 0.9255, 0.9640)),
+]
+
+
+@xtb_skip
+def test_xtb_quasi_rrho_on_real_floppy_molecule(tmp_path):
+    # a real molecule with engine-generated sub-100 cm^-1 modes; quasi-RRHO tames
+    # them, and (optimising once) leaves enthalpy and Cv invariant.
+    import numpy as np
+    from ase import Atoms
+    from ThermoScreening.calculator.xtb import optimise_and_frequencies, xtb_calculator
+    from ThermoScreening.thermo.api import run_thermo
+
+    octane = Atoms([s for s, _ in _OCTANE], positions=[p for _, p in _OCTANE])
+    octane.info["charge"] = 0
+    octane.info["spin"] = 0
+    opt, energy, freqs = optimise_and_frequencies(
+        octane.copy(), xtb_calculator("GFN2-xTB"), fmax=0.02
+    )
+
+    kept = np.sort(freqs)[-(3 * len(opt) - 6):]
+    assert np.any(kept < 100.0)  # genuine engine-generated low-frequency modes
+
+    common = dict(atoms=opt, energy=energy, engine="xtb", pressure=101325)
+    harmonic = run_thermo(freqs, quasi_rrho=False, **common)
+    qrrho = run_thermo(freqs, quasi_rrho=True, **common)
+
+    assert qrrho.total_entropy("cal/(mol*K)") < harmonic.total_entropy("cal/(mol*K)") - 2.0
+    # invariance (same optimised freqs -> only entropy differs)
+    assert qrrho.total_enthalpy("H") == pytest.approx(harmonic.total_enthalpy("H"), abs=1e-12)
+    assert qrrho.total_heat_capacity("cal/(mol*K)") == pytest.approx(
+        harmonic.total_heat_capacity("cal/(mol*K)"), abs=1e-12
+    )
+
+
+@xtb_skip
+def test_xtb_thermo_is_reproducible_single_thread(monkeypatch, tmp_path):
+    # same inputs + OMP_NUM_THREADS=1 -> bit-identical energy, thread-invariant S
+    from ase import Atoms
+    from ThermoScreening.thermo.api import xtb_thermo
+
+    monkeypatch.setenv("OMP_NUM_THREADS", "1")
+    oh = lambda: Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]])
+    a = xtb_thermo(oh(), directory=str(tmp_path / "a"), fmax=1e-3)
+    b = xtb_thermo(oh(), directory=str(tmp_path / "b"), fmax=1e-3)
 
-    # experimental standard molar entropy of gaseous water ~ 45.1 cal/mol/K
-    assert entropy == pytest.approx(45.1, abs=2.0)
+    assert abs(a.electronic_energy() - b.electronic_energy()) < 1e-9
+    assert abs(a.total_entropy("cal/(mol*K)") - b.total_entropy("cal/(mol*K)")) < 1e-6

tests/cli/test_dftb_setup.py

@@ -157,6 +157,26 @@ def test_install_gbsa_param_downloads_file(tmp_path):
     assert installed.read_text(encoding="utf-8") == "solvent parameters"
 
 
+@pytest.mark.parametrize(
+    "solvent,stem",
+    [("water", "h2o"), ("dmso", "dmso"), ("acetonitrile", "acetonitrile"), ("thf", "thf")],
+)
+def test_gbsa_plumbing_for_multiple_solvents(monkeypatch, tmp_path, solvent, stem):
+    assert _solvent_stem(solvent) == stem
+
+    path = gbsa_param_path(solvent, install_root=tmp_path)
+    assert path.name == f"param_gbsa_{stem}.txt"
+    assert path.parent.name == "gfn2-0-1"
+
+    # a mocked download lands the file at the resolved path
+    def fake_download(url, destination):
+        destination.write_text("params", encoding="utf-8")
+
+    monkeypatch.setattr(dftb_setup, "_download_file", fake_download)
+    installed = install_gbsa_param(solvent, install_root=tmp_path)
+    assert installed == path.resolve()
+
+
 def test_install_gbsa_param_requires_downloaded_file(monkeypatch, tmp_path):
     # a "download" that writes nothing must not silently succeed
     monkeypatch.setattr(dftb_setup, "_download_file", lambda url, destination: None)

tests/thermo/test_heat_capacity_cv.py

@@ -0,0 +1,54 @@
+"""The reported heat capacity is the ideal-gas constant-volume Cv.
+
+Cross-checked non-circularly against the closed-form ideal-gas Cv (re-derived
+here from the partition functions) and pinned to reference values. ASE's
+IdealGasThermo exposes no Cv accessor, so this is validated analytically rather
+than against ASE.
+"""
+
+import numpy as np
+import pytest
+from ase.build import molecule
+
+from ThermoScreening.thermo.api import run_thermo
+
+_h = 6.62607015e-34
+_c = 299792458.0
+_kB = 1.380649e-23
+_R = 8.314462618 / 4.184  # cal/(mol*K)
+_T = 298.15
+
+
+def _analytic_cv(wavenumbers, n_rot):
+    """Closed-form ideal-gas Cv (cal/mol/K) from translation, rotation, vibration."""
+    x = _h * _c * np.asarray(wavenumbers, float) * 100.0 / (_kB * _T)
+    vib = _R * np.sum(x**2 * np.exp(x) / np.expm1(x) ** 2)
+    return 1.5 * _R + 0.5 * n_rot * _R + vib
+
+
+# molecule, vibrational wavenumbers (cm^-1), rotational dof, reference Cv
+_CASES = [
+    ("H2O", [1538.86, 3642.96, 3651.45], 3, 6.026981),
+    ("CH4", [1384.87, 1384.87, 1384.87, 1556.89, 1556.89,
+             3096.01, 3109.71, 3109.71, 3109.71], 3, 6.418981),
+    ("CO2", [599.92, 599.92, 1425.30, 2594.02], 2, 7.130121),
+    ("N2", [2437.38], 2, 4.970154),
+]
+
+
+@pytest.mark.parametrize("name,wavenumbers,n_rot,expected_cv", _CASES)
+def test_heat_capacity_is_constant_volume_cv(name, wavenumbers, n_rot, expected_cv):
+    thermo = run_thermo(
+        np.array(wavenumbers),
+        atoms=molecule(name),
+        temperature=_T,
+        pressure=101325,
+        energy=0.0,
+        engine="dftb+",
+    )
+    cv = thermo.total_heat_capacity("cal/(mol*K)")
+
+    assert thermo._n_rot == n_rot
+    # matches the independently re-derived analytic Cv, and the pinned reference
+    assert cv == pytest.approx(_analytic_cv(wavenumbers, n_rot), abs=1e-3)
+    assert cv == pytest.approx(expected_cv, abs=1e-3)