Support Turbomole's multi-file output in read_cclib/cclib_thermo

ThermoScreening/calculator/qm.py

@@ -47,17 +47,40 @@ def _best_energy_ev(data):
     return None
 
 
+def _normalize_source(path):
+    """
+    Coerce ``path`` into what ``cclib.io.ccread`` expects: a single path string,
+    or a list of path strings for a multi-file program.
+    """
+    if isinstance(path, (list, tuple)):
+        return [str(p) for p in path]
+    return str(path)
+
+
+def _describe_source(path):
+    """A short, readable description of ``path`` for error messages."""
+    if isinstance(path, (list, tuple)):
+        return ", ".join(str(p) for p in path)
+    return str(path)
+
+
 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, ...).
+    Uses cclib to parse any supported program's output (Gaussian, ORCA, Psi4,
+    NWChem, ...). Most programs write one logfile per job; pass its path.
+
+    **Turbomole** splits a job's output across many small files instead of one
+    logfile (``control``, ``coord``, ``aoforce.out``, ...); pass a list of every
+    relevant file's path (cclib's own multi-file mode) rather than a single path.
 
     Parameters
     ----------
-    path : str
-        Path to a QM frequency-calculation output file.
+    path : str or list of str
+        Path to a QM frequency-calculation output file, or (for a multi-file
+        program such as Turbomole) a list of paths to every relevant file from
+        the same job.
 
     Returns
     -------
@@ -81,17 +104,18 @@ def read_cclib(path):
     Raises
     ------
     TSValueError
