Add Grimme's quasi-RRHO vibrational entropy option (#67)

Adds a `quasi_rrho` option (CLI `--quasi-rrho`) for the vibrational
entropy.

## Why
The rigid-rotor-harmonic-oscillator entropy of a mode **diverges as its
frequency → 0**, so low-frequency modes (floppy torsions, weak/hindered
rotors, non-covalent complexes) get a spuriously large entropy — a
well-known RRHO weakness. This is also the dominant error left in
implicit-solvent free energies (#66).

## What
Grimme's quasi-RRHO ([Chem. Eur. J. 2012, 18,
9955](https://doi.org/10.1002/chem.201200497)): each mode's entropy is
interpolated between the harmonic-oscillator and free-rotor values,
weighted by `w = 1 / (1 + (100/ν)^4)`. High-frequency modes stay
harmonic; low-frequency modes approach the finite free-rotor limit.

- **`thermo/thermo.py`** — `Thermo(..., quasi_rrho=False)`;
`_compute_vibrational_entropy` uses the new `_quasi_rrho_entropy`
interpolation when enabled. The harmonic path is numerically unchanged,
so the **default preserves existing results**.
- **`thermo/api.py`, `thermo/screening.py`, `cli/thermo.py`** — thread
`quasi_rrho` through `run_thermo`, `dftbplus_thermo`, `screen`, and
`screen --quasi-rrho`.

Engine-independent (it acts on the frequencies), so it also applies to a
future xTB engine.

## Usage
```
thermo screen mols/ --solvent water --quasi-rrho
```

## Validation
- **Physics**: a molecule with 25/40 cm⁻¹ modes has its entropy reduced
by ~3.4 cal/mol/K (≈ +1 kcal/mol in G) — the spurious low-mode entropy
is tamed; a high-frequency-only case is unchanged (ΔS ≈ 0.01).
- **Real DFTB+**: ethane (lowest real vibration 277 cm⁻¹, all > 100)
gives qRRHO ≈ harmonic, and the harmonic S = 54.5 cal/mol/K matches the
experimental ~54.8 — confirming the correct mode set is used.
- Offline suite green (206 passed); real-DFTB+ integration green (54
passed, 0 skips); codecov patch fully covered.

Author: galjos

ThermoScreening/cli/thermo.py

@@ -106,6 +106,11 @@ def _command_parser():
         help="GBSA/ALPB implicit-solvation solvent (e.g. 'water') applied to every "
         "molecule. Default gas phase. Install with 'setup-dftb --solvent <name>'.",
     )
+    screen_parser.add_argument(
+        "--quasi-rrho", action="store_true",
+        help="Use Grimme's quasi-RRHO vibrational entropy (better for low-frequency "
+        "modes) instead of the pure harmonic oscillator.",
+    )
 
     return parser
 
@@ -187,6 +192,7 @@ def run_screen(parser_args):
         directory=parser_args.directory,
         parameter_set=parser_args.parameter_set,
         solvent=parser_args.solvent,
+        quasi_rrho=parser_args.quasi_rrho,
     )
 
     failed = sum(1 for record in results if record["status"] != "ok")

ThermoScreening/thermo/api.py

@@ -346,6 +346,7 @@ def run_thermo(
     charge=0.0,
     atoms=None,
     spin=None,
+    quasi_rrho=False,
 ):
     """
     Run the thermo calculation. Returns thermo calculation object.
@@ -369,6 +370,9 @@ def run_thermo(
         The system charge.
     atoms : ase.Atoms, optional
         An optimized geometry to use directly instead of reading ``coord_file``.
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO treatment for the vibrational entropy
+        instead of the pure harmonic oscillator. Default False.
 
     The coordinate file should be in xyz format.
 
@@ -411,7 +415,11 @@ def run_thermo(
     )
 
     thermo_setup = Thermo(
-        system=system_info, temperature=temperature, pressure=pressure, engine=engine
+        system=system_info,
+        temperature=temperature,
+        pressure=pressure,
+        engine=engine,
+        quasi_rrho=quasi_rrho,
     )
 
     thermo_setup.run()
@@ -507,6 +515,7 @@ def dftbplus_thermo(
         spin_constants=None,
         solvent=None,
         solvation_param_file=None,
+        quasi_rrho=False,
         **kwargs
     ):
     """
