Import QM outputs (Gaussian, Turbomole, ORCA, ...) via cclib

ThermoScreening/calculator/__init__.py

@@ -2,3 +2,4 @@
 
 from .dftbplus import Geoopt, Hessian, Modes
 from .orca import read_orca_hess
+from .qm import read_cclib

ThermoScreening/calculator/qm.py

@@ -0,0 +1,103 @@
+"""Import QM outputs (Gaussian, Turbomole, ORCA, Psi4, NWChem, ...) via cclib.
+
+`cclib <https://cclib.github.io/>`_ auto-detects the program that produced an
+output file and exposes a uniform data model, so a single reader consumes a
+frequency calculation from any of the programs it supports. cclib is an optional
+dependency; install it with ``pip install thermoscreening[qm]``.
+"""
+
+import numpy as np
+from ase import Atoms
+
+from ..exceptions import TSValueError
+
+# 1 Hartree in eV (cclib 1.x reports energies in eV; the qm extra pins < 2.0).
+_EV_PER_HARTREE = 27.211386245988
+
+
+def _import_cclib():
+    """Import cclib, or raise a helpful error if the optional dep is missing."""
+    try:
+        import cclib  # noqa: F401
+        import cclib.io
+    except ImportError as exc:  # pragma: no cover - exercised via monkeypatch
+        raise TSValueError(
+            "cclib is required to import QM outputs; install it with "
+            "'pip install thermoscreening[qm]'."
+        ) from exc
+    return cclib
+
+
+def _best_energy_ev(data):
+    """
+    The highest-level electronic energy available (eV), or ``None``.
+
+    Prefers coupled-cluster, then Moller-Plesset, then SCF -- matching the level
+    the Hessian was most likely computed at.
+    """
+    ccenergies = getattr(data, "ccenergies", None)
+    if ccenergies is not None and len(ccenergies):
+        return float(ccenergies[-1])
+    mpenergies = getattr(data, "mpenergies", None)
+    if mpenergies is not None and len(mpenergies):
+        return float(mpenergies[-1][-1])  # last step, highest MP order
+    scfenergies = getattr(data, "scfenergies", None)
+    if scfenergies is not None and len(scfenergies):
+        return float(scfenergies[-1])
+    return None
+
+
+def read_cclib(path):
+    """
+    Read geometry, vibrational frequencies and energy from a QM output file.
+
+    Uses cclib to parse any supported program's output (Gaussian, Turbomole,
+    ORCA, Psi4, NWChem, ...).
+
+    Parameters
+    ----------
+    path : str
+        Path to a QM frequency-calculation output file.
+
+    Returns
+    -------
+    atoms : ase.Atoms
+        The final geometry (cclib reports coordinates in Angstrom).
+    frequencies : np.ndarray
+        The vibrational frequencies in cm^-1 (cclib gives the real vibrational
+        modes; imaginary modes appear as negatives).
+    energy : float or None
+        The electronic energy in Hartree (converted from cclib's eV), or ``None``
+        if the output has no parseable energy.
+
+    Notes
+    -----
+    cclib returns only the real vibrational modes (3N-6, or 3N-5 for a linear
+    molecule). The thermochemistry keeps the top ``dof`` modes for the geometry's
+    own linearity classification, so this matches for clearly (non)linear
+    molecules; a near-linear geometry where the classification disagrees with the
+    QM program's mode count could drop or miss a mode.
+
+    Raises
+    ------
+    TSValueError
+        If cclib is not installed, the file cannot be parsed, or it has no
+        vibrational frequencies.
+    """
+    cclib = _import_cclib()
+
+    data = cclib.io.ccread(str(path))
+    if data is None:
+        raise TSValueError(f"cclib could not parse '{path}' as a QM output.")
+    if getattr(data, "vibfreqs", None) is None or not len(data.vibfreqs):
+        raise TSValueError(
+            f"'{path}' has no vibrational frequencies (run a frequency calculation)."
+        )
+
+    atoms = Atoms(numbers=np.asarray(data.atomnos), positions=np.asarray(data.atomcoords[-1]))
+    frequencies = np.asarray(data.vibfreqs, dtype=float)
+
+    energy_ev = _best_energy_ev(data)
+    energy = energy_ev / _EV_PER_HARTREE if energy_ev is not None else None
+
+    return atoms, frequencies, energy

ThermoScreening/thermo/api.py

@@ -20,6 +20,7 @@
 from ..calculator import Geoopt, Hessian, Modes
 from ..calculator.dftbplus import _spin_kwargs, _solvation_kwargs, _dispersion_kwargs, SPIN_CONSTANTS_3OB
 from ..calculator.orca import read_orca_hess
+from ..calculator.qm import read_cclib
 from ..calculator.xtb import optimise_and_frequencies, xtb_calculator
 from ..calculator.xtb_cli import run_xtb
 
