From b4a56d2b27ee630419a9c5e77a6eaf14c052f4a2 Mon Sep 17 00:00:00 2001
From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com>
Date: Tue, 7 Jul 2026 13:58:04 +0200
Subject: [PATCH 1/2] Import QM outputs (Gaussian, Turbomole, ORCA, ...) via
cclib
Add an optional `qm` extra (cclib) and a cclib-backed reader that consumes a
frequency calculation from any cclib-supported program (Gaussian, Turbomole,
ORCA, Psi4, NWChem, ...):
- read_cclib(path) -> (ase.Atoms, frequencies cm^-1, energy Hartree), using
cclib's atomcoords (Angstrom), vibfreqs (cm^-1) and best available energy
(ccenergies/mpenergies/scfenergies, converted eV -> Hartree).
- cclib_thermo(output_file, energy=None, ...) runs the RRHO thermochemistry
via run_thermo (energy defaults to cclib's, overridable).
cclib is imported lazily, so the core install (and the conda-forge recipe) is
unchanged; a clear error asks for `thermoscreening[qm]` when it's absent.
Closes #100
---
ThermoScreening/calculator/__init__.py | 1 +
ThermoScreening/calculator/qm.py | 95 +++++++++++++++++++++++
ThermoScreening/thermo/api.py | 69 +++++++++++++++++
docs/api.rst | 2 +
docs/usage.rst | 14 ++++
pyproject.toml | 7 +-
tests/calculator/test_qm.py | 101 +++++++++++++++++++++++++
7 files changed, 288 insertions(+), 1 deletion(-)
create mode 100644 ThermoScreening/calculator/qm.py
create mode 100644 tests/calculator/test_qm.py
diff --git a/ThermoScreening/calculator/__init__.py b/ThermoScreening/calculator/__init__.py
index 7d3689d..aa54b07 100644
--- a/ThermoScreening/calculator/__init__.py
+++ b/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
diff --git a/ThermoScreening/calculator/qm.py b/ThermoScreening/calculator/qm.py
new file mode 100644
index 0000000..f2e8fa2
--- /dev/null
+++ b/ThermoScreening/calculator/qm.py
@@ -0,0 +1,95 @@
+"""Import QM outputs (Gaussian, Turbomole, ORCA, Psi4, NWChem, ...) via cclib.
+
+`cclib `_ 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 reports energies in eV).
+_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.
+
+ 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
diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py
index 66bdb9e..aefdfc6 100644
--- a/ThermoScreening/thermo/api.py
+++ b/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.
diff --git a/docs/api.rst b/docs/api.rst
index 11f7dee..d19c187 100644
--- a/docs/api.rst
+++ b/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
-------------
diff --git a/docs/usage.rst b/docs/usage.rst
index 624e996..dfb034c 100644
--- a/docs/usage.rst
+++ b/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 `_:
+
+.. 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
-----------------
diff --git a/pyproject.toml b/pyproject.toml
index 31bf9d7..88d824c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -43,10 +43,15 @@ 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
+qm = [
+ "cclib >= 1.8",
+]
test = [
"pytest",
"coverage",
- "pytest-cov"
+ "pytest-cov",
+ "cclib >= 1.8",
]
lint = [
"pylint >= 4.0, < 5",
diff --git a/tests/calculator/test_qm.py b/tests/calculator/test_qm.py
new file mode 100644
index 0000000..83727e8
--- /dev/null
+++ b/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"))
From d32977af17c0e6d51d8f1ecbac1ff049347bbaaa Mon Sep 17 00:00:00 2001
From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com>
Date: Tue, 7 Jul 2026 14:04:05 +0200
Subject: [PATCH 2/2] Pin cclib < 2.0 and document the vibfreqs/linearity
assumption
cclib 2.0 is a breaking rearchitecture (data model / units); cap the qm
(and test) extra below it since the reader targets the 1.x model. Note in
read_cclib that cclib supplies only the real modes, so the geometry's
linearity classification must match the QM program's mode count.
---
ThermoScreening/calculator/qm.py | 10 +++++++++-
pyproject.toml | 8 +++++---
2 files changed, 14 insertions(+), 4 deletions(-)
diff --git a/ThermoScreening/calculator/qm.py b/ThermoScreening/calculator/qm.py
index f2e8fa2..d0141fe 100644
--- a/ThermoScreening/calculator/qm.py
+++ b/ThermoScreening/calculator/qm.py
@@ -11,7 +11,7 @@
from ..exceptions import TSValueError
-# 1 Hartree in eV (cclib reports energies in eV).
+# 1 Hartree in eV (cclib 1.x reports energies in eV; the qm extra pins < 2.0).
_EV_PER_HARTREE = 27.211386245988
@@ -70,6 +70,14 @@ def read_cclib(path):
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
diff --git a/pyproject.toml b/pyproject.toml
index 88d824c..da94e9e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -43,15 +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
+# 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",
+ "cclib >= 1.8, < 2.0",
]
test = [
"pytest",
"coverage",
"pytest-cov",
- "cclib >= 1.8",
+ "cclib >= 1.8, < 2.0",
]
lint = [
"pylint >= 4.0, < 5",