diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 7732e2c..dd61a4b 100644 --- a/ThermoScreening/thermo/api.py +++ b/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 diff --git a/ThermoScreening/thermo/screening.py b/ThermoScreening/thermo/screening.py index 3ea444c..be5df6f 100644 --- a/ThermoScreening/thermo/screening.py +++ b/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 ``). + must be installed (``thermo setup-dftb --solvent ``). 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 diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py index 3565da0..af3a1a5 100644 --- a/ThermoScreening/thermo/thermo.py +++ b/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)": diff --git a/docs/benchmarks/entropy_accuracy.md b/docs/benchmarks/entropy_accuracy.md new file mode 100644 index 0000000..9e03ccd --- /dev/null +++ b/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. diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index e51289d..72833bc 100644 --- a/tests/calculator/test_dftbplus.py +++ b/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 ).", + ) + 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 diff --git a/tests/calculator/test_xtb.py b/tests/calculator/test_xtb.py index d186772..8a7907e 100644 --- a/tests/calculator/test_xtb.py +++ b/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 diff --git a/tests/cli/test_dftb_setup.py b/tests/cli/test_dftb_setup.py index be27468..19b26cf 100644 --- a/tests/cli/test_dftb_setup.py +++ b/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) diff --git a/tests/thermo/test_heat_capacity_cv.py b/tests/thermo/test_heat_capacity_cv.py new file mode 100644 index 0000000..2b7d76c --- /dev/null +++ b/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) diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py index 6ce591a..ade2b6b 100644 --- a/tests/thermo/test_thermo.py +++ b/tests/thermo/test_thermo.py @@ -130,15 +130,16 @@ def test_thermo_rejects_negative_pressure(): _CM_TO_EV = 1.23984198e-4 _EVK_TO_CALMOLK = 1.602176634e-19 * 6.02214076e23 / 4.184 +_EV_TO_KCAL = 1.602176634e-19 * 6.02214076e23 / 4184.0 _T, _P = 298.15, 101325.0 _R_CALMOLK = 8.314462618 / 4.184 -def _ts_thermo(symbols, positions, real_freqs, dof, quasi_rrho=False): +def _ts_thermo(symbols, positions, real_freqs, dof, quasi_rrho=False, charge=0, spin=None): atoms = [Atom(symbol=s, position=np.array(p, float)) for s, p in zip(symbols, positions)] pad = np.concatenate([np.zeros(3 * len(atoms) - dof), np.asarray(real_freqs, float)]) system = System( - atoms, periodicity=False, cell=None, charge=0, + atoms, periodicity=False, cell=None, charge=charge, spin=spin, electronic_energy=0.0, vibrational_frequencies=pad, ) thermo = Thermo( @@ -179,6 +180,125 @@ def test_total_entropy_matches_ase(symbols, positions, freqs, dof, geometry, sig assert ts_total == pytest.approx(ase_total, abs=0.05) +from ase.build import molecule # noqa: E402 + + +def _nh4_cation(): + s = 1.03 / np.sqrt(3.0) + return Atoms( + "NH4", + positions=[[0, 0, 0], [s, s, s], [s, -s, -s], [-s, s, -s], [-s, -s, s]], + ) + + +# name/factory, freqs (cm^-1), dof, ASE geometry, symmetry number, charge +_SHG_CASES = [ + ("H2O", lambda: molecule("H2O"), [1595, 3657, 3756], 3, "nonlinear", 2, 0), + ("CO2", lambda: molecule("CO2"), [667, 667, 1333, 2349], 4, "linear", 2, 0), + ("N2", lambda: molecule("N2"), [2359], 1, "linear", 2, 0), + ("NH3", lambda: molecule("NH3"), [950, 1627, 1627, 3337, 3444, 3444], 6, "nonlinear", 3, 0), + ("CH4", lambda: molecule("CH4"), + [1306, 1306, 1306, 1534, 1534, 2917, 3019, 3019, 3019], 9, "nonlinear", 12, 0), + ("CO", lambda: molecule("CO"), [2143], 1, "linear", 1, 0), + ("HCN", lambda: molecule("HCN"), [712, 712, 2089, 3311], 4, "linear", 1, 0), + ("C2H2", lambda: molecule("C2H2"), + [612, 612, 730, 730, 1974, 3289, 3374], 7, "linear", 2, 0), + ("C2H6", lambda: molecule("C2H6"), + [289, 822, 822, 995, 1190, 1190, 1379, 1388, 1468, 1468, 1469, 1469, + 2896, 2954, 2969, 2969, 2985, 2985], 18, "nonlinear", 6, 0), + ("CH3OH", lambda: molecule("CH3OH"), + [295, 1033, 1060, 1165, 1345, 1455, 1477, 1477, 2844, 2960, 3000, 3681], + 12, "nonlinear", 1, 0), + ("Ar", lambda: Atoms("Ar", positions=[[0, 0, 0]]), [], 0, "monatomic", 1, 0), + # charged species exercise the charge path (10 electrons each -> spin 0, so + # thermochemically identical to their neutral RRHO) + ("OH-", lambda: Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]), + [3700], 1, "linear", 1, -1), + ("NH4+", _nh4_cation, + [1400, 1400, 1400, 1680, 1680, 3040, 3145, 3145, 3145], 9, "nonlinear", 12, 1), +] + + +@pytest.mark.parametrize( + "name,factory,freqs,dof,geometry,sigma,charge", _SHG_CASES, + ids=[c[0] for c in _SHG_CASES], +) +def test_S_H_G_match_ase(name, factory, freqs, dof, geometry, sigma, charge): + atoms = factory() + symbols = list(atoms.get_chemical_symbols()) + positions = atoms.get_positions() + + thermo = _ts_thermo(symbols, positions, freqs, dof, charge=charge) + + # ASE reference with masses pinned to the tool's atomic masses + ase_atoms = Atoms(symbols=symbols, positions=positions) + ase_atoms.set_masses( + [Atom(symbol=s, position=np.zeros(3)).mass for s in symbols] + ) + igt = IdealGasThermo( + vib_energies=[f * _CM_TO_EV for f in freqs], + geometry=geometry, atoms=ase_atoms, + symmetrynumber=sigma, spin=0, potentialenergy=0.0, + ) + s_ase = igt.get_entropy(_T, _P, verbose=False) * _EVK_TO_CALMOLK + h_ase = igt.get_enthalpy(_T, verbose=False) * _EV_TO_KCAL + g_ase = igt.get_gibbs_energy(_T, _P, verbose=False) * _EV_TO_KCAL + + assert thermo._system.rotational_symmetry_number == sigma + assert thermo.total_entropy("cal/(mol*K)") == pytest.approx(s_ase, abs=2e-3) + assert thermo.total_enthalpy("kcal") == pytest.approx(h_ase, abs=2e-3) + assert thermo.total_gibbs_free_energy("kcal") == pytest.approx(g_ase, abs=2e-3) + + +def _grimme_vib_entropy(freqs_cm, quasi): + """Independent Grimme quasi-RRHO / harmonic vibrational entropy (cal/mol/K).""" + h, kB, c, R, cal = 6.62607015e-34, 1.380649e-23, 299792458.0, 8.314462618, 4.184 + Bav, nu0 = 1.0e-44, 100.0 + total = 0.0 + for nu_cm in freqs_cm: + nu_hz = nu_cm * c * 100.0 + x = h * nu_hz / (kB * _T) + s_ho = x / (np.exp(x) - 1) - np.log(1 - np.exp(-x)) + if not quasi: + total += s_ho + continue + mu = h / (8 * np.pi**2 * nu_hz) + mu_eff = mu * Bav / (mu + Bav) + s_fr = 0.5 + np.log(np.sqrt(8 * np.pi**3 * mu_eff * kB * _T / h**2)) + w = 1.0 / (1.0 + (nu0 / nu_cm) ** 4) + total += w * s_ho + (1 - w) * s_fr + return R * total / cal + + +def test_quasi_rrho_matches_independent_grimme_and_is_invariant(): + # ethane geometry with genuine low modes (25/40/90 cm^-1) + symbols = ["C", "C", "H", "H", "H", "H", "H", "H"] + positions = [[0, 0, 0], [1.5, 0, 0], [-0.4, 1.0, 0], [-0.4, -0.5, 0.87], + [-0.4, -0.5, -0.87], [1.9, 1.0, 0], [1.9, -0.5, 0.87], [1.9, -0.5, -0.87]] + freqs = [25.0, 40.0, 90.0, 300.0, 820.0, 995.0, 995.0, 1206.0, 1206.0, + 1388.0, 1469.0, 1479.0, 1486.0, 2896.0, 2915.0, 2954.0, 2969.0, 2985.0] + + harmonic = _ts_thermo(symbols, positions, freqs, 18) + qrrho = _ts_thermo(symbols, positions, freqs, 18, quasi_rrho=True) + modes = harmonic._system.real_vibrational_frequencies + + s_ho = harmonic.total_entropy("cal/(mol*K)") + s_q = qrrho.total_entropy("cal/(mol*K)") + + # quasi-RRHO lowers the entropy (low modes tamed) + assert s_q < s_ho - 2.0 + # non-circular: the tool's HO->qRRHO shift equals the independent Grimme shift + tool_shift = s_q - s_ho + indep_shift = _grimme_vib_entropy(modes, quasi=True) - _grimme_vib_entropy(modes, quasi=False) + assert tool_shift == pytest.approx(indep_shift, abs=0.02) + + # quasi_rrho changes ONLY the entropy: enthalpy and Cv are invariant + 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 + ) + + def test_quasi_rrho_matches_harmonic_for_high_frequencies(): # water: all modes are high (>1500 cm^-1) -> weight ~ 1 -> qRRHO == harmonic args = (["O", "H", "H"],