-        If cclib is not installed, the file cannot be parsed, or it has no
+        If cclib is not installed, the file(s) cannot be parsed, or there are no
         vibrational frequencies.
     """
     cclib = _import_cclib()
 
-    data = cclib.io.ccread(str(path))
+    data = cclib.io.ccread(_normalize_source(path))
     if data is None:
-        raise TSValueError(f"cclib could not parse '{path}' as a QM output.")
+        raise TSValueError(f"cclib could not parse '{_describe_source(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)."
+            f"'{_describe_source(path)}' has no vibrational frequencies "
+            "(run a frequency calculation)."
         )
 
     atoms = Atoms(numbers=np.asarray(data.atomnos), positions=np.asarray(data.atomcoords[-1]))

ThermoScreening/thermo/api.py

@@ -20,7 +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.qm import read_cclib, _describe_source
 from ..calculator.xtb import optimise_and_frequencies, xtb_calculator
 from ..calculator.xtb_cli import run_xtb
 
@@ -540,11 +540,14 @@ def cclib_thermo(
 
     Parameters
     ----------
-    output_file : str
-        Path to a QM frequency-calculation output file.
+    output_file : str or list of str
+        Path to a QM frequency-calculation output file. Turbomole splits a job's
+        output across many small files instead of one logfile (``control``,
+        ``coord``, ``aoforce.out``, ...); pass a list of every relevant file's
+        path for it (cclib's multi-file mode) rather than a single path.
     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
+        from the file(s); 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.
@@ -578,8 +581,8 @@ def cclib_thermo(
         energy = file_energy
     if energy is None:
         raise TSValueError(
-            f"No energy for '{output_file}': cclib parsed none, so pass energy=... "
-            "explicitly."
+            f"No energy for '{_describe_source(output_file)}': cclib parsed none, "
+            "so pass energy=... explicitly."
         )
 
     return run_thermo(

docs/usage.rst

@@ -131,9 +131,13 @@ which auto-detects the program via `cclib <https://cclib.github.io/>`_:
 
    from ThermoScreening.thermo.api import cclib_thermo
 
-   thermo = cclib_thermo("freq.log")        # Gaussian, Turbomole, ORCA, ...
+   thermo = cclib_thermo("freq.log")        # Gaussian, ORCA, Psi4, NWChem, ...
    print(thermo.total_gibbs_free_energy())  # energy read from the output, in Hartree
 
+   # Turbomole splits a job's output across many files instead of one logfile;
+   # pass every relevant file's path as a list (cclib's multi-file mode)
+   ts_thermo = cclib_thermo(["control", "coord", "aoforce.out"])
+
 ``read_cclib`` returns the parsed ``(atoms, frequencies, energy)`` if you want
 them directly; ``energy=`` overrides the parsed energy.
 

tests/calculator/test_qm.py

@@ -1,18 +1,30 @@
 import math
 import types
+from pathlib import Path
 
 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.calculator.qm import (
+    read_cclib,
+    _best_energy_ev,
+    _normalize_source,
+    _describe_source,
+)
 from ThermoScreening.thermo.api import cclib_thermo
 from ThermoScreening.exceptions import TSValueError
 
 _H_TO_EV = 27.211386245988
 
+_REAL_TURBOMOLE_FILES = sorted(
+    str(p) for p in
+    (Path(__file__).resolve().parents[1] / "data" / "calculator" / "turbomole").glob("*")
+    if p.suffix != ".md"
+)
+
 # water: geometry (Angstrom), three real modes, SCF energy in eV (~ -76.4 Ha)
 _WATER = dict(
     atomnos=np.array([8, 1, 1]),
@@ -106,3 +118,72 @@ def test_cclib_thermo_requires_energy(monkeypatch, tmp_path):
     _fake_ccread(monkeypatch, **data)
     with pytest.raises(TSValueError, match="No energy"):
         cclib_thermo(str(tmp_path / "water.log"))
+
+
+def test_cclib_thermo_requires_energy_multi_file_error_is_readable(monkeypatch):
+    data = dict(_WATER)
+    del data["scfenergies"]
+    _fake_ccread(monkeypatch, **data)
+    with pytest.raises(TSValueError, match=r"No energy for 'control, coord, aoforce\.out'"):
+        cclib_thermo(["control", "coord", "aoforce.out"])
+
+
+# --- multi-file (Turbomole) support --- #
+#
+# Turbomole splits a job's output across many small files instead of one
+# logfile, so cclib.io.ccread must receive a LIST of paths for it, not a
+# single path. read_cclib previously always did cclib.io.ccread(str(path)),
+# which silently could not support this at all.
+
+
+def test_normalize_source_single_path_is_a_string():
+    assert _normalize_source("a.log") == "a.log"
+    assert _normalize_source(Path("a.log")) == "a.log"
+
+
+def test_normalize_source_list_stays_a_list_of_strings():
+    assert _normalize_source(["control", Path("coord"), "aoforce.out"]) == [
+        "control", "coord", "aoforce.out",
+    ]
+
+
+def test_describe_source_formats_both_forms():
+    assert _describe_source("a.log") == "a.log"
+    assert _describe_source(["a", "b"]) == "a, b"
+
+
+def test_read_cclib_passes_a_list_through_to_ccread_unmodified(monkeypatch):
+    captured = {}
+
+    def fake_ccread(source):
+        captured["source"] = source
+        return types.SimpleNamespace(**_WATER)
+
+    monkeypatch.setattr(cclib.io, "ccread", fake_ccread)
+    files = ["control", "coord", "aoforce.out"]
+    read_cclib(files)
+
+    assert captured["source"] == files  # not str(files) / a single joined string
+
+
+def test_read_cclib_real_turbomole_output():
+    # a genuine Turbomole 7.2 aoforce (frequency) calculation; see
+    # tests/data/calculator/turbomole/README.md for provenance
+    atoms, freqs, energy = read_cclib(_REAL_TURBOMOLE_FILES)
+
+    assert list(atoms.get_chemical_symbols()) == ["Cl", "Au", "N", "N", "Au", "N", "N"]
+    assert len(freqs) == 15  # 3*7 - 6, all real (a genuine minimum)
+    assert freqs.min() > 0
+    assert freqs.min() == pytest.approx(17.75)
+    assert freqs.max() == pytest.approx(2303.92)
+    # cclib's scfenergies (eV) converted to Hartree, cross-checked against the
+    # raw eV value cclib itself reports for this file (-25880.26134295 eV)
+    assert energy == pytest.approx(-25880.26134295 / _H_TO_EV)
+
+
+def test_cclib_thermo_real_turbomole_output():
+    thermo = cclib_thermo(_REAL_TURBOMOLE_FILES)
+
+    assert thermo.electronic_energy() == pytest.approx(-951.0820621320888)
+    assert math.isfinite(thermo.total_EeGtot())
+    assert thermo.total_entropy("cal/(mol*K)") > 0

tests/data/calculator/turbomole/README.md

@@ -0,0 +1,7 @@
+Real Turbomole 7.2 `aoforce` output (a frequency calculation on a 7-atom
+gold/azide/chloride complex), used to test `read_cclib`/`cclib_thermo`'s
+multi-file mode against a genuine, non-single-logfile QM program.
+
+Source: [cclib/cclib-data](https://github.com/cclib/cclib-data), the cclib
+project's own public regression-test data,
+`Turbomole/Turbomole7.2/au2_n22_cl_+_bp86-d3bj-deftzvp/`.