@@ -499,6 +500,74 @@ def orca_thermo(
     )
 
 
+def cclib_thermo(
+    output_file,
+    energy=None,
+    temperature=298.15,
+    pressure=101325,
+    charge=0.0,
+    spin=None,
+    quasi_rrho=False,
+):
+    """
+    Run the thermochemistry from a QM output file via cclib.
+
+    Consumes a frequency calculation from any cclib-supported program (Gaussian,
+    Turbomole, ORCA, Psi4, NWChem, ...) and evaluates the RRHO thermochemistry on
+    its DFT-quality geometry, frequencies and energy. Requires the optional
+    ``qm`` extra (``pip install thermoscreening[qm]``).
+
+    Parameters
+    ----------
+    output_file : str
+        Path to a QM frequency-calculation output file.
+    energy : float, optional
+        Electronic energy in Hartree. Defaults to the best energy cclib parses
+        from the file; pass this to override it (e.g. a higher-level single
+        point) or when cclib finds no energy.
+    temperature : float
+        Temperature in K. Default 298.15.
+    pressure : float
+        Pressure in Pa. Default 101325.
+    charge : float
+        System charge. Default 0.0.
+    spin : float, optional
+        Spin quantum number S. Defaults to the minimum-spin electron-count guess.
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO vibrational entropy. Default False.
+
+    Returns
+    -------
+    Thermo
+        The thermo calculation object.
+
+    Raises
+    ------
+    TSValueError
+        If cclib is missing, the file cannot be parsed, or no energy is available.
+    """
+    atoms, frequencies, file_energy = read_cclib(output_file)
+    if energy is None:
+        energy = file_energy
+    if energy is None:
+        raise TSValueError(
+            f"No energy for '{output_file}': cclib parsed none, so pass energy=... "
+            "explicitly."
+        )
+
+    return run_thermo(
+        vibrational_frequencies=frequencies,
+        atoms=atoms,
+        energy=energy,
+        temperature=temperature,
+        pressure=pressure,
+        charge=charge,
+        spin=spin,
+        engine="dftb+",
+        quasi_rrho=quasi_rrho,
+    )
+
+
 def execute(input_file: str) -> Thermo:
     """
     Execute the thermo calculation. Returns thermo calculation object.

docs/api.rst

@@ -8,6 +8,7 @@ High-level pipeline
 .. autofunction:: ThermoScreening.thermo.api.xtb_thermo
 .. autofunction:: ThermoScreening.thermo.api.xtb_cli_thermo
 .. autofunction:: ThermoScreening.thermo.api.orca_thermo
+.. autofunction:: ThermoScreening.thermo.api.cclib_thermo
 .. autofunction:: ThermoScreening.thermo.api.run_thermo
 .. autofunction:: ThermoScreening.thermo.api.execute
 
@@ -56,6 +57,7 @@ Coordinate and frequency readers
 .. autofunction:: ThermoScreening.thermo.api.read_coord
 .. autofunction:: ThermoScreening.thermo.api.read_vibrational
 .. autofunction:: ThermoScreening.calculator.orca.read_orca_hess
+.. autofunction:: ThermoScreening.calculator.qm.read_cclib
 
 Backend setup
 -------------

docs/usage.rst

@@ -121,6 +121,20 @@ Pass ``energy=`` to override the file's ``$act_energy`` (e.g. a higher-level
 single-point). ``read_orca_hess`` exposes the parsed geometry/frequencies/energy
 directly if you need them.
 
+For **Gaussian, Turbomole, Psi4, NWChem** (and ORCA) in one call, install the
+``qm`` extra (``pip install thermoscreening[qm]``) and use ``cclib_thermo``,
+which auto-detects the program via `cclib <https://cclib.github.io/>`_:
+
+.. code-block:: python
+
+   from ThermoScreening.thermo.api import cclib_thermo
+
+   thermo = cclib_thermo("freq.log")        # Gaussian, Turbomole, ORCA, ...
+   print(thermo.total_gibbs_free_energy())  # energy read from the output, in Hartree
+
+``read_cclib`` returns the parsed ``(atoms, frequencies, energy)`` if you want
+them directly; ``energy=`` overrides the parsed energy.
+
 Temperature scans
 -----------------
 

pyproject.toml

@@ -43,10 +43,17 @@ Repository = "https://github.com/MolarVerse/ThermoScreening"
 Issues = "https://github.com/MolarVerse/ThermoScreening/issues"
 
 [project.optional-dependencies]
+# QM-output import (Gaussian, Turbomole, ORCA, Psi4, NWChem, ...) via cclib.
+# Capped below the 2.0 rearchitecture (breaking data model / units handling);
+# the reader targets the cclib 1.x model (plain arrays, energies in eV).
+qm = [
+    "cclib >= 1.8, < 2.0",
+]
 test = [
     "pytest",
     "coverage",
-    "pytest-cov"
+    "pytest-cov",
+    "cclib >= 1.8, < 2.0",
 ]
 lint = [
     "pylint >= 4.0, < 5",

tests/calculator/test_qm.py

@@ -0,0 +1,101 @@
+import math
+import types
+
+import numpy as np
+import pytest
+
+import cclib.io
+
+from ThermoScreening.calculator import qm
+from ThermoScreening.calculator.qm import read_cclib, _best_energy_ev
+from ThermoScreening.thermo.api import cclib_thermo
+from ThermoScreening.exceptions import TSValueError
+
+_H_TO_EV = 27.211386245988
+
+# water: geometry (Angstrom), three real modes, SCF energy in eV (~ -76.4 Ha)
+_WATER = dict(
+    atomnos=np.array([8, 1, 1]),
+    atomcoords=np.array([[[0.0, 0.0, 0.1173],
+                          [0.0, 0.7572, -0.4692],
+                          [0.0, -0.7572, -0.4692]]]),
+    vibfreqs=np.array([1600.0, 3700.0, 3800.0]),
+    scfenergies=np.array([-76.4 * _H_TO_EV]),
+)
+
+
+def _fake_ccread(monkeypatch, **data):
+    monkeypatch.setattr(cclib.io, "ccread", lambda path: types.SimpleNamespace(**data))
+
+
+def _fake_ccread_returns(monkeypatch, value):
+    monkeypatch.setattr(cclib.io, "ccread", lambda path: value)
+
+
+def test_read_cclib_atoms_frequencies_energy(monkeypatch, tmp_path):
+    _fake_ccread(monkeypatch, **_WATER)
+    atoms, freqs, energy = read_cclib(str(tmp_path / "water.log"))
+
+    assert list(atoms.get_chemical_symbols()) == ["O", "H", "H"]
+    assert atoms.positions[1, 1] == pytest.approx(0.7572)
+    assert list(freqs) == [1600.0, 3700.0, 3800.0]
+    assert energy == pytest.approx(-76.4)  # eV -> Hartree
+
+
+def test_best_energy_prefers_cc_then_mp_then_scf():
+    scf = types.SimpleNamespace(scfenergies=np.array([-1.0]))
+    assert _best_energy_ev(scf) == -1.0
+
+    mp = types.SimpleNamespace(scfenergies=np.array([-1.0]),
+                               mpenergies=np.array([[-2.0, -2.5]]))
+    assert _best_energy_ev(mp) == -2.5  # last step, highest MP order
+
+    cc = types.SimpleNamespace(scfenergies=np.array([-1.0]),
+                               mpenergies=np.array([[-2.0]]),
+                               ccenergies=np.array([-3.0]))
+    assert _best_energy_ev(cc) == -3.0
+
+    assert _best_energy_ev(types.SimpleNamespace()) is None
+
+
+def test_read_cclib_unparseable_raises(monkeypatch, tmp_path):
+    _fake_ccread_returns(monkeypatch, None)
+    with pytest.raises(TSValueError, match="could not parse"):
+        read_cclib(str(tmp_path / "x.log"))
+
+
+def test_read_cclib_without_frequencies_raises(monkeypatch, tmp_path):
+    _fake_ccread(monkeypatch, atomnos=np.array([1]),
+                 atomcoords=np.array([[[0.0, 0.0, 0.0]]]), vibfreqs=np.array([]))
+    with pytest.raises(TSValueError, match="no vibrational frequencies"):
+        read_cclib(str(tmp_path / "x.log"))
+
+
+def test_read_cclib_missing_dependency(monkeypatch, tmp_path):
+    def _raise():
+        raise TSValueError("cclib is required to import QM outputs; install ...")
+
+    monkeypatch.setattr(qm, "_import_cclib", _raise)
+    with pytest.raises(TSValueError, match="cclib is required"):
+        read_cclib(str(tmp_path / "x.log"))
+
+
+def test_cclib_thermo_uses_file_energy(monkeypatch, tmp_path):
+    _fake_ccread(monkeypatch, **_WATER)
+    thermo = cclib_thermo(str(tmp_path / "water.log"))
+    assert thermo.electronic_energy() == pytest.approx(-76.4)
+    assert math.isfinite(thermo.total_EeGtot())
+
+
+def test_cclib_thermo_energy_override(monkeypatch, tmp_path):
+    _fake_ccread(monkeypatch, **_WATER)
+    thermo = cclib_thermo(str(tmp_path / "water.log"), energy=-77.0)
+    assert thermo.electronic_energy() == pytest.approx(-77.0)
+
+
+def test_cclib_thermo_requires_energy(monkeypatch, tmp_path):
+    data = dict(_WATER)
+    del data["scfenergies"]  # no energy anywhere and none passed
+    _fake_ccread(monkeypatch, **data)
+    with pytest.raises(TSValueError, match="No energy"):
+        cclib_thermo(str(tmp_path / "water.log"))