@@ -544,6 +553,9 @@ def dftbplus_thermo(
     solvation_param_file : str, optional
         Explicit path to a GBSA parameter file, overriding ``solvent`` (use a
         method-consistent set instead of the default GFN-fit one).
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO treatment for the vibrational entropy,
+        which tames the entropy of low-frequency modes. Default False.
 
     Other Parameters
     ----------------
@@ -596,6 +608,7 @@ def dftbplus_thermo(
             engine='dftb+',
             charge=charge,
             spin=spin,
+            quasi_rrho=quasi_rrho,
         )
 
     return thermo

ThermoScreening/thermo/screening.py

@@ -143,6 +143,7 @@ def screen(
     spin=None,
     parameter_set="3ob",
     solvent=None,
+    quasi_rrho=False,
 ):
     """
     Run a thermochemistry screen over a set of molecules.
@@ -175,6 +176,10 @@ def screen(
         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>``).
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO treatment for the vibrational entropy
+        (recommended for flexible molecules with low-frequency modes). Default
+        False (pure harmonic oscillator).
 
     Returns
     -------
@@ -207,6 +212,7 @@ def screen(
                 spin=job.spin,
                 spin_constants=spin_constants,
                 solvent=solvent,
+                quasi_rrho=quasi_rrho,
                 **parameters,
             )
             record.update(_thermo_summary(thermo))

ThermoScreening/thermo/thermo.py

@@ -49,8 +49,19 @@ class Thermo:
     logger = logging.getLogger(__package_name__).getChild(__name__)
     logger = setup_logger(logger)
 
+    # Quasi-RRHO constants (Grimme, Chem. Eur. J. 2012, 18, 9955): the damping
+    # frequency and the average molecular moment of inertia used for the
+    # free-rotor limit of low-frequency modes.
+    _QRRHO_FREQ_CM = 100.0  # cm^-1
+    _QRRHO_BAV = 1.0e-44  # kg m^2
+
     def __init__(
-        self, temperature: float, pressure: float, system: System, engine: str
+        self,
+        temperature: float,
+        pressure: float,
+        system: System,
+        engine: str,
+        quasi_rrho: bool = False,
     ):
         """
         Initializes the Thermo class with the temperature, pressure, system information
@@ -66,6 +77,10 @@ def __init__(
         engine : str
             The engine used for the calculation to compute the thermochemical
             properties with correct units.
+        quasi_rrho : bool
+            If True, use Grimme's quasi-RRHO treatment for the vibrational
+            entropy (interpolating low-frequency modes towards a free rotor)
+            instead of the pure rigid-rotor-harmonic-oscillator model.
 
         Raises
         ------
@@ -96,6 +111,7 @@ def __init__(
         self._pressure = pressure
         self._system = system
         self._engine = engine
+        self._quasi_rrho = quasi_rrho
 
         if self._engine != "dftb+":
             raise TSValueError("The engine is not supported.")
@@ -352,26 +368,74 @@ def _compute_vibrational_entropy(self):
         """
         Computes the vibrational entropy of the system.
 
+        Uses the harmonic-oscillator model, or Grimme's quasi-RRHO interpolation
+        towards a free rotor for low-frequency modes when ``quasi_rrho`` is set.
+
         Returns
         -------
         None
         """
 
-        self._vibrational_entropy = np.subtract(
-            np.divide(
-                (self._vib_temp_K / self._temperature),
-                (np.exp(self._vib_temp_K / self._temperature) - 1),
-            ),
-            np.log(1 - np.exp(-self._vib_temp_K / self._temperature)),
+        # Per-mode harmonic-oscillator entropy (in units of R).
+        x = self._vib_temp_K / self._temperature
+        harmonic = np.subtract(
+            np.divide(x, (np.exp(x) - 1)),
+            np.log(1 - np.exp(-x)),
+        )
+
+        entropy_per_mode = (
+            self._quasi_rrho_entropy(harmonic)
+            if self._quasi_rrho
+            else harmonic
         )
 
         self._vibrational_entropy = (
             PhysicalConstants["R"]
-            * np.sum(self._vibrational_entropy)
+            * np.sum(entropy_per_mode)
             / PhysicalConstants["cal"]
         )
 
 
+    def _quasi_rrho_entropy(self, harmonic):
+        """
+        Grimme's quasi-RRHO per-mode entropy (Chem. Eur. J. 2012, 18, 9955).
+
+        Each mode's entropy is interpolated between the harmonic-oscillator value
+        and the free-rotor value, weighted by ``w = 1 / (1 + (nu0/nu)^4)`` so that
+        high-frequency modes stay harmonic while low-frequency modes (whose HO
+        entropy diverges as nu -> 0) approach the finite free-rotor limit.
+
+        Parameters
+        ----------
+        harmonic : np.ndarray
+            Per-mode harmonic-oscillator entropy in units of R.
+
+        Returns
+        -------
+        np.ndarray
+            Per-mode quasi-RRHO entropy in units of R.
+        """
+        h = PhysicalConstants["h"]
+        kB = PhysicalConstants["kB"]
+        temperature = self._temperature
+
+        # Moment of inertia of each mode (h*nu = kB*theta, so mu = h^2 / (8 pi^2 kB theta)),
+        # then damped towards an average molecular moment B_av for the free rotor.
+        moment = h**2 / (8 * np.pi**2 * kB * self._vib_temp_K)
+        effective_moment = moment * self._QRRHO_BAV / (moment + self._QRRHO_BAV)
+
+        free_rotor = 0.5 + np.log(
+            np.sqrt(8 * np.pi**3 * effective_moment * kB * temperature / h**2)
+        )
+
+        weight = 1.0 / (
+            1.0
+            + (self._QRRHO_FREQ_CM / self._system.real_vibrational_frequencies) ** 4
+        )
+
+        return weight * harmonic + (1.0 - weight) * free_rotor
+
+
     def _compute_vibrational_energy(self):
         """
         Computes the vibrational energy of the system and the zero point energy correction.

tests/thermo/test_api.py

@@ -408,8 +408,9 @@ class FakeModes:
         def __init__(self):
             self.wave_numbers = np.array([1.0, 2.0, 3.0])
 
-    def fake_run_thermo(frequencies, atoms=None, spin=None, **kwargs):
+    def fake_run_thermo(frequencies, atoms=None, spin=None, quasi_rrho=False, **kwargs):
         seen["run_thermo_spin"] = spin
+        seen["run_thermo_quasi_rrho"] = quasi_rrho
         return "thermo-result"
 
     monkeypatch.setattr(api, "Geoopt", FakeGeoopt)
@@ -498,5 +499,16 @@ def test_dftbplus_thermo_gas_phase_has_no_solvation(monkeypatch, tmp_path):
     assert "Hamiltonian_Solvation" not in seen["geoopt_kwargs"]
 
 
+def test_dftbplus_thermo_forwards_quasi_rrho(monkeypatch, tmp_path):
+    api, seen = _mock_pipeline(monkeypatch)
+    api.dftbplus_thermo(
+        Atoms("OH2", positions=[[0, 0, 0.12], [0, 0.76, -0.48], [0, -0.76, -0.48]]),
+        directory=str(tmp_path / "j"),
+        quasi_rrho=True,
+    )
+
+    assert seen["run_thermo_quasi_rrho"] is True
+
+
 if __name__ == "__main__":
     unittest.main()

tests/thermo/test_screening.py

@@ -146,6 +146,26 @@ def fake_thermo(atoms, solvent=None, **kwargs):
     assert captured["solvent"] == "water"
 
 
+def test_screen_passes_quasi_rrho_to_dftbplus_thermo(monkeypatch, tmp_path):
+    _write_xyz(tmp_path / "mol.xyz")
+
+    captured = {}
+
+    def fake_thermo(atoms, quasi_rrho=False, **kwargs):
+        captured["quasi_rrho"] = quasi_rrho
+        return _FakeThermo()
+
+    monkeypatch.setattr(screening, "dftbplus_thermo", fake_thermo)
+    screening.screen(
+        str(tmp_path),
+        out=str(tmp_path / "r"),
+        directory=str(tmp_path / "runs"),
+        quasi_rrho=True,
+    )
+
+    assert captured["quasi_rrho"] is True
+
+
 def test_load_jobs_rejects_unknown_source(tmp_path):
     bad = tmp_path / "thing.txt"
     bad.write_text("x", encoding="utf-8")
@@ -253,6 +273,10 @@ def test_cli_parse_args_routes_screen():
 
     solv_args = cli.parse_args(["screen", "molecules.csv", "--solvent", "water"])
     assert solv_args.solvent == "water"
+    assert solv_args.quasi_rrho is False  # harmonic by default
+
+    qrrho_args = cli.parse_args(["screen", "molecules.csv", "--quasi-rrho"])
+    assert qrrho_args.quasi_rrho is True
 
 
 def test_cli_run_screen_returns_failure_count(monkeypatch):
@@ -265,7 +289,7 @@ def test_cli_run_screen_returns_failure_count(monkeypatch):
     args = Namespace(
         source="x", out="res", charge=0.0, temperature=298.15,
         pressure=101325.0, directory="screening", parameter_set="3ob",
-        solvent=None,
+        solvent=None, quasi_rrho=False,
     )
 
     assert cli.run_screen(args) == 1  # one molecule failed

tests/thermo/test_thermo.py

@@ -126,14 +126,17 @@ def test_thermo_rejects_negative_pressure():
 _R_CALMOLK = 8.314462618 / 4.184
 
 
-def _ts_thermo(symbols, positions, real_freqs, dof):
+def _ts_thermo(symbols, positions, real_freqs, dof, quasi_rrho=False):
     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,
         electronic_energy=0.0, vibrational_frequencies=pad,
     )
-    thermo = Thermo(temperature=_T, pressure=_P, system=system, engine="dftb+")
+    thermo = Thermo(
+        temperature=_T, pressure=_P, system=system, engine="dftb+",
+        quasi_rrho=quasi_rrho,
+    )
     thermo.run()
     return thermo
 
@@ -168,6 +171,34 @@ def test_total_entropy_matches_ase(symbols, positions, freqs, dof, geometry, sig
     assert ts_total == pytest.approx(ase_total, abs=0.05)
 
 
+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"],
+            [[0, 0, 0.119], [0, 0.763, -0.477], [0, -0.763, -0.477]],
+            [1595.0, 3657.0, 3756.0], 3)
+    harmonic = _ts_thermo(*args).total_entropy("cal/(mol*K)")
+    qrrho = _ts_thermo(*args, quasi_rrho=True).total_entropy("cal/(mol*K)")
+
+    assert qrrho == pytest.approx(harmonic, abs=0.05)
+
+
+def test_quasi_rrho_reduces_low_frequency_entropy():
+    # a floppy molecule with very low modes: harmonic overestimates their entropy
+    # (S_HO -> inf as nu -> 0), quasi-RRHO tames it towards the free rotor
+    args = (["C", "C", "H", "H", "H", "H", "H", "H"],
+            [[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]],
+            [25.0, 40.0, 820.0, 995.0, 1206.0, 1388.0, 1469.0, 1479.0,
+             2896.0, 2915.0, 2954.0, 2969.0, 2985.0, 2985.0, 1206.0, 1486.0, 995.0, 300.0],
+            18)
+    harmonic = _ts_thermo(*args).total_entropy("cal/(mol*K)")
+    qrrho = _ts_thermo(*args, quasi_rrho=True).total_entropy("cal/(mol*K)")
+
+    # the low modes' spurious harmonic entropy is removed -> qRRHO is lower, finite
+    assert qrrho < harmonic - 1.0
+    assert np.isfinite(qrrho)
+
+
 def test_thermo_rejects_imaginary_vibrational_mode():
     # an imaginary (negative) mode in the kept dof set used to give NaN; it now
     # raises (the geometry is not a minimum)