From 1d6e7f64d753e51165ab2ac05478b017ce1cd233 Mon Sep 17 00:00:00 2001 From: rlaplaza Date: Thu, 15 Oct 2026 10:08:49 +0200 Subject: [PATCH 1/5] Add ORCA plugin for vmol2 with view:freq|opt|sp|scan|irc selectors Multi-job ORCA logs can be parsed resiliently (ORCA 6 / PyPI cclib) and opened selectively instead of dumping frequencies and trajectories together. Co-authored-by: Cursor --- .github/CHANGELOG.md | 1 + README.md | 8 +- python/README.md | 25 ++- python/pyproject.toml | 2 +- python/ruff.toml | 1 + python/tests/test_orca_plugin.py | 257 +++++++++++++++++++++ python/vmol/plugins/__init__.py | 12 + python/vmol/plugins/orca.py | 361 ++++++++++++++++++++++++++++++ python/vmol/plugins/orca_parse.py | 246 ++++++++++++++++++++ python/vmol/vmol2.py | 300 +++++++++++++++++++++---- 10 files changed, 1168 insertions(+), 45 deletions(-) create mode 100644 python/tests/test_orca_plugin.py create mode 100644 python/vmol/plugins/__init__.py create mode 100644 python/vmol/plugins/orca.py create mode 100644 python/vmol/plugins/orca_parse.py diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 8d46c78..cb5f061 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -37,6 +37,7 @@ ### Coming in the next version: * high-symmetry determination bugs (#21) * headless fallback build by @JamesOBrien2 (#44) +* ORCA plugin for `vmol2`: prompt / `view:` selectors for frequencies, optimization, single point, scan, and IRC sidecar trajectories **Full Changelog**: https://github.com/briling/v/compare/v2.0...v3.0rc7 diff --git a/README.md b/README.md index 4d6a6f7..2bd89d0 100644 --- a/README.md +++ b/README.md @@ -10,14 +10,18 @@ A simple X11 molecular viewer. - [Priroda](http://rad.chem.msu.ru/~laikov) input and output files - [`.xyz`](https://en.wikipedia.org/wiki/XYZ_file_format) files - [extended `.xyz`](https://github.com/libAtoms/extxyz) files (currently the extra columns are ignored) -- various quantum-chemical outputs with [cclib](https://github.com/cclib/cclib), see the Python wrapper [page](python/README.md). +- various quantum-chemical outputs with [cclib](https://github.com/cclib/cclib), see the Python wrapper [page](python/README.md) + (`vmol2`; ORCA gets a plugin for frequencies, optimization, single point, scan, and IRC) ## Python package (wrapper / API) available See python package page [here](python/README.md). Provides wrapper scripts with a simple installation and -allows to open unsupported file formats with `cclib`. +allows to open unsupported file formats with `cclib` (`vmol2`). +For ORCA multi-job logs, `vmol2` can prompt (or take `view:` / `vib:`) for +frequencies, optimization, single point, IRC, or scan +(see [python/README.md](python/README.md)). ## Download ```bash diff --git a/python/README.md b/python/README.md index 784d69c..d760b6d 100644 --- a/python/README.md +++ b/python/README.md @@ -46,12 +46,28 @@ vmol ../mol/MOL_3525.xyz cell:8.93,0.0,0.0,4.2,8.9,0.0,0.48,2.32,10 python -m vmol../mol/periodic.in bonds:0 ``` -Unsupported quantum-chemical outputs, e.g., Orca: +Unsupported quantum-chemical outputs, e.g., ORCA: ```bash vmol2 ../mol/CEHZOF_1_SPE.out python -m vmol.vmol2 ../mol/AJALIH_5_SPE.out ``` -or Gaussian: +ORCA logs use a dedicated plugin (resilient to ORCA 6 / PyPI cclib). When a +file has several of frequencies, optimization, trailing single point, scan, or +IRC, `vmol2` prompts on a TTY (each choice opens in its own `v` session), or +accepts: +```bash +vmol2 job.out view:freq # normal modes (on the opt minimum) +vmol2 job.out view:opt # geometry optimization frames +vmol2 job.out view:sp # final single-point geometry / energy +vmol2 job.out view:irc # IRC (needs basename_IRC_Full_trj.xyz beside the .out) +vmol2 job.out view:scan # relaxed surface scan +vmol2 job.out vib:0 # geometry path (opt, else SP, else IRC, else scan) +vmol2 job.out vib:1 # frequencies +``` +IRC frames come from `basename_IRC_Full_trj.xyz` (or the path in the log's +`Storing full IRC trajectory in ....` line), not from the `.out` itself. + +Gaussian and other packages use the generic cclib path: ```bash vmol2 mol/Cat1_TSOA_CfA_Cf0.log vib:0 # geometry optimization vmol2 mol/Cat1_TSOA_CfA_Cf0.log # normal modes (default for files containing both) @@ -154,8 +170,9 @@ out = vmol.capture(mols=mols) ``` File formats which are not supported natively can be read with `cclib` and passed -(see [example](https://github.com/briling/v/tree/master/python/examples/ex_cclib.py) -and [vmol2 source](https://github.com/briling/v/tree/master/python/vmol/vmol2.py)). +(see [example](https://github.com/briling/v/tree/master/python/examples/ex_cclib.py), +[vmol2](https://github.com/briling/v/tree/master/python/vmol/vmol2.py), and the +ORCA [plugin](https://github.com/briling/v/tree/master/python/vmol/plugins/orca.py)). One can also pass normal modes data with ```python diff --git a/python/pyproject.toml b/python/pyproject.toml index 947ca37..668e681 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -46,4 +46,4 @@ Repository = "https://github.com/briling/v.git" "Bug Tracker" = "https://github.com/briling/v/issues" [tool.setuptools] -packages = ["vmol"] +packages = ["vmol", "vmol.plugins"] diff --git a/python/ruff.toml b/python/ruff.toml index 6bb3f0f..19a6a94 100644 --- a/python/ruff.toml +++ b/python/ruff.toml @@ -53,3 +53,4 @@ ignore = [ [lint.per-file-ignores] "vmol/__main__.py" = ["D"] "examples/ex*.py" = ["D"] +"tests/*.py" = ["SLF", "PLC2701"] diff --git a/python/tests/test_orca_plugin.py b/python/tests/test_orca_plugin.py new file mode 100644 index 0000000..e1c55e2 --- /dev/null +++ b/python/tests/test_orca_plugin.py @@ -0,0 +1,257 @@ +"""Unit tests for the ORCA vmol2 plugin (no viewer / X11).""" + +from __future__ import annotations + +import os +import tempfile +import unittest +import warnings +from pathlib import Path +from types import SimpleNamespace + +import numpy as np + +from vmol import vmol2 +from vmol.plugins.orca import Dataset, OrcaPlugin, load_irc_xyz, resolve_irc_path + + +def _coords(*frames): + return np.asarray(frames, dtype=float) + + +def _base_data(nframes=1): + frames = [[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74 + 0.01 * i]] for i in range(nframes)] + return SimpleNamespace( + atomnos=np.array([1, 1]), + atomcoords=_coords(*frames), + scfenergies=np.array([-27.2 - 0.01 * i for i in range(nframes)]), + ) + + +class TestOrcaInventory(unittest.TestCase): + """Inventory detection for frequencies, optimization, scan, and IRC.""" + + def setUp(self): + """Create a temporary directory and plugin instance.""" + self.plugin = OrcaPlugin() + self.tmpdir = tempfile.TemporaryDirectory() + self.addCleanup(self.tmpdir.cleanup) + self.out_path = os.path.join(self.tmpdir.name, "job.out") + + def _write_out(self, text: str) -> str: + Path(self.out_path).write_text(text, encoding="utf-8") + return self.out_path + + def test_frequency_only(self): + """A frequency job yields a single freq dataset.""" + data = _base_data(1) + data.vibfreqs = np.array([4401.0]) + data.vibdisps = np.zeros((1, 2, 3)) + data.vibdisps[0, 0, 2] = 0.7 + data.vibdisps[0, 1, 2] = -0.7 + data.vibirs = np.array([0.0]) + data.vibrmasses = np.array([1.0]) + self._write_out("* O R C A *\n") + + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["freq"]) + self.assertEqual(len(datasets[0].mols), 1) + self.assertIsNotNone(datasets[0].vib) + self.assertEqual(len(datasets[0].vib["freq"]), 1) + + def test_opt_plus_freq(self): + """Opt+freq offers both datasets and truncates at optdone.""" + data = _base_data(4) + data.geovalues = np.ones((3, 4)) + data.optdone = [3] + data.vibfreqs = np.array([100.0, 200.0]) + data.vibdisps = np.zeros((2, 2, 3)) + data.vibdisps[:, 0, 0] = 1.0 + self._write_out("* O R C A *\n") + + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["freq", "opt"]) + opt = datasets[1] + self.assertEqual(len(opt.mols), 3) + self.assertIsNone(opt.vib) + self.assertIn("step 1/3", opt.mols[0]["name"]) + + def test_opt_truncates_when_optdone_bool(self): + """Bool optdone uses geovalues (+ final eval), dropping a later job frame.""" + data = _base_data(5) + data.geovalues = np.ones((3, 4)) + data.optdone = True + self._write_out("* O R C A *\n") + + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["opt"]) + self.assertEqual(len(datasets[0].mols), 4) + + def test_opt_freq_and_sp(self): + """Multi-job files expose frequencies, opt, and a selectable SP.""" + data = _base_data(5) + data.geovalues = np.ones((3, 4)) + data.optdone = True + data.vibfreqs = np.array([100.0]) + data.vibdisps = np.zeros((1, 2, 3)) + data.vibdisps[0, 0, 0] = 1.0 + data.sp_frame_index = 4 + data.sp_energy_eh = -100.5 + self._write_out("* O R C A *\n") + + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["freq", "opt", "sp"]) + self.assertEqual(len(datasets[1].mols), 4) + self.assertIn("E=-100.50000000 Eh", datasets[2].label) + np.testing.assert_allclose(datasets[0].mols[0]["r"], data.atomcoords[3]) + + def test_scan_not_labeled_as_opt(self): + """Relaxed scans are labeled scan, not opt.""" + data = _base_data(3) + data.scannames = ["R"] + data.scanparm = [(0.7,), (0.8,), (0.9,)] + data.geovalues = np.ones((3, 4)) + self._write_out("* O R C A *\n") + + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["scan"]) + self.assertIn("R=0.7", datasets[0].mols[0]["name"]) + + def test_single_point_geometry(self): + """A single-point job yields one geom dataset.""" + data = _base_data(1) + self._write_out("* O R C A *\n") + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["geom"]) + self.assertEqual(len(datasets[0].mols), 1) + + def test_missing_irc_file_warns(self): + """Mentioned but missing IRC xyz warns and is omitted.""" + data = _base_data(1) + self._write_out( + "* O R C A *\n" + "Intrinsic Reaction Coordinate Calculation\n" + "Storing full IRC trajectory in .... job_IRC_Full_trj.xyz\n", + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["geom"]) + self.assertTrue(any("IRC trajectory" in str(w.message) for w in caught)) + + def test_irc_sidecar_loaded(self): + """IRC xyz comments become frame names via cclib.""" + data = _base_data(1) + xyz_path = os.path.join(self.tmpdir.name, "job_IRC_Full_trj.xyz") + Path(xyz_path).write_text( + "2\n" + "Coordinates from ORCA-job job_IRC_Full E -1.100000000000\n" + "H 0.0 0.0 0.0\n" + "H 0.0 0.0 0.74\n" + "2\n" + "Coordinates from ORCA-job job_IRC_Full E -1.200000000000\n" + "H 0.0 0.0 0.0\n" + "H 0.0 0.0 0.80\n", + encoding="utf-8", + ) + self._write_out( + "* O R C A *\n" + "Storing full IRC trajectory in .... job_IRC_Full_trj.xyz\n", + ) + datasets = self.plugin.inventory(self.out_path, data) + self.assertEqual([d.key for d in datasets], ["irc"]) + irc = datasets[0] + self.assertEqual(len(irc.mols), 2) + self.assertIn("E -1.100000000000", irc.mols[0]["name"]) + self.assertIn("E -1.200000000000", irc.mols[1]["name"]) + + +class TestIrcHelpers(unittest.TestCase): + """IRC path resolution and xyz loading via cclib.""" + + def test_resolve_from_log_line(self): + """Log line filename is preferred when the file exists.""" + with tempfile.TemporaryDirectory() as tmp: + out_path = os.path.join(tmp, "calc.out") + xyz_name = "custom_IRC_Full_trj.xyz" + xyz_path = os.path.join(tmp, xyz_name) + Path(out_path).write_text( + f"Storing full IRC trajectory in .... {xyz_name}\n", + encoding="utf-8", + ) + Path(xyz_path).touch() + self.assertEqual(resolve_irc_path(out_path), xyz_path) + + def test_resolve_fallback_stem(self): + """Fallback uses _IRC_Full_trj.xyz next to the .out.""" + with tempfile.TemporaryDirectory() as tmp: + out_path = os.path.join(tmp, "myjob.out") + xyz_path = os.path.join(tmp, "myjob_IRC_Full_trj.xyz") + Path(out_path).write_text("* O R C A *\n", encoding="utf-8") + Path(xyz_path).touch() + self.assertEqual(resolve_irc_path(out_path), xyz_path) + + def test_load_irc_xyz_comments(self): + """XYZ comment lines from cclib become molecule names.""" + with tempfile.TemporaryDirectory() as tmp: + xyz_path = os.path.join(tmp, "trj.xyz") + Path(xyz_path).write_text( + "2\n" + "Coordinates from ORCA-job trj E -99.5\n" + "H 0 0 0\n" + "H 0 0 1\n", + encoding="utf-8", + ) + ds = load_irc_xyz(xyz_path) + self.assertEqual(ds.key, "irc") + self.assertEqual(len(ds.mols), 1) + self.assertEqual(ds.mols[0]["name"], "Coordinates from ORCA-job trj E -99.5") + + +class TestViewSelection(unittest.TestCase): + """Non-interactive dataset selection helpers in vmol2.""" + + def test_view_and_vib_selectors(self): + """View and vib selectors pick the matching dataset keys.""" + datasets = [ + Dataset( + key="freq", + label="f", + mols=[{"q": [1], "r": [[0, 0, 0]]}], + vib={"freq": [1], "disp": [[[1, 0, 0]]]}, + ), + Dataset( + key="opt", + label="o", + mols=[{"q": [1], "r": [[0, 0, 0]]}, {"q": [1], "r": [[0, 0, 1]]}], + ), + Dataset(key="sp", label="s", mols=[{"q": [1], "r": [[0, 0, 0]]}]), + Dataset(key="irc", label="i", mols=[{"q": [1], "r": [[0, 0, 0]]}]), + ] + self.assertEqual( + [d.key for d in vmol2._choose_datasets(datasets, ["opt"], None)], + ["opt"], + ) + self.assertEqual( + [d.key for d in vmol2._choose_datasets(datasets, ["sp"], None)], + ["sp"], + ) + self.assertEqual( + [d.key for d in vmol2._choose_datasets(datasets, [], True)], + ["freq"], + ) + self.assertEqual( + [d.key for d in vmol2._choose_datasets(datasets, [], False)], + ["opt"], + ) + + keys, vib, pops = vmol2._collect_view_selectors( + ["file.out", "view:sp", "vib:0", "bonds:0"], + ) + self.assertEqual(keys, ["sp"]) + self.assertIs(vib, False) + self.assertEqual(pops, [1]) + + +if __name__ == "__main__": + unittest.main() diff --git a/python/vmol/plugins/__init__.py b/python/vmol/plugins/__init__.py new file mode 100644 index 0000000..7c8f7a4 --- /dev/null +++ b/python/vmol/plugins/__init__.py @@ -0,0 +1,12 @@ +"""Format-specific plugins for reading quantum-chemical outputs.""" + +from vmol.plugins.orca import OrcaPlugin + +__all__ = ["OrcaPlugin", "plugin_for"] + + +def plugin_for(parser): + """Return a plugin for *parser*, or None if no dedicated plugin exists.""" + if getattr(parser, "logname", None) == "ORCA": + return OrcaPlugin() + return None diff --git a/python/vmol/plugins/orca.py b/python/vmol/plugins/orca.py new file mode 100644 index 0000000..017af55 --- /dev/null +++ b/python/vmol/plugins/orca.py @@ -0,0 +1,361 @@ +"""ORCA-specific dataset inventory for vmol2.""" + +from __future__ import annotations + +import os +import re +import warnings +from dataclasses import dataclass + +import numpy as np +from cclib.io.xyzreader import XYZ +from cclib.parser import utils + +from vmol.plugins.orca_parse import parse_orca + +# On-screen molecule name buffer in v (src/v/v.h STRLEN). +_NAME_MAX = 255 + +_IRC_STORE_RE = re.compile( + r"Storing full IRC trajectory in\s+\.+\s+(\S+)", + re.IGNORECASE, +) +_IRC_MENTION_RE = re.compile( + r"Intrinsic Reaction Coordinate|" + r"Storing (?:full|forward|backward) IRC|" + r"(?:FORWARD|BACKWARD) IRC|" + r"IRC PATH SUMMARY", + re.IGNORECASE, +) + +CCLIB2VMOL = { + "vibdisps": "disp", + "vibfreqs": "freq", + "vibrmasses": "mass", + "vibirs": "ints", +} + + +@dataclass(frozen=True) +class Dataset: + """One viewable dataset extracted from an ORCA output.""" + + key: str + label: str + mols: list + vib: dict | None = None + + +def _clip_name(name: str) -> str: + if len(name) <= _NAME_MAX: + return name + return name[: _NAME_MAX - 3] + "..." + + +def _has_frequencies(data) -> bool: + freqs = getattr(data, "vibfreqs", None) + disps = getattr(data, "vibdisps", None) + return freqs is not None and disps is not None and len(freqs) and len(disps) + + +def _vib_from_data(data) -> dict: + return {dst: getattr(data, src, None) for src, dst in CCLIB2VMOL.items()} + + +def _energy_suffix(data, index: int) -> str: + energies = getattr(data, "scfenergies", None) + if energies is None or index >= len(energies): + return "" + val = float(energies[index]) + if np.isnan(val): + return "" + eh = utils.convertor(val, "eV", "hartree") + return f" E={eh:.8f} Eh" + + +def _frame_mols(data, coords, *, name_prefix: str, name_for_index=None) -> list: + atomnos = data.atomnos + n = len(coords) + mols = [] + for i, r in enumerate(coords): + if name_for_index is not None: + name = name_for_index(i, r) + else: + name = f"{name_prefix} step {i + 1}/{n}{_energy_suffix(data, i)}" + mols.append({"q": atomnos, "r": r, "name": _clip_name(name)}) + return mols + + +def _opt_end_index(data) -> int | None: + optdone = getattr(data, "optdone", None) + n = len(data.atomcoords) + if isinstance(optdone, (list, tuple)) and optdone: + last = max(int(x) for x in optdone) + if 1 <= last <= n: + return last + return None + if optdone is True: + geovalues = getattr(data, "geovalues", None) + if geovalues is not None and len(geovalues): + end = len(geovalues) + 1 + return min(end, n) + return None + + +def _scan_irc_filename(out_path: str) -> tuple[str | None, bool]: + stored = None + mentioned = False + with open(out_path, encoding="utf-8", errors="replace") as handle: + for line in handle: + match = _IRC_STORE_RE.search(line) + if match: + stored = os.path.basename(match.group(1)) + mentioned = True + break + if not mentioned and _IRC_MENTION_RE.search(line): + mentioned = True + return stored, mentioned + + +def resolve_irc_path(out_path: str, log_text: str | None = None) -> str | None: + """Find the ORCA IRC full trajectory xyz next to *out_path*. + + Returns: + Absolute path to the trajectory file, or None if missing. + """ + directory = os.path.dirname(os.path.abspath(out_path)) + stem = os.path.splitext(os.path.basename(out_path))[0] + + if log_text is None: + stored, mentioned = _scan_irc_filename(out_path) + else: + match = _IRC_STORE_RE.search(log_text) + stored = os.path.basename(match.group(1)) if match else None + mentioned = stored is not None or bool(_IRC_MENTION_RE.search(log_text)) + + candidates = [] + if stored: + candidates.append(os.path.join(directory, stored)) + candidates.append(os.path.join(directory, f"{stem}_IRC_Full_trj.xyz")) + + for candidate in candidates: + if os.path.isfile(candidate): + return candidate + + if mentioned or stored: + warnings.warn( + f"ORCA IRC trajectory mentioned in {out_path} but xyz not found " + f"(tried: {', '.join(candidates)}).", + RuntimeWarning, + stacklevel=3, + ) + return None + + +def load_irc_xyz(xyz_path: str) -> Dataset: + """Load an ORCA IRC multi-xyz as a Dataset (via cclib). + + Returns: + Dataset: One frame per xyz image; names from comment lines. + """ + data = XYZ(xyz_path).parse() + comments = data.metadata.get("comments", []) + nframes = len(data.atomcoords) + + def name_for_index(i, _r): + if i < len(comments) and comments[i].strip(): + return comments[i].strip() + return f"IRC frame {i + 1}/{nframes}" + + return Dataset( + key="irc", + label=f"IRC trajectory ({nframes} frames, {os.path.basename(xyz_path)})", + mols=_frame_mols( + data, data.atomcoords, name_prefix="IRC", name_for_index=name_for_index, + ), + vib=None, + ) + + +def _freq_geometry_index(data) -> int: + end = _opt_end_index(data) + if end is not None and end > 0: + return end - 1 + return len(data.atomcoords) - 1 + + +def _frequencies_dataset(data, name: str) -> Dataset: + idx = _freq_geometry_index(data) + return Dataset( + key="freq", + label=f"Frequencies ({len(data.vibfreqs)} modes)", + mols=[{ + "q": data.atomnos, + "r": data.atomcoords[idx], + "name": _clip_name(f"{name} (frequencies)"), + }], + vib=_vib_from_data(data), + ) + + +def _is_scan(data) -> bool: + return ( + getattr(data, "scanparm", None) is not None + or getattr(data, "scannames", None) is not None + ) + + +def _opt_coords(data): + if _is_scan(data): + return None + has_opt = ( + getattr(data, "geovalues", None) is not None + or getattr(data, "optdone", None) is not None + ) + if not has_opt or len(data.atomcoords) <= 1: + return None + end = _opt_end_index(data) + coords = data.atomcoords[:end] if end is not None else data.atomcoords + if len(coords) <= 1: + return None + return coords + + +def _optimization_dataset(data, name: str) -> Dataset | None: + coords = _opt_coords(data) + if coords is None: + return None + mols = _frame_mols(data, coords, name_prefix=f"{name} opt") + return Dataset( + key="opt", + label=f"Optimization trajectory ({len(mols)} frames)", + mols=mols, + vib=None, + ) + + +def _scan_dataset(data, name: str) -> Dataset | None: + if not _is_scan(data) or not len(data.atomcoords): + return None + + scanparm = getattr(data, "scanparm", None) + scannames = getattr(data, "scannames", None) or [] + n = len(data.atomcoords) + + def name_for_index(i, _r): + base = f"{name} scan step {i + 1}/{n}" + if scanparm is not None and i < len(scanparm): + params = scanparm[i] + if scannames and len(scannames) == len(params): + parts = [f"{n_}={v}" for n_, v in zip(scannames, params, strict=False)] + else: + parts = [str(v) for v in params] + base += " " + ", ".join(parts) + return base + _energy_suffix(data, i) + + mols = _frame_mols( + data, data.atomcoords, name_prefix=f"{name} scan", name_for_index=name_for_index, + ) + return Dataset( + key="scan", + label=f"Relaxed scan ({len(mols)} frames)", + mols=mols, + vib=None, + ) + + +def _single_point_dataset(data, name: str) -> Dataset | None: + idx = getattr(data, "sp_frame_index", None) + if idx is None: + return None + if _opt_coords(data) is None and not _has_frequencies(data): + return None + + energy = getattr(data, "sp_energy_eh", None) + if energy is None: + suffix = _energy_suffix(data, idx) + else: + suffix = f" E={float(energy):.8f} Eh" + label_energy = suffix.strip() or "final single point" + return Dataset( + key="sp", + label=f"Single point ({label_energy})", + mols=[{ + "q": data.atomnos, + "r": data.atomcoords[idx], + "name": _clip_name(f"{name} (single point){suffix}"), + }], + vib=None, + ) + + +def _geometry_dataset(data, name: str) -> Dataset: + mols = [ + {"q": data.atomnos, "r": r, "name": _clip_name(name)} + for r in data.atomcoords + ] + n = len(mols) + return Dataset( + key="geom", + label=f"Geometry ({n} frame{'s' if n != 1 else ''})", + mols=mols, + vib=None, + ) + + +class OrcaPlugin: + """Detect frequencies, optimization, SP, scan, and IRC datasets in an ORCA log.""" + + parse = staticmethod(parse_orca) + + @staticmethod + def inventory(path: str, data) -> list[Dataset]: + """List viewable datasets for *path* / *data*. + + Returns: + list[Dataset]: Non-empty list of datasets. + + Raises: + RuntimeError: If no atomic coordinates are present in *data*. + """ + if not ( + hasattr(data, "atomcoords") + and hasattr(data, "atomnos") + and len(data.atomcoords) + ): + msg = f"No coordinates found in {path}." + raise RuntimeError(msg) + + name = f"ORCA {os.path.basename(path)}" + datasets: list[Dataset] = [] + + if _has_frequencies(data): + datasets.append(_frequencies_dataset(data, name)) + + opt = _optimization_dataset(data, name) + if opt is not None: + datasets.append(opt) + + sp = _single_point_dataset(data, name) + if sp is not None: + datasets.append(sp) + + scan = _scan_dataset(data, name) + if scan is not None: + datasets.append(scan) + + irc_path = resolve_irc_path(path) + if irc_path is not None: + try: + datasets.append(load_irc_xyz(irc_path)) + except (OSError, ValueError, AssertionError) as exc: + warnings.warn( + f"Failed to load IRC trajectory {irc_path}: {exc}", + RuntimeWarning, + stacklevel=2, + ) + + if not datasets: + datasets.append(_geometry_dataset(data, name)) + + return datasets diff --git a/python/vmol/plugins/orca_parse.py b/python/vmol/plugins/orca_parse.py new file mode 100644 index 0000000..a2636f1 --- /dev/null +++ b/python/vmol/plugins/orca_parse.py @@ -0,0 +1,246 @@ +"""Resilient ORCA parsing for visualization (coords, freqs, multi-job SP).""" + +from __future__ import annotations + +import contextlib +import warnings +from types import SimpleNamespace + +import numpy as np +from cclib.io import ccopen +from cclib.parser import orcaparser, utils +from cclib.parser.utils import PeriodicTable + +_PT = PeriodicTable() +_PATCHED = False +_CART_MIN_TOKENS = 4 + + +def _patch_cclib_orca() -> None: + """Harden cclib's ORCA parser for ORCA 6 / incomplete SCF blocks. + + PyPI cclib 1.8.1 aborts on empty ``scfvalues`` and on some ORCA 6 sections. + Skip those sections so geometries and frequencies still get collected. + """ + global _PATCHED + if _PATCHED: + return + + orig_append = orcaparser.ORCA._append_scfvalues_scftargets + + def patched_append(self, inputfile, line): + if not getattr(self, "scfvalues", None): + self.scfvalues = [[]] + elif len(self.scfvalues) == 0: + self.scfvalues.append([]) + try: + return orig_append(self, inputfile, line) + except Exception as exc: # noqa: BLE001 + warnings.warn( + f"ORCA SCF convergence summary skipped ({type(exc).__name__}: {exc})", + RuntimeWarning, + stacklevel=2, + ) + return None + + orig_extract = orcaparser.ORCA.extract + + def safe_extract(self, inputfile, line): + try: + return orig_extract(self, inputfile, line) + except Exception as exc: # noqa: BLE001 + warnings.warn( + f"ORCA section skipped ({type(exc).__name__}: {exc})", + RuntimeWarning, + stacklevel=2, + ) + return None + + orcaparser.ORCA._append_scfvalues_scftargets = patched_append + orcaparser.ORCA.extract = safe_extract + _PATCHED = True + + +def _read_cartesian_block(lines, start: int) -> tuple[list[int], list[list[float]], int]: + i = start + 1 + if i < len(lines) and set(lines[i].strip()) <= {"-"}: + i += 1 + atomnos: list[int] = [] + coords: list[list[float]] = [] + while i < len(lines) and lines[i].strip(): + parts = lines[i].split() + if len(parts) < _CART_MIN_TOKENS: + break + sym = parts[0] + if sym.endswith(">"): + i += 1 + continue + if sym not in _PT.number: + break + atomnos.append(int(_PT.number[sym])) + coords.append([float(parts[1]), float(parts[2]), float(parts[3])]) + i += 1 + return atomnos, coords, i + + +def _scan_jobs(path: str) -> list[dict]: + with open(path, encoding="utf-8", errors="replace") as handle: + lines = handle.readlines() + + jobs: list[dict] = [] + current = { + "index": 0, + "kind": "unknown", + "frames": [], + "energy": None, + } + + def flush(): + nonlocal current + if current["frames"] or current["kind"] != "unknown": + jobs.append(current) + current = { + "index": len(jobs), + "kind": "unknown", + "frames": [], + "energy": None, + } + + i = 0 + while i < len(lines): + line = lines[i] + if "JOB NUMBER" in line: + if current["frames"] or current["kind"] != "unknown": + flush() + current["index"] = len(jobs) + elif "Geometry Optimization Run" in line: + current["kind"] = "opt" + elif "Single Point Calculation" in line: + current["kind"] = "sp" + elif "VIBRATIONAL FREQUENCIES" in line and current["kind"] == "opt": + current["kind"] = "optfreq" + elif line.startswith("CARTESIAN COORDINATES (ANGSTROEM)"): + atomnos, coords, i = _read_cartesian_block(lines, i) + if coords: + current["frames"].append((atomnos, coords)) + continue + elif "FINAL SINGLE POINT ENERGY" in line: + parts = line.split() + if parts: + with contextlib.suppress(ValueError): + current["energy"] = float(parts[-1]) + i += 1 + + if current["frames"] or current["kind"] != "unknown": + jobs.append(current) + return jobs + + +def _as_namespace(data) -> SimpleNamespace: + ns = SimpleNamespace() + for key in ( + "natom", + "atomnos", + "atomcoords", + "scfenergies", + "geovalues", + "optdone", + "vibfreqs", + "vibdisps", + "vibirs", + "vibrmasses", + "scanparm", + "scannames", + "metadata", + ): + if hasattr(data, key): + setattr(ns, key, getattr(data, key)) + return ns + + +def _append_sp_from_jobs(ns: SimpleNamespace, jobs: list[dict]) -> None: + """Ensure the last single-point job geometry is available on *ns*.""" + sp_jobs = [j for j in jobs if j["kind"] == "sp" and j["frames"]] + if not sp_jobs: + return + + sp = sp_jobs[-1] + atomnos, coords = sp["frames"][-1] + coords_arr = np.asarray(coords, dtype=float) + ns.sp_energy_eh = sp["energy"] + + if not hasattr(ns, "atomcoords") or len(ns.atomcoords) == 0: + ns.atomnos = np.asarray(atomnos, dtype=int) + ns.atomcoords = coords_arr[None, ...] + ns.natom = len(atomnos) + ns.sp_frame_index = 0 + return + + last = np.asarray(ns.atomcoords[-1], dtype=float) + if last.shape == coords_arr.shape and np.allclose(last, coords_arr, atol=1e-6): + ns.sp_frame_index = len(ns.atomcoords) - 1 + return + + ns.atomcoords = np.concatenate([ns.atomcoords, coords_arr[None, ...]], axis=0) + ns.sp_frame_index = len(ns.atomcoords) - 1 + if sp["energy"] is None: + return + energies = list(getattr(ns, "scfenergies", []) or []) + ev = utils.convertor(float(sp["energy"]), "hartree", "eV") + while len(energies) < ns.sp_frame_index: + energies.append(np.nan) + energies.append(ev) + ns.scfenergies = np.asarray(energies, dtype=float) + + +def parse_orca(path: str): + """Parse an ORCA logfile for visualization. + + Returns: + tuple: ``(parser, data)`` where *data* has at least atomcoords/atomnos + when the file contains geometries. Frequencies are included when present. + ``data.sp_frame_index`` is set when a trailing single-point job exists. + + Raises: + RuntimeError: If the file cannot be opened or has no coordinates. + """ + _patch_cclib_orca() + + parser = ccopen(path) + if parser is None: + msg = f"Could not open {path} with cclib." + raise RuntimeError(msg) + + try: + raw = parser.parse() + except Exception as exc: # noqa: BLE001 + warnings.warn( + f"Error parsing {path} with cclib ({type(exc).__name__}: {exc}). " + "Using partial data.", + RuntimeWarning, + stacklevel=2, + ) + raw = parser + + data = _as_namespace(raw) + jobs = _scan_jobs(path) + data.orca_jobs = jobs + _append_sp_from_jobs(data, jobs) + + coords = getattr(data, "atomcoords", None) + if coords is not None and len(coords) > 0: + return parser, data + + frames = [] + atomnos = None + for job in jobs: + for nos, xyz in job["frames"]: + atomnos = nos + frames.append(xyz) + if not frames: + msg = f"No coordinates found in {path}." + raise RuntimeError(msg) + data.atomnos = np.asarray(atomnos, dtype=int) + data.atomcoords = np.asarray(frames, dtype=float) + data.natom = len(atomnos) + return parser, data diff --git a/python/vmol/vmol2.py b/python/vmol/vmol2.py index 9cfb38d..11ed6f2 100755 --- a/python/vmol/vmol2.py +++ b/python/vmol/vmol2.py @@ -4,7 +4,10 @@ import sys import os import warnings + from vmol import vmol, _paths +from vmol.plugins import plugin_for +from vmol.plugins.orca import OrcaPlugin try: import cclib @@ -13,48 +16,273 @@ raise ModuleNotFoundError(msg) from None -def read_mols_from_cclib(path): - """Read molecular data from a file using cclib. +VIEW_KEYS = ("freq", "opt", "sp", "irc", "scan", "geom") - Arguments: - path (str): Path to the file to read. +_CCLIB2VMOL = { + "vibdisps": "disp", + "vibfreqs": "freq", + "vibrmasses": "mass", + "vibirs": "ints", +} + + +def _vib_from_data(data): + vib = {dst: getattr(data, src, None) for src, dst in _CCLIB2VMOL.items()} + if vib["disp"] is None or vib["freq"] is None: + return None + return vib + + +def _mols_from_data(data, parser): + return [{"q": data.atomnos, "r": r, "name": str(parser)} for r in data.atomcoords] + + +def _peek_parser(path): + """Open *path* with cclib without parsing contents. Returns: - tuple containing: - list[dict]: Dictionaries containing atomic numbers ('q'), coordinates ('r'), - and a name for each molecule, i.e. suitable argument for vmol.run(). - dict or None: Dictionary containing vibrational frequencies ('freq') and modes ('disp'). + cclib parser instance. Raises: - FileNotFoundError: If the specified file is not found. - RuntimeError: If cclib cannot parse the file or if no coordinates are found in the parsed data. + FileNotFoundError: If *path* does not exist. + RuntimeError: If cclib cannot identify the file. """ if not os.path.isfile(path): msg = f"File {path} not found." raise FileNotFoundError(msg) - parser = cclib.io.ccopen(path) if parser is None: msg = f"Could not open {path} with cclib." raise RuntimeError(msg) + return parser + + +def _parse_file(path): + """Open *path* with cclib as ``(parser, data)``. + Returns: + tuple: ``(parser, data)``. + """ + parser = _peek_parser(path) try: data = parser.parse() - except Exception: # ruff: ignore[blind-except] + except Exception: # noqa: BLE001 msg = f"Error parsing {path} with cclib. Using partial data." warnings.warn(msg, RuntimeWarning, stacklevel=2) data = parser + return parser, data + +def read_mols_from_cclib(path): + """Read molecular data from a file using cclib. + + Arguments: + path (str): Path to the file to read. + + Returns: + tuple containing: + list[dict]: Dictionaries containing atomic numbers ('q'), coordinates ('r'), + and a name for each molecule, i.e. suitable argument for vmol.run(). + dict or None: Dictionary containing vibrational frequencies ('freq') and modes ('disp'). + + Raises: + RuntimeError: If cclib cannot parse the file or if no coordinates are found + in the parsed data. + """ + parser, data = _parse_file(path) if not (hasattr(data, "atomcoords") and hasattr(data, "atomnos") and len(data.atomcoords)): msg = f"No coordinates found in {path}." raise RuntimeError(msg) + return _mols_from_data(data, parser), _vib_from_data(data) - cclib2vmol = {"vibdisps":"disp", "vibfreqs":"freq", "vibrmasses":"mass", "vibirs":"ints"} - vib = {key2: getattr(data, key1, None) for key1, key2 in cclib2vmol.items()} - if vib["disp"] is None or vib["freq"] is None: - vib = None - return [{'q': data.atomnos, 'r': r, 'name': str(parser)} for r in data.atomcoords], vib +def _collect_view_selectors(argv): + """Parse ``view:`` / ``vib:`` tokens from *argv*. + + Returns: + tuple: ``(view_keys, vib_mode, indices_to_pop)``. + + Raises: + ValueError: If a ``view:`` key is unknown. + """ + view_keys = [] + vib_mode = None + to_pop = [] + for i, arg in enumerate(argv): + if arg.startswith("view:"): + key = arg.split(":", 1)[1].strip().lower() + if key not in VIEW_KEYS: + msg = ( + f"Unknown view selector '{arg}'. Expected one of: " + + ", ".join(f"view:{k}" for k in VIEW_KEYS) + ) + raise ValueError(msg) + view_keys.append(key) + to_pop.append(i) + elif arg.startswith("vib:"): + try: + vib_mode = bool(int(arg.split(":", 1)[1])) + except ValueError: + continue + return view_keys, vib_mode, to_pop + + +def _default_datasets(datasets): + for ds in datasets: + if ds.key == "freq": + return [ds] + return [datasets[0]] + + +def _prompt_datasets(datasets): + print(f"Found {len(datasets)} dataset(s) in this file:") + for i, ds in enumerate(datasets, start=1): + print(f" {i}) [{ds.key}] {ds.label}") + print("Enter number(s), comma-separated list, or 'all': ", end="", flush=True) + try: + line = sys.stdin.readline() + except EOFError: + line = "" + + choice = line.strip().lower() + if not choice: + return _default_datasets(datasets) + if choice in {"all", "*"}: + return list(datasets) + + selected = [] + for part in choice.split(","): + part = part.strip() + if not part: + continue + try: + idx = int(part) + except ValueError: + warnings.warn(f"Ignoring invalid choice '{part}'.", RuntimeWarning, stacklevel=2) + continue + if not 1 <= idx <= len(datasets): + warnings.warn(f"Ignoring out-of-range choice {idx}.", RuntimeWarning, stacklevel=2) + continue + selected.append(datasets[idx - 1]) + return selected or _default_datasets(datasets) + + +def _select_by_keys(datasets, keys): + by_key = {ds.key: ds for ds in datasets} + selected = [] + for key in keys: + if key not in by_key: + warnings.warn( + f"Requested view:{key} but that dataset is not available " + f"(have: {', '.join(by_key)}).", + RuntimeWarning, + stacklevel=2, + ) + continue + selected.append(by_key[key]) + return selected + + +def _select_by_vib(datasets, vib_mode): + if vib_mode: + for ds in datasets: + if ds.key == "freq": + return [ds] + warnings.warn("vib:1 requested but no frequencies found.", RuntimeWarning, stacklevel=2) + return _default_datasets(datasets) + + for key in ("opt", "sp", "irc", "scan", "geom"): + for ds in datasets: + if ds.key == key: + return [ds] + return _default_datasets(datasets) + + +def _choose_datasets(datasets, view_keys, vib_mode): + if len(datasets) == 1: + return list(datasets) + if view_keys: + return _select_by_keys(datasets, view_keys) or _default_datasets(datasets) + if vib_mode is not None: + return _select_by_vib(datasets, vib_mode) + if sys.stdin.isatty(): + return _prompt_datasets(datasets) + return _default_datasets(datasets) + + +def _run_datasets(datasets, args): + ret = 0 + for ds in datasets: + ret = vmol.run(mols=ds.mols, vib=ds.vib, args=args, with_arg0=True) + return ret + + +def _handle_orca(path, view_keys, vib_mode, args): + _parser, data = OrcaPlugin.parse(path) + datasets = OrcaPlugin.inventory(path, data) + chosen = _choose_datasets(datasets, view_keys, vib_mode) + if not chosen: + msg = f"No datasets selected for {path}." + raise RuntimeError(msg) + return _run_datasets(chosen, args) + + +def _load_cli_files(view_keys, vib_mode): + """Parse CLI molecule arguments into a viewer payload. + + Returns: + tuple: ``(to_pop, mols, vib, last_ret, orca_handled)``. + """ + to_pop = [] + mols = [] + vib = None + last_ret = 0 + orca_handled = False + + for i, arg in enumerate(sys.argv[1:], start=1): + if ":" in arg: + continue + to_pop.append(i) + popped = set(to_pop) + + try: + peek = _peek_parser(arg) + except (FileNotFoundError, RuntimeError) as e: + warnings.warn(str(e), RuntimeWarning, stacklevel=2) + continue + + if plugin_for(peek) is not None: + remaining = [sys.argv[j] for j in range(len(sys.argv)) if j not in popped] + try: + last_ret = _handle_orca(arg, view_keys, vib_mode, remaining) + except RuntimeError as e: + warnings.warn(str(e), RuntimeWarning, stacklevel=2) + continue + orca_handled = True + continue + + try: + parser, data = _parse_file(arg) + except (FileNotFoundError, RuntimeError) as e: + warnings.warn(str(e), RuntimeWarning, stacklevel=2) + continue + + if not ( + hasattr(data, "atomcoords") + and hasattr(data, "atomnos") + and len(data.atomcoords) + ): + warnings.warn(f"No coordinates found in {arg}.", RuntimeWarning, stacklevel=2) + continue + + m = _mols_from_data(data, parser) + v = _vib_from_data(data) + if not mols and v is not None and vib_mode is not False: + mols, vib = m, v + break + mols.extend(m) + + return to_pop, mols, vib, last_ret, orca_handled def main(): @@ -63,6 +291,10 @@ def main(): Usage: vmol2.py file [options] + ORCA outputs may offer frequencies, optimization, single point, scan, + and/or IRC. With several datasets and a TTY, a menu is shown. Otherwise + use ``view:freq|opt|sp|irc|scan`` or ``vib:0`` / ``vib:1``. + Returns: int: Exit code (0 for success). @@ -71,32 +303,25 @@ def main(): RuntimeError: If no valid molecular data is found in any provided files. """ if vmol is None: - paths = '\n'.join(_paths) - msg = f"vmol shared library not found in paths. Check your installation.\nSearched paths:\n{paths}" + paths = "\n".join(_paths) + msg = ( + f"vmol shared library not found in paths. Check your installation.\n" + f"Searched paths:\n{paths}" + ) raise ImportError(msg) - if len(sys.argv)==1: + if len(sys.argv) == 1: vmol.run(args=sys.argv, with_arg0=True) return 0 - to_pop = [] - mols = [] - vib = None - vib_disabled = 'vib:0' in sys.argv[1:] + view_keys, vib_mode, view_pops = _collect_view_selectors(sys.argv[1:]) + for i in reversed([j + 1 for j in view_pops]): + sys.argv.pop(i) - for i, arg in enumerate(sys.argv[1:], start=1): - if ':' in arg: - continue - to_pop.append(i) - try: - m, v = read_mols_from_cclib(arg) - except (FileNotFoundError, RuntimeError) as e: - m, v = [], None - warnings.warn(str(e), RuntimeWarning, stacklevel=2) - if len(mols)==0 and v is not None and not vib_disabled: - mols, vib = m, v - break - mols.extend(m) + to_pop, mols, vib, last_ret, orca_handled = _load_cli_files(view_keys, vib_mode) + + if orca_handled and not mols: + return last_ret if not mols: msg = "No valid molecular data found in any provided files." @@ -104,7 +329,6 @@ def main(): for i in reversed(to_pop): sys.argv.pop(i) - return vmol.run(mols=mols, vib=vib, args=sys.argv, with_arg0=True) From dc0a52671cd47fec402f3b60848ebd937d7648de Mon Sep 17 00:00:00 2001 From: rlaplaza Date: Thu, 15 Oct 2026 11:25:32 +0200 Subject: [PATCH 2/5] Add Gaussian vmol2 plugin and share inventory with ORCA Standardize both packages around a common parse/inventory layer so opt/freq/IRC/scan selectors behave the same across log formats. Co-authored-by: Cursor --- .github/CHANGELOG.md | 2 +- README.md | 6 +- python/README.md | 29 +-- python/ruff.toml | 2 +- python/tests/_plugin_fixtures.py | 20 ++ python/tests/test_gaussian_plugin.py | 112 +++++++++ python/tests/test_orca_plugin.py | 62 ++--- python/vmol/plugins/__init__.py | 15 +- python/vmol/plugins/common.py | 341 +++++++++++++++++++++++++++ python/vmol/plugins/gaussian.py | 55 +++++ python/vmol/plugins/orca.py | 281 +++------------------- python/vmol/plugins/orca_parse.py | 50 ++-- python/vmol/vmol2.py | 93 ++++---- 13 files changed, 664 insertions(+), 404 deletions(-) create mode 100644 python/tests/_plugin_fixtures.py create mode 100644 python/tests/test_gaussian_plugin.py create mode 100644 python/vmol/plugins/common.py create mode 100644 python/vmol/plugins/gaussian.py diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index cb5f061..0dfc9f2 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -37,7 +37,7 @@ ### Coming in the next version: * high-symmetry determination bugs (#21) * headless fallback build by @JamesOBrien2 (#44) -* ORCA plugin for `vmol2`: prompt / `view:` selectors for frequencies, optimization, single point, scan, and IRC sidecar trajectories +* ORCA / Gaussian plugins for `vmol2`: prompt / `view:` selectors for frequencies, optimization, single point, scan, and IRC trajectories **Full Changelog**: https://github.com/briling/v/compare/v2.0...v3.0rc7 diff --git a/README.md b/README.md index 2bd89d0..fe9ba21 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ A simple X11 molecular viewer. - [`.xyz`](https://en.wikipedia.org/wiki/XYZ_file_format) files - [extended `.xyz`](https://github.com/libAtoms/extxyz) files (currently the extra columns are ignored) - various quantum-chemical outputs with [cclib](https://github.com/cclib/cclib), see the Python wrapper [page](python/README.md) - (`vmol2`; ORCA gets a plugin for frequencies, optimization, single point, scan, and IRC) + (`vmol2`; ORCA and Gaussian plugins for frequencies, optimization, single point, scan, and IRC) ## Python package (wrapper / API) available @@ -19,8 +19,8 @@ See python package page [here](python/README.md). Provides wrapper scripts with a simple installation and allows to open unsupported file formats with `cclib` (`vmol2`). -For ORCA multi-job logs, `vmol2` can prompt (or take `view:` / `vib:`) for -frequencies, optimization, single point, IRC, or scan +For ORCA / Gaussian multi-job logs, `vmol2` can prompt (or take `view:` / `vib:`) +for frequencies, optimization, single point, IRC, or scan (see [python/README.md](python/README.md)). ## Download diff --git a/python/README.md b/python/README.md index d760b6d..adfc0b3 100644 --- a/python/README.md +++ b/python/README.md @@ -46,32 +46,28 @@ vmol ../mol/MOL_3525.xyz cell:8.93,0.0,0.0,4.2,8.9,0.0,0.48,2.32,10 python -m vmol../mol/periodic.in bonds:0 ``` -Unsupported quantum-chemical outputs, e.g., ORCA: +Unsupported quantum-chemical outputs, e.g., ORCA or Gaussian: ```bash vmol2 ../mol/CEHZOF_1_SPE.out python -m vmol.vmol2 ../mol/AJALIH_5_SPE.out +vmol2 mol/Cat1_TSOA_CfA_Cf0.log ``` -ORCA logs use a dedicated plugin (resilient to ORCA 6 / PyPI cclib). When a -file has several of frequencies, optimization, trailing single point, scan, or -IRC, `vmol2` prompts on a TTY (each choice opens in its own `v` session), or -accepts: +ORCA and Gaussian logs use dedicated plugins. When a file has several of +frequencies, optimization, trailing single point (ORCA), scan, or IRC, +`vmol2` prompts on a TTY (each choice opens in its own `v` session), or accepts: ```bash vmol2 job.out view:freq # normal modes (on the opt minimum) vmol2 job.out view:opt # geometry optimization frames -vmol2 job.out view:sp # final single-point geometry / energy -vmol2 job.out view:irc # IRC (needs basename_IRC_Full_trj.xyz beside the .out) +vmol2 job.out view:sp # final single-point geometry / energy (ORCA multi-job) +vmol2 job.out view:irc # IRC trajectory vmol2 job.out view:scan # relaxed surface scan vmol2 job.out vib:0 # geometry path (opt, else SP, else IRC, else scan) vmol2 job.out vib:1 # frequencies ``` -IRC frames come from `basename_IRC_Full_trj.xyz` (or the path in the log's -`Storing full IRC trajectory in ....` line), not from the `.out` itself. - -Gaussian and other packages use the generic cclib path: -```bash -vmol2 mol/Cat1_TSOA_CfA_Cf0.log vib:0 # geometry optimization -vmol2 mol/Cat1_TSOA_CfA_Cf0.log # normal modes (default for files containing both) -``` +ORCA IRC frames come from `basename_IRC_Full_trj.xyz` (or the path in the log's +`Storing full IRC trajectory in ....` line). Gaussian IRC frames are read from +the `.log` itself when an IRC job is detected. Other packages still use the +generic cclib path. ### Library @@ -172,7 +168,8 @@ out = vmol.capture(mols=mols) File formats which are not supported natively can be read with `cclib` and passed (see [example](https://github.com/briling/v/tree/master/python/examples/ex_cclib.py), [vmol2](https://github.com/briling/v/tree/master/python/vmol/vmol2.py), and the -ORCA [plugin](https://github.com/briling/v/tree/master/python/vmol/plugins/orca.py)). +ORCA / Gaussian plugins under +[python/vmol/plugins](https://github.com/briling/v/tree/master/python/vmol/plugins)). One can also pass normal modes data with ```python diff --git a/python/ruff.toml b/python/ruff.toml index 19a6a94..3df1a61 100644 --- a/python/ruff.toml +++ b/python/ruff.toml @@ -53,4 +53,4 @@ ignore = [ [lint.per-file-ignores] "vmol/__main__.py" = ["D"] "examples/ex*.py" = ["D"] -"tests/*.py" = ["SLF", "PLC2701"] +"tests/*.py" = ["SLF", "PLC2701", "D"] diff --git a/python/tests/_plugin_fixtures.py b/python/tests/_plugin_fixtures.py new file mode 100644 index 0000000..8b743c2 --- /dev/null +++ b/python/tests/_plugin_fixtures.py @@ -0,0 +1,20 @@ +"""Shared fixtures for ORCA / Gaussian plugin tests.""" + +from __future__ import annotations + +from types import SimpleNamespace + +import numpy as np + + +def coords(*frames): + return np.asarray(frames, dtype=float) + + +def base_data(nframes=1): + frames = [[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74 + 0.01 * i]] for i in range(nframes)] + return SimpleNamespace( + atomnos=np.array([1, 1]), + atomcoords=coords(*frames), + scfenergies=np.array([-27.2 - 0.01 * i for i in range(nframes)]), + ) diff --git a/python/tests/test_gaussian_plugin.py b/python/tests/test_gaussian_plugin.py new file mode 100644 index 0000000..27431cd --- /dev/null +++ b/python/tests/test_gaussian_plugin.py @@ -0,0 +1,112 @@ +"""Unit tests for the Gaussian vmol2 plugin (no viewer / X11).""" + +from __future__ import annotations + +import os +import tempfile +import unittest +from pathlib import Path + +import numpy as np + +from _plugin_fixtures import base_data +from vmol import vmol2 +from vmol.plugins import plugin_for +from vmol.plugins.gaussian import GaussianPlugin, _log_has_irc + + +class TestGaussianInventory(unittest.TestCase): + """Inventory detection for Gaussian frequencies / optimization / IRC.""" + + def setUp(self): + """Create a temporary directory for fake log paths.""" + self.tmpdir = tempfile.TemporaryDirectory() + self.addCleanup(self.tmpdir.cleanup) + self.log_path = os.path.join(self.tmpdir.name, "job.log") + + def _write_log(self, text: str) -> str: + Path(self.log_path).write_text(text, encoding="utf-8") + return self.log_path + + def test_opt_plus_freq(self): + """Opt+freq Link1 jobs expose both datasets.""" + data = base_data(4) + data.geovalues = np.ones((4, 4)) + data.optdone = True + data.vibfreqs = np.array([100.0, 200.0]) + data.vibdisps = np.zeros((2, 2, 3)) + data.vibdisps[:, 0, 0] = 1.0 + self._write_log("Entering Gaussian System\n") + + datasets = GaussianPlugin.inventory(self.log_path, data) + self.assertEqual([d.key for d in datasets], ["freq", "opt"]) + self.assertEqual(len(datasets[1].mols), 4) + np.testing.assert_allclose(datasets[0].mols[0]["r"], data.atomcoords[3]) + + def test_irc_log_labeled_as_irc(self): + """IRC marker lines label the trajectory as irc, not opt.""" + data = base_data(5) + data.geovalues = np.ones((5, 4)) + data.optdone = True + self._write_log( + "Entering Gaussian System\nIRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC\n", + ) + datasets = GaussianPlugin.inventory(self.log_path, data) + self.assertEqual([d.key for d in datasets], ["irc"]) + self.assertEqual(len(datasets[0].mols), 5) + + def test_single_point_geometry(self): + """A lone single-point job yields one geom dataset.""" + data = base_data(1) + self._write_log("Entering Gaussian System\n") + datasets = GaussianPlugin.inventory(self.log_path, data) + self.assertEqual([d.key for d in datasets], ["geom"]) + + def test_log_has_irc_route(self): + """Route lines containing irc are detected.""" + path = self._write_log("#p b3lyp/6-31g(d) irc=(calcfc,maxpoints=20)\n") + self.assertTrue(_log_has_irc(path)) + path = self._write_log("#p b3lyp/6-31g(d) opt freq\n") + self.assertFalse(_log_has_irc(path)) + + +class TestGaussianSampleFile(unittest.TestCase): + """Integration check against the repo Gaussian opt+freq sample.""" + + @classmethod + def setUpClass(cls): + """Resolve the sample Gaussian logfile paths.""" + mol = os.path.abspath( + os.path.join(os.path.dirname(__file__), "..", "..", "mol"), + ) + cls.path = os.path.join(mol, "Cat1_TSOA_CfA_Cf0.log") + cls.spc = os.path.join(mol, "Cat1_TSOA_CfA_Cf0_SPC.log") + + def test_parse_and_inventory_opt_freq(self): + """Sample opt+freq log exposes freq and opt selectors.""" + if not os.path.isfile(self.path): + self.skipTest("sample Gaussian logfile not present") + + parser, data = GaussianPlugin.parse(self.path) + self.assertEqual(parser.logname, "Gaussian") + self.assertIs(plugin_for(parser), GaussianPlugin) + datasets = GaussianPlugin.inventory(self.path, data) + self.assertEqual([d.key for d in datasets], ["freq", "opt"]) + by_key = {d.key: d for d in datasets} + self.assertEqual(len(by_key["freq"].vib["freq"]), 123) + self.assertEqual(len(by_key["opt"].mols), 10) + chosen = vmol2._choose_datasets(datasets, ["opt"], None) + self.assertEqual(chosen[0].key, "opt") + + def test_parse_spc(self): + """Separate single-point log is a single geometry.""" + if not os.path.isfile(self.spc): + self.skipTest("sample Gaussian SPC logfile not present") + _parser, data = GaussianPlugin.parse(self.spc) + datasets = GaussianPlugin.inventory(self.spc, data) + self.assertEqual([d.key for d in datasets], ["geom"]) + self.assertEqual(len(datasets[0].mols), 1) + + +if __name__ == "__main__": + unittest.main() diff --git a/python/tests/test_orca_plugin.py b/python/tests/test_orca_plugin.py index e1c55e2..7702745 100644 --- a/python/tests/test_orca_plugin.py +++ b/python/tests/test_orca_plugin.py @@ -7,33 +7,20 @@ import unittest import warnings from pathlib import Path -from types import SimpleNamespace import numpy as np +from _plugin_fixtures import base_data from vmol import vmol2 -from vmol.plugins.orca import Dataset, OrcaPlugin, load_irc_xyz, resolve_irc_path - - -def _coords(*frames): - return np.asarray(frames, dtype=float) - - -def _base_data(nframes=1): - frames = [[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74 + 0.01 * i]] for i in range(nframes)] - return SimpleNamespace( - atomnos=np.array([1, 1]), - atomcoords=_coords(*frames), - scfenergies=np.array([-27.2 - 0.01 * i for i in range(nframes)]), - ) +from vmol.plugins.common import Dataset +from vmol.plugins.orca import OrcaPlugin, load_irc_xyz, resolve_irc_path class TestOrcaInventory(unittest.TestCase): """Inventory detection for frequencies, optimization, scan, and IRC.""" def setUp(self): - """Create a temporary directory and plugin instance.""" - self.plugin = OrcaPlugin() + """Create a temporary directory for fake ORCA paths.""" self.tmpdir = tempfile.TemporaryDirectory() self.addCleanup(self.tmpdir.cleanup) self.out_path = os.path.join(self.tmpdir.name, "job.out") @@ -44,7 +31,7 @@ def _write_out(self, text: str) -> str: def test_frequency_only(self): """A frequency job yields a single freq dataset.""" - data = _base_data(1) + data = base_data(1) data.vibfreqs = np.array([4401.0]) data.vibdisps = np.zeros((1, 2, 3)) data.vibdisps[0, 0, 2] = 0.7 @@ -53,7 +40,7 @@ def test_frequency_only(self): data.vibrmasses = np.array([1.0]) self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["freq"]) self.assertEqual(len(datasets[0].mols), 1) self.assertIsNotNone(datasets[0].vib) @@ -61,7 +48,7 @@ def test_frequency_only(self): def test_opt_plus_freq(self): """Opt+freq offers both datasets and truncates at optdone.""" - data = _base_data(4) + data = base_data(4) data.geovalues = np.ones((3, 4)) data.optdone = [3] data.vibfreqs = np.array([100.0, 200.0]) @@ -69,7 +56,7 @@ def test_opt_plus_freq(self): data.vibdisps[:, 0, 0] = 1.0 self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["freq", "opt"]) opt = datasets[1] self.assertEqual(len(opt.mols), 3) @@ -78,18 +65,18 @@ def test_opt_plus_freq(self): def test_opt_truncates_when_optdone_bool(self): """Bool optdone uses geovalues (+ final eval), dropping a later job frame.""" - data = _base_data(5) + data = base_data(5) data.geovalues = np.ones((3, 4)) data.optdone = True self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["opt"]) self.assertEqual(len(datasets[0].mols), 4) def test_opt_freq_and_sp(self): """Multi-job files expose frequencies, opt, and a selectable SP.""" - data = _base_data(5) + data = base_data(5) data.geovalues = np.ones((3, 4)) data.optdone = True data.vibfreqs = np.array([100.0]) @@ -99,7 +86,7 @@ def test_opt_freq_and_sp(self): data.sp_energy_eh = -100.5 self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["freq", "opt", "sp"]) self.assertEqual(len(datasets[1].mols), 4) self.assertIn("E=-100.50000000 Eh", datasets[2].label) @@ -107,27 +94,27 @@ def test_opt_freq_and_sp(self): def test_scan_not_labeled_as_opt(self): """Relaxed scans are labeled scan, not opt.""" - data = _base_data(3) + data = base_data(3) data.scannames = ["R"] data.scanparm = [(0.7,), (0.8,), (0.9,)] data.geovalues = np.ones((3, 4)) self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["scan"]) self.assertIn("R=0.7", datasets[0].mols[0]["name"]) def test_single_point_geometry(self): """A single-point job yields one geom dataset.""" - data = _base_data(1) + data = base_data(1) self._write_out("* O R C A *\n") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["geom"]) self.assertEqual(len(datasets[0].mols), 1) def test_missing_irc_file_warns(self): """Mentioned but missing IRC xyz warns and is omitted.""" - data = _base_data(1) + data = base_data(1) self._write_out( "* O R C A *\n" "Intrinsic Reaction Coordinate Calculation\n" @@ -135,13 +122,13 @@ def test_missing_irc_file_warns(self): ) with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always") - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["geom"]) self.assertTrue(any("IRC trajectory" in str(w.message) for w in caught)) def test_irc_sidecar_loaded(self): """IRC xyz comments become frame names via cclib.""" - data = _base_data(1) + data = base_data(1) xyz_path = os.path.join(self.tmpdir.name, "job_IRC_Full_trj.xyz") Path(xyz_path).write_text( "2\n" @@ -158,7 +145,7 @@ def test_irc_sidecar_loaded(self): "* O R C A *\n" "Storing full IRC trajectory in .... job_IRC_Full_trj.xyz\n", ) - datasets = self.plugin.inventory(self.out_path, data) + datasets = OrcaPlugin.inventory(self.out_path, data) self.assertEqual([d.key for d in datasets], ["irc"]) irc = datasets[0] self.assertEqual(len(irc.mols), 2) @@ -196,16 +183,15 @@ def test_load_irc_xyz_comments(self): with tempfile.TemporaryDirectory() as tmp: xyz_path = os.path.join(tmp, "trj.xyz") Path(xyz_path).write_text( - "2\n" - "Coordinates from ORCA-job trj E -99.5\n" - "H 0 0 0\n" - "H 0 0 1\n", + "2\nCoordinates from ORCA-job trj E -99.5\nH 0 0 0\nH 0 0 1\n", encoding="utf-8", ) ds = load_irc_xyz(xyz_path) self.assertEqual(ds.key, "irc") self.assertEqual(len(ds.mols), 1) - self.assertEqual(ds.mols[0]["name"], "Coordinates from ORCA-job trj E -99.5") + self.assertEqual( + ds.mols[0]["name"], "Coordinates from ORCA-job trj E -99.5", + ) class TestViewSelection(unittest.TestCase): diff --git a/python/vmol/plugins/__init__.py b/python/vmol/plugins/__init__.py index 7c8f7a4..95470ed 100644 --- a/python/vmol/plugins/__init__.py +++ b/python/vmol/plugins/__init__.py @@ -1,12 +1,17 @@ """Format-specific plugins for reading quantum-chemical outputs.""" +from vmol.plugins.gaussian import GaussianPlugin from vmol.plugins.orca import OrcaPlugin -__all__ = ["OrcaPlugin", "plugin_for"] +__all__ = ["GaussianPlugin", "OrcaPlugin", "plugin_for"] + +_PLUGINS = { + "ORCA": OrcaPlugin, + "Gaussian": GaussianPlugin, +} def plugin_for(parser): - """Return a plugin for *parser*, or None if no dedicated plugin exists.""" - if getattr(parser, "logname", None) == "ORCA": - return OrcaPlugin() - return None + """Return a plugin class for *parser*, or None if none is registered.""" + name = getattr(parser, "logname", None) + return _PLUGINS.get(name) diff --git a/python/vmol/plugins/common.py b/python/vmol/plugins/common.py new file mode 100644 index 0000000..b77d43a --- /dev/null +++ b/python/vmol/plugins/common.py @@ -0,0 +1,341 @@ +"""Shared dataset inventory helpers for quantum-chemistry log plugins.""" + +from __future__ import annotations + +import os +import warnings +from collections.abc import Callable, Sequence +from dataclasses import dataclass + +import numpy as np +from cclib.io import ccopen +from cclib.parser import utils + +# On-screen molecule name buffer in v (src/v/v.h STRLEN). +_NAME_MAX = 255 + +CCLIB2VMOL = { + "vibdisps": "disp", + "vibfreqs": "freq", + "vibrmasses": "mass", + "vibirs": "ints", +} + + +@dataclass(frozen=True) +class Dataset: + """One viewable dataset extracted from a quantum-chemistry log.""" + + key: str + label: str + mols: list + vib: dict | None = None + + +def has_coords(data) -> bool: + """Check whether *data* has usable atomcoords and atomnos. + + Returns: + bool: True when coordinates are present. + """ + return ( + hasattr(data, "atomcoords") + and hasattr(data, "atomnos") + and len(data.atomcoords) > 0 + ) + + +def parse_cclib(path: str): + """Open and parse *path* with cclib. + + Returns: + tuple: ``(parser, data)``. On parse errors, *data* may be the partial + parser object (with a warning). + + Raises: + RuntimeError: If cclib cannot open the file. + """ + parser = ccopen(path) + if parser is None: + msg = f"Could not open {path} with cclib." + raise RuntimeError(msg) + try: + data = parser.parse() + except Exception as exc: # noqa: BLE001 + warnings.warn( + f"Error parsing {path} with cclib ({type(exc).__name__}: {exc}). " + "Using partial data.", + RuntimeWarning, + stacklevel=2, + ) + data = parser + return parser, data + + +def _clip_name(name: str) -> str: + if len(name) <= _NAME_MAX: + return name + return name[: _NAME_MAX - 3] + "..." + + +def _has_frequencies(data) -> bool: + freqs = getattr(data, "vibfreqs", None) + disps = getattr(data, "vibdisps", None) + return freqs is not None and disps is not None and len(freqs) and len(disps) + + +def vib_from_data(data) -> dict | None: + """Map cclib vibration attributes to a vmol vib dict. + + Returns: + dict | None: Vibration payload, or None if modes are missing. + """ + vib = {dst: getattr(data, src, None) for src, dst in CCLIB2VMOL.items()} + if vib["disp"] is None or vib["freq"] is None: + return None + return vib + + +def _energy_suffix(data, index: int) -> str: + energies = getattr(data, "scfenergies", None) + if energies is None or index >= len(energies): + return "" + val = float(energies[index]) + if np.isnan(val): + return "" + eh = utils.convertor(val, "eV", "hartree") + return f" E={eh:.8f} Eh" + + +def _frame_mols(data, coords, *, name_prefix: str, name_for_index=None) -> list: + atomnos = data.atomnos + n = len(coords) + mols = [] + for i, r in enumerate(coords): + if name_for_index is not None: + name = name_for_index(i, r) + else: + name = f"{name_prefix} step {i + 1}/{n}{_energy_suffix(data, i)}" + mols.append({"q": atomnos, "r": r, "name": _clip_name(name)}) + return mols + + +def _opt_end_index(data) -> int | None: + optdone = getattr(data, "optdone", None) + n = len(data.atomcoords) + if isinstance(optdone, (list, tuple)) and optdone: + last = max(int(x) for x in optdone) + if 1 <= last <= n: + return last + return None + if optdone is True: + geovalues = getattr(data, "geovalues", None) + if geovalues is not None and len(geovalues): + return min(len(geovalues) + 1, n) + return None + + +def _freq_geometry_index(data) -> int: + end = _opt_end_index(data) + if end is not None and end > 0: + return end - 1 + return len(data.atomcoords) - 1 + + +def _is_scan(data) -> bool: + return ( + getattr(data, "scanparm", None) is not None + or getattr(data, "scannames", None) is not None + ) + + +def _opt_coords(data): + if _is_scan(data): + return None + has_opt = ( + getattr(data, "geovalues", None) is not None + or getattr(data, "optdone", None) is not None + ) + if not has_opt or len(data.atomcoords) <= 1: + return None + end = _opt_end_index(data) + coords = data.atomcoords[:end] if end is not None else data.atomcoords + if len(coords) <= 1: + return None + return coords + + +def _frequencies_dataset(data, name: str) -> Dataset: + idx = _freq_geometry_index(data) + return Dataset( + key="freq", + label=f"Frequencies ({len(data.vibfreqs)} modes)", + mols=[ + { + "q": data.atomnos, + "r": data.atomcoords[idx], + "name": _clip_name(f"{name} (frequencies)"), + }, + ], + vib=vib_from_data(data), + ) + + +def _optimization_dataset(data, name: str) -> Dataset | None: + coords = _opt_coords(data) + if coords is None: + return None + mols = _frame_mols(data, coords, name_prefix=f"{name} opt") + return Dataset( + key="opt", + label=f"Optimization trajectory ({len(mols)} frames)", + mols=mols, + vib=None, + ) + + +def _scan_dataset(data, name: str) -> Dataset | None: + if not _is_scan(data) or not len(data.atomcoords): + return None + + scanparm = getattr(data, "scanparm", None) + scannames = getattr(data, "scannames", None) or [] + n = len(data.atomcoords) + + def name_for_index(i, _r): + base = f"{name} scan step {i + 1}/{n}" + if scanparm is not None and i < len(scanparm): + params = scanparm[i] + if scannames and len(scannames) == len(params): + parts = [f"{n_}={v}" for n_, v in zip(scannames, params, strict=False)] + else: + parts = [str(v) for v in params] + base += " " + ", ".join(parts) + return base + _energy_suffix(data, i) + + mols = _frame_mols( + data, + data.atomcoords, + name_prefix=f"{name} scan", + name_for_index=name_for_index, + ) + return Dataset( + key="scan", + label=f"Relaxed scan ({len(mols)} frames)", + mols=mols, + vib=None, + ) + + +def _single_point_dataset(data, name: str) -> Dataset | None: + idx = getattr(data, "sp_frame_index", None) + if idx is None: + return None + if _opt_coords(data) is None and not _has_frequencies(data): + return None + + energy = getattr(data, "sp_energy_eh", None) + if energy is None: + suffix = _energy_suffix(data, idx) + else: + suffix = f" E={float(energy):.8f} Eh" + label_energy = suffix.strip() or "final single point" + return Dataset( + key="sp", + label=f"Single point ({label_energy})", + mols=[ + { + "q": data.atomnos, + "r": data.atomcoords[idx], + "name": _clip_name(f"{name} (single point){suffix}"), + }, + ], + vib=None, + ) + + +def _geometry_dataset(data, name: str) -> Dataset: + mols = [ + {"q": data.atomnos, "r": r, "name": _clip_name(name)} for r in data.atomcoords + ] + n = len(mols) + return Dataset( + key="geom", + label=f"Geometry ({n} frame{'s' if n != 1 else ''})", + mols=mols, + vib=None, + ) + + +def _irc_dataset(data, name: str, *, source: str | None = None) -> Dataset | None: + if len(data.atomcoords) <= 1: + return None + n = len(data.atomcoords) + label = f"IRC trajectory ({n} frames" + if source: + label += f", {os.path.basename(source)}" + label += ")" + return Dataset( + key="irc", + label=label, + mols=_frame_mols(data, data.atomcoords, name_prefix=f"{name} IRC"), + vib=None, + ) + + +def require_coords(path: str, data) -> None: + """Ensure *data* has coordinates. + + Raises: + RuntimeError: If *data* has no coordinates. + """ + if not has_coords(data): + msg = f"No coordinates found in {path}." + raise RuntimeError(msg) + + +def build_inventory( + path: str, + data, + *, + name: str, + prefer_irc: bool = False, + extras: Sequence[Callable[[], Dataset | None]] = (), +) -> list[Dataset]: + """Assemble selectable datasets from parsed *data*. + + Returns: + list[Dataset]: Non-empty list of datasets. + """ + require_coords(path, data) + datasets: list[Dataset] = [] + + if _has_frequencies(data): + datasets.append(_frequencies_dataset(data, name)) + + if prefer_irc: + irc = _irc_dataset(data, name) + if irc is not None: + datasets.append(irc) + else: + opt = _optimization_dataset(data, name) + if opt is not None: + datasets.append(opt) + + sp = _single_point_dataset(data, name) + if sp is not None: + datasets.append(sp) + + scan = _scan_dataset(data, name) + if scan is not None: + datasets.append(scan) + + for builder in extras: + extra = builder() + if extra is not None: + datasets.append(extra) + + if not datasets: + datasets.append(_geometry_dataset(data, name)) + + return datasets diff --git a/python/vmol/plugins/gaussian.py b/python/vmol/plugins/gaussian.py new file mode 100644 index 0000000..61257d0 --- /dev/null +++ b/python/vmol/plugins/gaussian.py @@ -0,0 +1,55 @@ +"""Gaussian plugin for vmol2.""" + +from __future__ import annotations + +import os + +from vmol.plugins.common import Dataset, build_inventory, parse_cclib, require_coords + +_IRC_MARKERS = ( + "IRC-IRC-IRC", + "INPUT DATA FOR L123", + "IRC Pt", +) + + +def _log_has_irc(path: str) -> bool: + with open(path, encoding="utf-8", errors="replace") as handle: + for line in handle: + if any(marker in line for marker in _IRC_MARKERS): + return True + stripped = line.lstrip() + if stripped.startswith("#") and " irc" in f" {stripped.lower()}": + return True + return False + + +def parse_gaussian(path: str): + """Parse a Gaussian logfile with cclib. + + Returns: + tuple: ``(parser, data)``. + """ + parser, data = parse_cclib(path) + require_coords(path, data) + return parser, data + + +class GaussianPlugin: + """Frequencies, optimization, scan, and IRC datasets from a Gaussian log.""" + + parse = staticmethod(parse_gaussian) + + @staticmethod + def inventory(path: str, data) -> list[Dataset]: + """List viewable datasets for *path* / *data*. + + Returns: + list[Dataset]: Non-empty list of datasets. + """ + return build_inventory( + path, + data, + name=f"Gaussian {os.path.basename(path)}", + prefer_irc=_log_has_irc(path), + ) diff --git a/python/vmol/plugins/orca.py b/python/vmol/plugins/orca.py index 017af55..14867b8 100644 --- a/python/vmol/plugins/orca.py +++ b/python/vmol/plugins/orca.py @@ -1,19 +1,17 @@ -"""ORCA-specific dataset inventory for vmol2.""" +"""ORCA plugin for vmol2 (resilient parse + IRC xyz sidecar).""" from __future__ import annotations import os import re import warnings -from dataclasses import dataclass -import numpy as np from cclib.io.xyzreader import XYZ -from cclib.parser import utils +from vmol.plugins.common import Dataset, build_inventory from vmol.plugins.orca_parse import parse_orca -# On-screen molecule name buffer in v (src/v/v.h STRLEN). +# Matches common._NAME_MAX / v STRLEN. _NAME_MAX = 255 _IRC_STORE_RE = re.compile( @@ -28,79 +26,6 @@ re.IGNORECASE, ) -CCLIB2VMOL = { - "vibdisps": "disp", - "vibfreqs": "freq", - "vibrmasses": "mass", - "vibirs": "ints", -} - - -@dataclass(frozen=True) -class Dataset: - """One viewable dataset extracted from an ORCA output.""" - - key: str - label: str - mols: list - vib: dict | None = None - - -def _clip_name(name: str) -> str: - if len(name) <= _NAME_MAX: - return name - return name[: _NAME_MAX - 3] + "..." - - -def _has_frequencies(data) -> bool: - freqs = getattr(data, "vibfreqs", None) - disps = getattr(data, "vibdisps", None) - return freqs is not None and disps is not None and len(freqs) and len(disps) - - -def _vib_from_data(data) -> dict: - return {dst: getattr(data, src, None) for src, dst in CCLIB2VMOL.items()} - - -def _energy_suffix(data, index: int) -> str: - energies = getattr(data, "scfenergies", None) - if energies is None or index >= len(energies): - return "" - val = float(energies[index]) - if np.isnan(val): - return "" - eh = utils.convertor(val, "eV", "hartree") - return f" E={eh:.8f} Eh" - - -def _frame_mols(data, coords, *, name_prefix: str, name_for_index=None) -> list: - atomnos = data.atomnos - n = len(coords) - mols = [] - for i, r in enumerate(coords): - if name_for_index is not None: - name = name_for_index(i, r) - else: - name = f"{name_prefix} step {i + 1}/{n}{_energy_suffix(data, i)}" - mols.append({"q": atomnos, "r": r, "name": _clip_name(name)}) - return mols - - -def _opt_end_index(data) -> int | None: - optdone = getattr(data, "optdone", None) - n = len(data.atomcoords) - if isinstance(optdone, (list, tuple)) and optdone: - last = max(int(x) for x in optdone) - if 1 <= last <= n: - return last - return None - if optdone is True: - geovalues = getattr(data, "geovalues", None) - if geovalues is not None and len(geovalues): - end = len(geovalues) + 1 - return min(end, n) - return None - def _scan_irc_filename(out_path: str) -> tuple[str | None, bool]: stored = None @@ -161,150 +86,40 @@ def load_irc_xyz(xyz_path: str) -> Dataset: data = XYZ(xyz_path).parse() comments = data.metadata.get("comments", []) nframes = len(data.atomcoords) - - def name_for_index(i, _r): + mols = [] + for i, r in enumerate(data.atomcoords): if i < len(comments) and comments[i].strip(): - return comments[i].strip() - return f"IRC frame {i + 1}/{nframes}" - + name = comments[i].strip() + else: + name = f"IRC frame {i + 1}/{nframes}" + if len(name) > _NAME_MAX: + name = name[: _NAME_MAX - 3] + "..." + mols.append({"q": data.atomnos, "r": r, "name": name}) return Dataset( key="irc", label=f"IRC trajectory ({nframes} frames, {os.path.basename(xyz_path)})", - mols=_frame_mols( - data, data.atomcoords, name_prefix="IRC", name_for_index=name_for_index, - ), - vib=None, - ) - - -def _freq_geometry_index(data) -> int: - end = _opt_end_index(data) - if end is not None and end > 0: - return end - 1 - return len(data.atomcoords) - 1 - - -def _frequencies_dataset(data, name: str) -> Dataset: - idx = _freq_geometry_index(data) - return Dataset( - key="freq", - label=f"Frequencies ({len(data.vibfreqs)} modes)", - mols=[{ - "q": data.atomnos, - "r": data.atomcoords[idx], - "name": _clip_name(f"{name} (frequencies)"), - }], - vib=_vib_from_data(data), - ) - - -def _is_scan(data) -> bool: - return ( - getattr(data, "scanparm", None) is not None - or getattr(data, "scannames", None) is not None - ) - - -def _opt_coords(data): - if _is_scan(data): - return None - has_opt = ( - getattr(data, "geovalues", None) is not None - or getattr(data, "optdone", None) is not None - ) - if not has_opt or len(data.atomcoords) <= 1: - return None - end = _opt_end_index(data) - coords = data.atomcoords[:end] if end is not None else data.atomcoords - if len(coords) <= 1: - return None - return coords - - -def _optimization_dataset(data, name: str) -> Dataset | None: - coords = _opt_coords(data) - if coords is None: - return None - mols = _frame_mols(data, coords, name_prefix=f"{name} opt") - return Dataset( - key="opt", - label=f"Optimization trajectory ({len(mols)} frames)", - mols=mols, - vib=None, - ) - - -def _scan_dataset(data, name: str) -> Dataset | None: - if not _is_scan(data) or not len(data.atomcoords): - return None - - scanparm = getattr(data, "scanparm", None) - scannames = getattr(data, "scannames", None) or [] - n = len(data.atomcoords) - - def name_for_index(i, _r): - base = f"{name} scan step {i + 1}/{n}" - if scanparm is not None and i < len(scanparm): - params = scanparm[i] - if scannames and len(scannames) == len(params): - parts = [f"{n_}={v}" for n_, v in zip(scannames, params, strict=False)] - else: - parts = [str(v) for v in params] - base += " " + ", ".join(parts) - return base + _energy_suffix(data, i) - - mols = _frame_mols( - data, data.atomcoords, name_prefix=f"{name} scan", name_for_index=name_for_index, - ) - return Dataset( - key="scan", - label=f"Relaxed scan ({len(mols)} frames)", mols=mols, vib=None, ) -def _single_point_dataset(data, name: str) -> Dataset | None: - idx = getattr(data, "sp_frame_index", None) - if idx is None: +def _load_irc_extra(path: str) -> Dataset | None: + irc_path = resolve_irc_path(path) + if irc_path is None: return None - if _opt_coords(data) is None and not _has_frequencies(data): + try: + return load_irc_xyz(irc_path) + except (OSError, ValueError, AssertionError) as exc: + warnings.warn( + f"Failed to load IRC trajectory {irc_path}: {exc}", + RuntimeWarning, + stacklevel=3, + ) return None - energy = getattr(data, "sp_energy_eh", None) - if energy is None: - suffix = _energy_suffix(data, idx) - else: - suffix = f" E={float(energy):.8f} Eh" - label_energy = suffix.strip() or "final single point" - return Dataset( - key="sp", - label=f"Single point ({label_energy})", - mols=[{ - "q": data.atomnos, - "r": data.atomcoords[idx], - "name": _clip_name(f"{name} (single point){suffix}"), - }], - vib=None, - ) - - -def _geometry_dataset(data, name: str) -> Dataset: - mols = [ - {"q": data.atomnos, "r": r, "name": _clip_name(name)} - for r in data.atomcoords - ] - n = len(mols) - return Dataset( - key="geom", - label=f"Geometry ({n} frame{'s' if n != 1 else ''})", - mols=mols, - vib=None, - ) - class OrcaPlugin: - """Detect frequencies, optimization, SP, scan, and IRC datasets in an ORCA log.""" + """Frequencies, optimization, SP, scan, and IRC datasets from an ORCA log.""" parse = staticmethod(parse_orca) @@ -314,48 +129,10 @@ def inventory(path: str, data) -> list[Dataset]: Returns: list[Dataset]: Non-empty list of datasets. - - Raises: - RuntimeError: If no atomic coordinates are present in *data*. """ - if not ( - hasattr(data, "atomcoords") - and hasattr(data, "atomnos") - and len(data.atomcoords) - ): - msg = f"No coordinates found in {path}." - raise RuntimeError(msg) - - name = f"ORCA {os.path.basename(path)}" - datasets: list[Dataset] = [] - - if _has_frequencies(data): - datasets.append(_frequencies_dataset(data, name)) - - opt = _optimization_dataset(data, name) - if opt is not None: - datasets.append(opt) - - sp = _single_point_dataset(data, name) - if sp is not None: - datasets.append(sp) - - scan = _scan_dataset(data, name) - if scan is not None: - datasets.append(scan) - - irc_path = resolve_irc_path(path) - if irc_path is not None: - try: - datasets.append(load_irc_xyz(irc_path)) - except (OSError, ValueError, AssertionError) as exc: - warnings.warn( - f"Failed to load IRC trajectory {irc_path}: {exc}", - RuntimeWarning, - stacklevel=2, - ) - - if not datasets: - datasets.append(_geometry_dataset(data, name)) - - return datasets + return build_inventory( + path, + data, + name=f"ORCA {os.path.basename(path)}", + extras=[lambda: _load_irc_extra(path)], + ) diff --git a/python/vmol/plugins/orca_parse.py b/python/vmol/plugins/orca_parse.py index a2636f1..6ec9dbc 100644 --- a/python/vmol/plugins/orca_parse.py +++ b/python/vmol/plugins/orca_parse.py @@ -7,10 +7,11 @@ from types import SimpleNamespace import numpy as np -from cclib.io import ccopen from cclib.parser import orcaparser, utils from cclib.parser.utils import PeriodicTable +from vmol.plugins.common import has_coords, parse_cclib + _PT = PeriodicTable() _PATCHED = False _CART_MIN_TOKENS = 4 @@ -61,7 +62,9 @@ def safe_extract(self, inputfile, line): _PATCHED = True -def _read_cartesian_block(lines, start: int) -> tuple[list[int], list[list[float]], int]: +def _read_cartesian_block( + lines, start: int, +) -> tuple[list[int], list[list[float]], int]: i = start + 1 if i < len(lines) and set(lines[i].strip()) <= {"-"}: i += 1 @@ -83,28 +86,22 @@ def _read_cartesian_block(lines, start: int) -> tuple[list[int], list[list[float return atomnos, coords, i +def _new_job(index: int) -> dict: + return {"index": index, "kind": "unknown", "frames": [], "energy": None} + + def _scan_jobs(path: str) -> list[dict]: with open(path, encoding="utf-8", errors="replace") as handle: lines = handle.readlines() jobs: list[dict] = [] - current = { - "index": 0, - "kind": "unknown", - "frames": [], - "energy": None, - } + current = _new_job(0) def flush(): nonlocal current if current["frames"] or current["kind"] != "unknown": jobs.append(current) - current = { - "index": len(jobs), - "kind": "unknown", - "frames": [], - "energy": None, - } + current = _new_job(len(jobs)) i = 0 while i < len(lines): @@ -159,7 +156,6 @@ def _as_namespace(data) -> SimpleNamespace: def _append_sp_from_jobs(ns: SimpleNamespace, jobs: list[dict]) -> None: - """Ensure the last single-point job geometry is available on *ns*.""" sp_jobs = [j for j in jobs if j["kind"] == "sp" and j["frames"]] if not sp_jobs: return @@ -197,38 +193,20 @@ def parse_orca(path: str): """Parse an ORCA logfile for visualization. Returns: - tuple: ``(parser, data)`` where *data* has at least atomcoords/atomnos - when the file contains geometries. Frequencies are included when present. + tuple: ``(parser, data)`` with atomcoords/atomnos when geometries exist. ``data.sp_frame_index`` is set when a trailing single-point job exists. Raises: RuntimeError: If the file cannot be opened or has no coordinates. """ _patch_cclib_orca() - - parser = ccopen(path) - if parser is None: - msg = f"Could not open {path} with cclib." - raise RuntimeError(msg) - - try: - raw = parser.parse() - except Exception as exc: # noqa: BLE001 - warnings.warn( - f"Error parsing {path} with cclib ({type(exc).__name__}: {exc}). " - "Using partial data.", - RuntimeWarning, - stacklevel=2, - ) - raw = parser - + parser, raw = parse_cclib(path) data = _as_namespace(raw) jobs = _scan_jobs(path) data.orca_jobs = jobs _append_sp_from_jobs(data, jobs) - coords = getattr(data, "atomcoords", None) - if coords is not None and len(coords) > 0: + if has_coords(data): return parser, data frames = [] diff --git a/python/vmol/vmol2.py b/python/vmol/vmol2.py index 11ed6f2..c3bb8fe 100755 --- a/python/vmol/vmol2.py +++ b/python/vmol/vmol2.py @@ -7,7 +7,7 @@ from vmol import vmol, _paths from vmol.plugins import plugin_for -from vmol.plugins.orca import OrcaPlugin +from vmol.plugins.common import has_coords, parse_cclib, vib_from_data try: import cclib @@ -18,20 +18,6 @@ VIEW_KEYS = ("freq", "opt", "sp", "irc", "scan", "geom") -_CCLIB2VMOL = { - "vibdisps": "disp", - "vibfreqs": "freq", - "vibrmasses": "mass", - "vibirs": "ints", -} - - -def _vib_from_data(data): - vib = {dst: getattr(data, src, None) for src, dst in _CCLIB2VMOL.items()} - if vib["disp"] is None or vib["freq"] is None: - return None - return vib - def _mols_from_data(data, parser): return [{"q": data.atomnos, "r": r, "name": str(parser)} for r in data.atomcoords] @@ -62,15 +48,14 @@ def _parse_file(path): Returns: tuple: ``(parser, data)``. + + Raises: + FileNotFoundError: If *path* does not exist. """ - parser = _peek_parser(path) - try: - data = parser.parse() - except Exception: # noqa: BLE001 - msg = f"Error parsing {path} with cclib. Using partial data." - warnings.warn(msg, RuntimeWarning, stacklevel=2) - data = parser - return parser, data + if not os.path.isfile(path): + msg = f"File {path} not found." + raise FileNotFoundError(msg) + return parse_cclib(path) def read_mols_from_cclib(path): @@ -90,10 +75,10 @@ def read_mols_from_cclib(path): in the parsed data. """ parser, data = _parse_file(path) - if not (hasattr(data, "atomcoords") and hasattr(data, "atomnos") and len(data.atomcoords)): + if not has_coords(data): msg = f"No coordinates found in {path}." raise RuntimeError(msg) - return _mols_from_data(data, parser), _vib_from_data(data) + return _mols_from_data(data, parser), vib_from_data(data) def _collect_view_selectors(argv): @@ -112,9 +97,8 @@ def _collect_view_selectors(argv): if arg.startswith("view:"): key = arg.split(":", 1)[1].strip().lower() if key not in VIEW_KEYS: - msg = ( - f"Unknown view selector '{arg}'. Expected one of: " - + ", ".join(f"view:{k}" for k in VIEW_KEYS) + msg = f"Unknown view selector '{arg}'. Expected one of: " + ", ".join( + f"view:{k}" for k in VIEW_KEYS ) raise ValueError(msg) view_keys.append(key) @@ -158,10 +142,14 @@ def _prompt_datasets(datasets): try: idx = int(part) except ValueError: - warnings.warn(f"Ignoring invalid choice '{part}'.", RuntimeWarning, stacklevel=2) + warnings.warn( + f"Ignoring invalid choice '{part}'.", RuntimeWarning, stacklevel=2, + ) continue if not 1 <= idx <= len(datasets): - warnings.warn(f"Ignoring out-of-range choice {idx}.", RuntimeWarning, stacklevel=2) + warnings.warn( + f"Ignoring out-of-range choice {idx}.", RuntimeWarning, stacklevel=2, + ) continue selected.append(datasets[idx - 1]) return selected or _default_datasets(datasets) @@ -188,7 +176,9 @@ def _select_by_vib(datasets, vib_mode): for ds in datasets: if ds.key == "freq": return [ds] - warnings.warn("vib:1 requested but no frequencies found.", RuntimeWarning, stacklevel=2) + warnings.warn( + "vib:1 requested but no frequencies found.", RuntimeWarning, stacklevel=2, + ) return _default_datasets(datasets) for key in ("opt", "sp", "irc", "scan", "geom"): @@ -217,9 +207,9 @@ def _run_datasets(datasets, args): return ret -def _handle_orca(path, view_keys, vib_mode, args): - _parser, data = OrcaPlugin.parse(path) - datasets = OrcaPlugin.inventory(path, data) +def _handle_plugin(plugin, path, view_keys, vib_mode, args): + _parser, data = plugin.parse(path) + datasets = plugin.inventory(path, data) chosen = _choose_datasets(datasets, view_keys, vib_mode) if not chosen: msg = f"No datasets selected for {path}." @@ -231,13 +221,13 @@ def _load_cli_files(view_keys, vib_mode): """Parse CLI molecule arguments into a viewer payload. Returns: - tuple: ``(to_pop, mols, vib, last_ret, orca_handled)``. + tuple: ``(to_pop, mols, vib, last_ret, plugin_handled)``. """ to_pop = [] mols = [] vib = None last_ret = 0 - orca_handled = False + plugin_handled = False for i, arg in enumerate(sys.argv[1:], start=1): if ":" in arg: @@ -251,14 +241,15 @@ def _load_cli_files(view_keys, vib_mode): warnings.warn(str(e), RuntimeWarning, stacklevel=2) continue - if plugin_for(peek) is not None: + plugin = plugin_for(peek) + if plugin is not None: remaining = [sys.argv[j] for j in range(len(sys.argv)) if j not in popped] try: - last_ret = _handle_orca(arg, view_keys, vib_mode, remaining) + last_ret = _handle_plugin(plugin, arg, view_keys, vib_mode, remaining) except RuntimeError as e: warnings.warn(str(e), RuntimeWarning, stacklevel=2) continue - orca_handled = True + plugin_handled = True continue try: @@ -267,22 +258,20 @@ def _load_cli_files(view_keys, vib_mode): warnings.warn(str(e), RuntimeWarning, stacklevel=2) continue - if not ( - hasattr(data, "atomcoords") - and hasattr(data, "atomnos") - and len(data.atomcoords) - ): - warnings.warn(f"No coordinates found in {arg}.", RuntimeWarning, stacklevel=2) + if not has_coords(data): + warnings.warn( + f"No coordinates found in {arg}.", RuntimeWarning, stacklevel=2, + ) continue m = _mols_from_data(data, parser) - v = _vib_from_data(data) + v = vib_from_data(data) if not mols and v is not None and vib_mode is not False: mols, vib = m, v break mols.extend(m) - return to_pop, mols, vib, last_ret, orca_handled + return to_pop, mols, vib, last_ret, plugin_handled def main(): @@ -291,9 +280,9 @@ def main(): Usage: vmol2.py file [options] - ORCA outputs may offer frequencies, optimization, single point, scan, - and/or IRC. With several datasets and a TTY, a menu is shown. Otherwise - use ``view:freq|opt|sp|irc|scan`` or ``vib:0`` / ``vib:1``. + ORCA and Gaussian outputs may offer frequencies, optimization, single + point, scan, and/or IRC. With several datasets and a TTY, a menu is shown. + Otherwise use ``view:freq|opt|sp|irc|scan`` or ``vib:0`` / ``vib:1``. Returns: int: Exit code (0 for success). @@ -318,9 +307,9 @@ def main(): for i in reversed([j + 1 for j in view_pops]): sys.argv.pop(i) - to_pop, mols, vib, last_ret, orca_handled = _load_cli_files(view_keys, vib_mode) + to_pop, mols, vib, last_ret, plugin_handled = _load_cli_files(view_keys, vib_mode) - if orca_handled and not mols: + if plugin_handled and not mols: return last_ret if not mols: From 3ec54a7c0470a88ce36d8ed574ac17c11e700c3d Mon Sep 17 00:00:00 2001 From: rlaplaza Date: Thu, 15 Oct 2026 12:34:48 +0200 Subject: [PATCH 3/5] Turn frequency view into a selectable mode-displacement trajectory Listing hundreds of native vib modes is hard to use; prompt for mode, scale, and frame count (or mode:/scale:/npoints:) and open geometry frames. Co-authored-by: Cursor --- .github/CHANGELOG.md | 2 +- python/README.md | 9 +- python/ruff.toml | 1 + python/tests/test_orca_plugin.py | 12 +- python/tests/test_vibtraj.py | 93 +++++++++++++++ python/vmol/plugins/vibtraj.py | 91 ++++++++++++++ python/vmol/vmol2.py | 199 ++++++++++++++++++++++++------- 7 files changed, 362 insertions(+), 45 deletions(-) create mode 100644 python/tests/test_vibtraj.py create mode 100644 python/vmol/plugins/vibtraj.py diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 0dfc9f2..8b4c160 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -37,7 +37,7 @@ ### Coming in the next version: * high-symmetry determination bugs (#21) * headless fallback build by @JamesOBrien2 (#44) -* ORCA / Gaussian plugins for `vmol2`: prompt / `view:` selectors for frequencies, optimization, single point, scan, and IRC trajectories +* ORCA / Gaussian plugins for `vmol2`: prompt / `view:` selectors for frequencies, optimization, single point, scan, and IRC trajectories; frequency view builds a selectable mode-displacement trajectory (`mode:` / `scale:` / `npoints:`) **Full Changelog**: https://github.com/briling/v/compare/v2.0...v3.0rc7 diff --git a/python/README.md b/python/README.md index adfc0b3..a9bf97e 100644 --- a/python/README.md +++ b/python/README.md @@ -56,13 +56,18 @@ ORCA and Gaussian logs use dedicated plugins. When a file has several of frequencies, optimization, trailing single point (ORCA), scan, or IRC, `vmol2` prompts on a TTY (each choice opens in its own `v` session), or accepts: ```bash -vmol2 job.out view:freq # normal modes (on the opt minimum) +vmol2 job.out view:freq # pick a normal mode, then show a displacement path vmol2 job.out view:opt # geometry optimization frames vmol2 job.out view:sp # final single-point geometry / energy (ORCA multi-job) vmol2 job.out view:irc # IRC trajectory vmol2 job.out view:scan # relaxed surface scan vmol2 job.out vib:0 # geometry path (opt, else SP, else IRC, else scan) -vmol2 job.out vib:1 # frequencies +vmol2 job.out vib:1 # frequencies (same as view:freq) +``` +Choosing frequencies (TTY) lists modes with their cm⁻¹ values, then asks for a +mode number, displacement magnitude (Å), and number of frames. Non-interactive: +```bash +vmol2 job.out view:freq mode:1 scale:0.5 npoints:21 ``` ORCA IRC frames come from `basename_IRC_Full_trj.xyz` (or the path in the log's `Storing full IRC trajectory in ....` line). Gaussian IRC frames are read from diff --git a/python/ruff.toml b/python/ruff.toml index 3df1a61..2dd3d31 100644 --- a/python/ruff.toml +++ b/python/ruff.toml @@ -54,3 +54,4 @@ ignore = [ "vmol/__main__.py" = ["D"] "examples/ex*.py" = ["D"] "tests/*.py" = ["SLF", "PLC2701", "D"] +"vmol/plugins/vibtraj.py" = ["PLC2701"] diff --git a/python/tests/test_orca_plugin.py b/python/tests/test_orca_plugin.py index 7702745..a9922f7 100644 --- a/python/tests/test_orca_plugin.py +++ b/python/tests/test_orca_plugin.py @@ -231,13 +231,23 @@ def test_view_and_vib_selectors(self): ["opt"], ) - keys, vib, pops = vmol2._collect_view_selectors( + keys, vib, mode, scale, npoints, pops = vmol2._collect_cli_selectors( ["file.out", "view:sp", "vib:0", "bonds:0"], ) self.assertEqual(keys, ["sp"]) self.assertIs(vib, False) + self.assertIsNone(mode) self.assertEqual(pops, [1]) + keys, vib, mode, scale, npoints, pops = vmol2._collect_cli_selectors( + ["file.out", "view:freq", "mode:2", "scale:0.3", "npoints:9"], + ) + self.assertEqual(keys, ["freq"]) + self.assertEqual(mode, 2) + self.assertEqual(scale, 0.3) + self.assertEqual(npoints, 9) + self.assertEqual(pops, [1, 2, 3, 4]) + if __name__ == "__main__": unittest.main() diff --git a/python/tests/test_vibtraj.py b/python/tests/test_vibtraj.py new file mode 100644 index 0000000..8cff504 --- /dev/null +++ b/python/tests/test_vibtraj.py @@ -0,0 +1,93 @@ +"""Unit tests for normal-mode displacement trajectories.""" + +from __future__ import annotations + +import unittest + +import numpy as np + +from vmol import vmol2 +from vmol.plugins.common import Dataset +from vmol.plugins.vibtraj import mode_labels, mode_trajectory + + +class TestModeTrajectory(unittest.TestCase): + """Geometry frames along a normalized mode.""" + + def test_h2_stretch_endpoints(self): + """Endpoints are ±A along the mode; midpoint is the equilibrium.""" + mol0 = { + "q": [1, 1], + "r": np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.74]]), + "name": "H2", + } + vib = { + "freq": np.array([4401.0]), + "disp": np.array([[[0.0, 0.0, -2.0], [0.0, 0.0, 2.0]]]), + "ints": np.array([0.0]), + } + scale = 0.5 + mols = mode_trajectory(mol0, vib, 0, scale=scale, npoints=5) + self.assertEqual(len(mols), 5) + + u = vib["disp"][0] / np.linalg.norm(vib["disp"][0]) + np.testing.assert_allclose(mols[0]["r"], mol0["r"] - scale * u) + np.testing.assert_allclose(mols[2]["r"], mol0["r"], atol=1e-12) + np.testing.assert_allclose(mols[4]["r"], mol0["r"] + scale * u) + self.assertIn("4401.00 cm-1", mols[0]["name"]) + self.assertIn("mode 1/1", mols[0]["name"]) + + def test_mode_labels_include_imag_and_ir(self): + """Labels show frequency, imag tag, and IR intensity when present.""" + vib = { + "freq": np.array([-66.0, 100.5]), + "disp": np.zeros((2, 1, 3)), + "ints": np.array([12.3, 0.0]), + } + labels = mode_labels(vib) + self.assertEqual(len(labels), 2) + self.assertIn("-66.00", labels[0]) + self.assertIn("imag", labels[0]) + self.assertIn("IR=12.3", labels[0]) + self.assertIn("100.50", labels[1]) + + +class TestFreqSelectors(unittest.TestCase): + """CLI selectors and freq expansion.""" + + def test_collect_mode_scale_npoints(self): + """mode:/scale:/npoints: are parsed and marked for removal.""" + keys, vib, mode, scale, npoints, pops = vmol2._collect_cli_selectors( + ["file.out", "view:freq", "mode:3", "scale:0.25", "npoints:11", "bonds:0"], + ) + self.assertEqual(keys, ["freq"]) + self.assertIsNone(vib) + self.assertEqual(mode, 3) + self.assertEqual(scale, 0.25) + self.assertEqual(npoints, 11) + self.assertEqual(pops, [1, 2, 3, 4]) + + def test_freq_to_trajectory_noninteractive(self): + """Non-TTY expansion uses CLI mode/scale/npoints without prompting.""" + mol0 = { + "q": [1, 1], + "r": np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]), + "name": "H2", + } + vib = { + "freq": np.array([100.0, 200.0]), + "disp": np.array([ + [[0.0, 0.0, -1.0], [0.0, 0.0, 1.0]], + [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0]], + ]), + } + ds = Dataset(key="freq", label="f", mols=[mol0], vib=vib) + out = vmol2._freq_to_trajectory(ds, (2, 0.4, 3)) + self.assertIsNone(out.vib) + self.assertEqual(len(out.mols), 3) + self.assertIn("mode 2/2", out.mols[0]["name"]) + self.assertIn("200.00 cm-1", out.mols[0]["name"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/python/vmol/plugins/vibtraj.py b/python/vmol/plugins/vibtraj.py new file mode 100644 index 0000000..30ab530 --- /dev/null +++ b/python/vmol/plugins/vibtraj.py @@ -0,0 +1,91 @@ +"""Build geometry trajectories along a selected normal mode.""" + +from __future__ import annotations + +import numpy as np + +from vmol.plugins.common import _clip_name + +_MIN_NPOINTS = 2 +DEFAULT_SCALE = 0.5 +DEFAULT_NPOINTS = 21 + + +def mode_labels(vib: dict) -> list[str]: + """Return one label per mode for interactive listing. + + Returns: + list[str]: Labels like ``1) -66.03 cm-1 IR=12.3``. + """ + freqs = np.asarray(vib["freq"], dtype=float) + ints = vib.get("ints") + if ints is not None: + ints = np.asarray(ints, dtype=float) + labels = [] + width = max(1, len(str(len(freqs)))) + for i, fq in enumerate(freqs, start=1): + line = f"{i:>{width}}) {fq:8.2f} cm-1" + if fq < 0: + line += " imag" + if ints is not None and i - 1 < len(ints): + line += f" IR={ints[i - 1]:.1f}" + labels.append(line) + return labels + + +def mode_trajectory( + mol0: dict, + vib: dict, + mode_index: int, + *, + scale: float = DEFAULT_SCALE, + npoints: int = DEFAULT_NPOINTS, +) -> list[dict]: + """Displace *mol0* along vib mode *mode_index* (0-based). + + Frames use ``r = r0 + A * t * u`` with ``t`` in ``linspace(-1, 1, npoints)`` + and ``u`` the L2-normalized cartesian displacement. + + Returns: + list[dict]: Molecule frames suitable for ``vmol.run(mols=...)``. + + Raises: + ValueError: If mode index / scale / npoints are invalid. + """ + freqs = np.asarray(vib["freq"], dtype=float) + disps = np.asarray(vib["disp"], dtype=float) + nmodes = len(freqs) + if not 0 <= mode_index < nmodes: + msg = f"mode index {mode_index} out of range (0..{nmodes - 1})" + raise ValueError(msg) + if npoints < _MIN_NPOINTS: + msg = f"npoints must be >= {_MIN_NPOINTS}, got {npoints}" + raise ValueError(msg) + if scale < 0: + msg = f"scale must be >= 0, got {scale}" + raise ValueError(msg) + + r0 = np.asarray(mol0["r"], dtype=float) + u = np.asarray(disps[mode_index], dtype=float) + if u.shape != r0.shape: + msg = f"displacement shape {u.shape} does not match coords {r0.shape}" + raise ValueError(msg) + norm = float(np.linalg.norm(u.reshape(-1))) + if norm == 0.0: + msg = f"mode {mode_index + 1} has zero displacement" + raise ValueError(msg) + u /= norm + + fq = float(freqs[mode_index]) + ts = np.linspace(-1.0, 1.0, npoints) + base = mol0.get("name", "mode") + mols = [] + for i, t in enumerate(ts): + r = r0 + scale * float(t) * u + name = ( + f"{base} mode {mode_index + 1}/{nmodes} " + f"{fq:.2f} cm-1 t={t:+.3f} A={scale:.3f} A " + f"frame {i + 1}/{npoints}" + ) + mols.append({"q": mol0["q"], "r": r, "name": _clip_name(name)}) + return mols diff --git a/python/vmol/vmol2.py b/python/vmol/vmol2.py index c3bb8fe..377b6c5 100755 --- a/python/vmol/vmol2.py +++ b/python/vmol/vmol2.py @@ -7,7 +7,13 @@ from vmol import vmol, _paths from vmol.plugins import plugin_for -from vmol.plugins.common import has_coords, parse_cclib, vib_from_data +from vmol.plugins.common import Dataset, has_coords, parse_cclib, vib_from_data +from vmol.plugins.vibtraj import ( + DEFAULT_NPOINTS, + DEFAULT_SCALE, + mode_labels, + mode_trajectory, +) try: import cclib @@ -81,17 +87,20 @@ def read_mols_from_cclib(path): return _mols_from_data(data, parser), vib_from_data(data) -def _collect_view_selectors(argv): - """Parse ``view:`` / ``vib:`` tokens from *argv*. +def _collect_cli_selectors(argv): + """Parse ``view:`` / ``vib:`` / ``mode:`` / ``scale:`` / ``npoints:`` from *argv*. Returns: - tuple: ``(view_keys, vib_mode, indices_to_pop)``. + tuple: ``(view_keys, vib_mode, mode, scale, npoints, indices_to_pop)``. Raises: - ValueError: If a ``view:`` key is unknown. + ValueError: If a ``view:`` key is unknown or a selector value is invalid. """ view_keys = [] vib_mode = None + mode = None + scale = None + npoints = None to_pop = [] for i, arg in enumerate(argv): if arg.startswith("view:"): @@ -108,7 +117,16 @@ def _collect_view_selectors(argv): vib_mode = bool(int(arg.split(":", 1)[1])) except ValueError: continue - return view_keys, vib_mode, to_pop + elif arg.startswith("mode:"): + mode = int(arg.split(":", 1)[1]) + to_pop.append(i) + elif arg.startswith("scale:"): + scale = float(arg.split(":", 1)[1]) + to_pop.append(i) + elif arg.startswith("npoints:"): + npoints = int(arg.split(":", 1)[1]) + to_pop.append(i) + return view_keys, vib_mode, mode, scale, npoints, to_pop def _default_datasets(datasets): @@ -122,13 +140,9 @@ def _prompt_datasets(datasets): print(f"Found {len(datasets)} dataset(s) in this file:") for i, ds in enumerate(datasets, start=1): print(f" {i}) [{ds.key}] {ds.label}") - print("Enter number(s), comma-separated list, or 'all': ", end="", flush=True) - try: - line = sys.stdin.readline() - except EOFError: - line = "" - - choice = line.strip().lower() + choice = _read_prompt_line( + "Enter number(s), comma-separated list, or 'all': ", + ).strip().lower() if not choice: return _default_datasets(datasets) if choice in {"all", "*"}: @@ -200,29 +214,138 @@ def _choose_datasets(datasets, view_keys, vib_mode): return _default_datasets(datasets) -def _run_datasets(datasets, args): +def _read_prompt_line(prompt: str) -> str: + print(prompt, end="", flush=True) + try: + return sys.stdin.readline() + except EOFError: + return "" + + +def _prompt_or_default(prompt: str, default, cast): + if not sys.stdin.isatty(): + return default + raw = _read_prompt_line(prompt).strip() + return cast(raw) if raw else default + + +def _resolve_freq_params(vib, mode, scale, npoints): + """Resolve mode (1-based), scale, and npoints via CLI values and/or TTY prompts. + + Returns: + tuple: ``(mode_0based, scale, npoints)``. + + Raises: + ValueError: If the mode index is out of range. + """ + nmodes = len(vib["freq"]) + if mode is None: + if sys.stdin.isatty(): + print(f"Found {nmodes} vibrational mode(s):") + for line in mode_labels(vib): + print(f" {line}") + mode = _prompt_or_default( + f"Enter mode number (1-{nmodes}): ", 1, int, + ) + else: + warnings.warn( + "No mode:N given for frequencies; using mode 1.", + RuntimeWarning, + stacklevel=2, + ) + mode = 1 + + if not 1 <= mode <= nmodes: + msg = f"mode:{mode} out of range (1-{nmodes})" + raise ValueError(msg) + + if scale is None: + scale = _prompt_or_default( + f"Displacement magnitude in Angstrom [{DEFAULT_SCALE}]: ", + DEFAULT_SCALE, + float, + ) + if npoints is None: + npoints = _prompt_or_default( + f"Number of frames [{DEFAULT_NPOINTS}]: ", + DEFAULT_NPOINTS, + int, + ) + return mode - 1, float(scale), int(npoints) + + +def _freq_to_trajectory(ds, freq_opts) -> Dataset: + if ds.vib is None or not ds.mols: + msg = "Frequency dataset is missing vibration or geometry data." + raise RuntimeError(msg) + mode0, scale, npoints = _resolve_freq_params(ds.vib, *freq_opts) + mols = mode_trajectory( + ds.mols[0], ds.vib, mode0, scale=scale, npoints=npoints, + ) + fq = float(ds.vib["freq"][mode0]) + return Dataset( + key="freq", + label=f"Mode {mode0 + 1} ({fq:.2f} cm-1, A={scale:g} A, {npoints} frames)", + mols=mols, + vib=None, + ) + + +def _run_datasets(datasets, args, freq_opts): ret = 0 for ds in datasets: + if ds.key == "freq": + ds = _freq_to_trajectory(ds, freq_opts) ret = vmol.run(mols=ds.mols, vib=ds.vib, args=args, with_arg0=True) return ret -def _handle_plugin(plugin, path, view_keys, vib_mode, args): +def _handle_plugin(plugin, path, select, args, freq_opts): + view_keys, vib_mode = select _parser, data = plugin.parse(path) datasets = plugin.inventory(path, data) chosen = _choose_datasets(datasets, view_keys, vib_mode) if not chosen: msg = f"No datasets selected for {path}." raise RuntimeError(msg) - return _run_datasets(chosen, args) + return _run_datasets(chosen, args, freq_opts) -def _load_cli_files(view_keys, vib_mode): +def _try_generic_file(arg, vib_mode, freq_opts, mols, vib): + """Load a non-plugin cclib file into the accumulating mol list. + + Returns: + tuple: ``(mols, vib, done)`` where *done* means stop loading files. + """ + try: + parser, data = _parse_file(arg) + except (FileNotFoundError, RuntimeError) as e: + warnings.warn(str(e), RuntimeWarning, stacklevel=3) + return mols, vib, False + + if not has_coords(data): + warnings.warn( + f"No coordinates found in {arg}.", RuntimeWarning, stacklevel=3, + ) + return mols, vib, False + + m = _mols_from_data(data, parser) + v = vib_from_data(data) + if not mols and v is not None and vib_mode is not False: + fake = Dataset(key="freq", label="frequencies", mols=m[:1], vib=v) + traj = _freq_to_trajectory(fake, freq_opts) + return traj.mols, None, True + mols.extend(m) + return mols, vib, False + + +def _load_cli_files(select, freq_opts): """Parse CLI molecule arguments into a viewer payload. Returns: tuple: ``(to_pop, mols, vib, last_ret, plugin_handled)``. """ + vib_mode = select[1] to_pop = [] mols = [] vib = None @@ -233,7 +356,6 @@ def _load_cli_files(view_keys, vib_mode): if ":" in arg: continue to_pop.append(i) - popped = set(to_pop) try: peek = _peek_parser(arg) @@ -243,33 +365,22 @@ def _load_cli_files(view_keys, vib_mode): plugin = plugin_for(peek) if plugin is not None: - remaining = [sys.argv[j] for j in range(len(sys.argv)) if j not in popped] + remaining = [ + sys.argv[j] for j in range(len(sys.argv)) if j not in set(to_pop) + ] try: - last_ret = _handle_plugin(plugin, arg, view_keys, vib_mode, remaining) - except RuntimeError as e: + last_ret = _handle_plugin( + plugin, arg, select, remaining, freq_opts, + ) + except (RuntimeError, ValueError) as e: warnings.warn(str(e), RuntimeWarning, stacklevel=2) continue plugin_handled = True continue - try: - parser, data = _parse_file(arg) - except (FileNotFoundError, RuntimeError) as e: - warnings.warn(str(e), RuntimeWarning, stacklevel=2) - continue - - if not has_coords(data): - warnings.warn( - f"No coordinates found in {arg}.", RuntimeWarning, stacklevel=2, - ) - continue - - m = _mols_from_data(data, parser) - v = vib_from_data(data) - if not mols and v is not None and vib_mode is not False: - mols, vib = m, v + mols, vib, done = _try_generic_file(arg, vib_mode, freq_opts, mols, vib) + if done: break - mols.extend(m) return to_pop, mols, vib, last_ret, plugin_handled @@ -282,7 +393,9 @@ def main(): ORCA and Gaussian outputs may offer frequencies, optimization, single point, scan, and/or IRC. With several datasets and a TTY, a menu is shown. - Otherwise use ``view:freq|opt|sp|irc|scan`` or ``vib:0`` / ``vib:1``. + For frequencies, pick a mode and displacement (or pass ``mode:`` / + ``scale:`` / ``npoints:``). Otherwise use ``view:freq|opt|sp|irc|scan`` + or ``vib:0`` / ``vib:1``. Returns: int: Exit code (0 for success). @@ -303,11 +416,15 @@ def main(): vmol.run(args=sys.argv, with_arg0=True) return 0 - view_keys, vib_mode, view_pops = _collect_view_selectors(sys.argv[1:]) + view_keys, vib_mode, mode, scale, npoints, view_pops = _collect_cli_selectors( + sys.argv[1:], + ) for i in reversed([j + 1 for j in view_pops]): sys.argv.pop(i) - to_pop, mols, vib, last_ret, plugin_handled = _load_cli_files(view_keys, vib_mode) + select = (view_keys, vib_mode) + freq_opts = (mode, scale, npoints) + to_pop, mols, vib, last_ret, plugin_handled = _load_cli_files(select, freq_opts) if plugin_handled and not mols: return last_ret From a541688d6eb18df72027221999a9553e41122fe6 Mon Sep 17 00:00:00 2001 From: rlaplaza Date: Thu, 15 Oct 2026 12:36:41 +0200 Subject: [PATCH 4/5] Fix ruff 0.16 noqa style and float-equality lint in plugins Match the CI `ruff check vmol/ examples/` command so the PR lint job passes. Co-authored-by: Cursor --- python/vmol/plugins/common.py | 2 +- python/vmol/plugins/orca_parse.py | 4 ++-- python/vmol/plugins/vibtraj.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/vmol/plugins/common.py b/python/vmol/plugins/common.py index b77d43a..46dcee3 100644 --- a/python/vmol/plugins/common.py +++ b/python/vmol/plugins/common.py @@ -61,7 +61,7 @@ def parse_cclib(path: str): raise RuntimeError(msg) try: data = parser.parse() - except Exception as exc: # noqa: BLE001 + except Exception as exc: # ruff: ignore[blind-except] warnings.warn( f"Error parsing {path} with cclib ({type(exc).__name__}: {exc}). " "Using partial data.", diff --git a/python/vmol/plugins/orca_parse.py b/python/vmol/plugins/orca_parse.py index 6ec9dbc..2e18ce2 100644 --- a/python/vmol/plugins/orca_parse.py +++ b/python/vmol/plugins/orca_parse.py @@ -36,7 +36,7 @@ def patched_append(self, inputfile, line): self.scfvalues.append([]) try: return orig_append(self, inputfile, line) - except Exception as exc: # noqa: BLE001 + except Exception as exc: # ruff: ignore[blind-except] warnings.warn( f"ORCA SCF convergence summary skipped ({type(exc).__name__}: {exc})", RuntimeWarning, @@ -49,7 +49,7 @@ def patched_append(self, inputfile, line): def safe_extract(self, inputfile, line): try: return orig_extract(self, inputfile, line) - except Exception as exc: # noqa: BLE001 + except Exception as exc: # ruff: ignore[blind-except] warnings.warn( f"ORCA section skipped ({type(exc).__name__}: {exc})", RuntimeWarning, diff --git a/python/vmol/plugins/vibtraj.py b/python/vmol/plugins/vibtraj.py index 30ab530..bb8c570 100644 --- a/python/vmol/plugins/vibtraj.py +++ b/python/vmol/plugins/vibtraj.py @@ -71,7 +71,7 @@ def mode_trajectory( msg = f"displacement shape {u.shape} does not match coords {r0.shape}" raise ValueError(msg) norm = float(np.linalg.norm(u.reshape(-1))) - if norm == 0.0: + if not norm: msg = f"mode {mode_index + 1} has zero displacement" raise ValueError(msg) u /= norm From 4db2744e7a7d6c84643d57fbabddafa9a6178607 Mon Sep 17 00:00:00 2001 From: rlaplaza Date: Thu, 15 Oct 2026 12:46:20 +0200 Subject: [PATCH 5/5] Silence expected ORCA 6 / cclib section-skip warnings Known parse gaps are skipped quietly so opening multi-job ORCA logs does not spam RuntimeWarnings on every incomplete section. Co-authored-by: Cursor --- python/vmol/plugins/orca_parse.py | 39 +++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/python/vmol/plugins/orca_parse.py b/python/vmol/plugins/orca_parse.py index 2e18ce2..d00c4b2 100644 --- a/python/vmol/plugins/orca_parse.py +++ b/python/vmol/plugins/orca_parse.py @@ -15,6 +15,15 @@ _PT = PeriodicTable() _PATCHED = False _CART_MIN_TOKENS = 4 +# Expected on ORCA 6 / incomplete blocks with PyPI cclib; skip quietly. +_BENIGN_PARSE_ERRORS = ( + ValueError, + AssertionError, + StopIteration, + IndexError, + AttributeError, + KeyError, +) def _patch_cclib_orca() -> None: @@ -22,6 +31,7 @@ def _patch_cclib_orca() -> None: PyPI cclib 1.8.1 aborts on empty ``scfvalues`` and on some ORCA 6 sections. Skip those sections so geometries and frequencies still get collected. + Expected ORCA 6 parse gaps are silent; unexpected errors still warn once. """ global _PATCHED if _PATCHED: @@ -36,12 +46,17 @@ def patched_append(self, inputfile, line): self.scfvalues.append([]) try: return orig_append(self, inputfile, line) + except _BENIGN_PARSE_ERRORS: + return None except Exception as exc: # ruff: ignore[blind-except] - warnings.warn( - f"ORCA SCF convergence summary skipped ({type(exc).__name__}: {exc})", - RuntimeWarning, - stacklevel=2, - ) + if not getattr(self, "_vmol_warned_scf", False): + warnings.warn( + f"ORCA SCF convergence summary skipped " + f"({type(exc).__name__}: {exc})", + RuntimeWarning, + stacklevel=2, + ) + self._vmol_warned_scf = True return None orig_extract = orcaparser.ORCA.extract @@ -49,12 +64,16 @@ def patched_append(self, inputfile, line): def safe_extract(self, inputfile, line): try: return orig_extract(self, inputfile, line) + except _BENIGN_PARSE_ERRORS: + return None except Exception as exc: # ruff: ignore[blind-except] - warnings.warn( - f"ORCA section skipped ({type(exc).__name__}: {exc})", - RuntimeWarning, - stacklevel=2, - ) + if not getattr(self, "_vmol_warned_extract", False): + warnings.warn( + f"ORCA section skipped ({type(exc).__name__}: {exc})", + RuntimeWarning, + stacklevel=2, + ) + self._vmol_warned_extract = True return None orcaparser.ORCA._append_scfvalues_scftargets = patched_append