diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b9a811a..ed62a39 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11"] + python-version: ["3.9", "3.10", "3.11"] steps: - uses: actions/checkout@v4 diff --git a/AUDIT_FINDINGS.md b/AUDIT_FINDINGS.md new file mode 100644 index 0000000..ed228b1 --- /dev/null +++ b/AUDIT_FINDINGS.md @@ -0,0 +1,185 @@ +# QuantUI Repository Audit — Findings + +**Date:** 2026-07-14 +**Scope:** Full read of the core calculation, persistence, resolver, and infrastructure modules (`quantui/*.py`, ~29k LOC), targeted review of the large UI modules (`app*.py`, `benchmarks.py`, visualization modules), plus config/CI/packaging. Static analysis (ruff, bug-focused rule set) came back clean apart from `raise … from` hygiene, so everything below was found by manual review. + +Findings are ordered by severity. File/line references are to the audited commit (`310a8cf`). + +**Status (2026-07-14):** all 5 High, all 14 Medium, and 13 of 14 Low findings below are fixed, regression-tested, and merged on `claude/repo-audit-findings-yolusq`. Only item #12 (`prune_events` malformed-entry retention) is left as-is, by explicit choice — its original "deliberate" behavior stays unchanged. + +--- + +## High — functional bugs + +### H1. Orbital isosurface fails for charged and all open-shell molecules — ✅ Fixed +`orbital_visualization.generate_cube_from_arrays()` builds the PySCF molecule as: + +```python +mol = gto.M(atom=mol_atom, basis=mol_basis, unit="Angstrom", spin=spin) # orbital_visualization.py:673 +``` + +No `charge` is passed, and the only call site passes neither `charge` nor `spin`: + +```python +generate_cube_from_arrays(mol_atom, mol_basis, mo_coeff, orb_idx, cube_path) # app_visualization.py:1207 +``` + +PySCF raises at `mol.build()` whenever the electron count implied by `charge=0` is inconsistent with `spin=0` — i.e. for **every odd-electron system**: radicals (multiplicity 2), and charged species like H₃O⁺, NH₄⁺, OH⁻. The isosurface panel therefore errors out for these systems even though the MO coefficients are available and valid. `generate_cube_file()` (line 597) has the same omission. + +**Fix:** thread `charge` and `spin` (`multiplicity − 1`) from the result object through to both cube generators. + +### H2. Frequency results store `pyscf_mol_atom` in **Bohr**; every consumer assumes **Ångström** — ✅ Fixed +`freq_calc.py:286` populates the field from PySCF internals: + +```python +pyscf_mol_atom = [(str(s), list(map(float, c))) for s, c in mol._atom] # mol._atom is in Bohr +``` + +whereas `session_calc.py:529` and `optimizer.py` populate the same field in Ångström (as documented at `session_calc.py:90`). Consumers all assume Å: + +- `results_storage.save_molden()` (line 303) builds `mol.atom = …` with PySCF's default Å unit → the `[Atoms]`/`[GTO]` sections of Molden files exported from **frequency runs have the geometry inflated by ×1.889** (Bohr values read as Å), while the `[FR-COORD]` block (which the Molden spec defines in Bohr, written verbatim at line 368) is coincidentally correct — so the exported file is internally inconsistent and renders wrong in Avogadro/IQmol. +- `results_storage.save_orbitals()` persists the Bohr coordinates to `orbitals_meta.json`, so **history replay** of a frequency result feeds Bohr coordinates into `generate_cube_from_arrays(..., unit="Angstrom")`. +- `app_visualization.py:968` caches the same field for live isosurface rendering. + +**Fix:** in `freq_calc`, build the field from `molecule.atoms`/`molecule.coordinates` (like `session_calc` does), or convert `mol._atom` with `param.BOHR`; and have `_append_molden_vibrations` do an explicit Å→Bohr conversion for `[FR-COORD]` so its input convention matches everything else. + +### H3. "Export Script" and the method-notes panel are broken for `wB97X-D` — ✅ Fixed +`calculator.PySCFCalculation.__init__` upper-cases the method then checks membership against the mixed-case list: + +```python +self.method = method.upper() # "wB97X-D" -> "WB97X-D" +if self.method not in config.SUPPORTED_METHODS: # list contains "wB97X-D" + raise ValueError(...) # calculator.py:47-53 +``` + +The UI passes display names straight from the dropdown (options = `SUPPORTED_METHODS`, `app_builders.py:579`), so: + +- **Export Script** (`app_exports.on_export`, line 17) fails with "Method 'wB97X-D' not supported". +- **Educational notes** (`app_runflow.update_notes`, line 1068) silently disappear for wB97X-D — the exception is swallowed at line 1089. + +Even if construction succeeded, the generated script's `_XC_ALIAS` / `_NEEDS_D3` tables (config.py:444-450) key on `'wB97X-D'` and would not match the upper-cased `'WB97X-D'`, producing a script that sets an invalid `mf.xc`. Same latent issue applies to any future mixed-case method name. + +**Fix:** compare case-insensitively while preserving the canonical display name (mirror `session_calc.resolve_xc`'s approach). + +### H4. `generate_2d_structure_svg()`'s XYZ input path can never work — ✅ Fixed +`pubchem.py:808-825` builds the molecule with the immutable class: + +```python +rdkit_mol = Chem.Mol() +... +rdkit_mol.AddAtom(atom) # AddAtom exists only on Chem.RWMol → AttributeError +``` + +The AttributeError is swallowed by the broad `except` at line 841, so `xyz_string=` input always returns `None`. Two further latent bugs in the same block: the conformer index uses the raw line index `i` (lines with `len(parts) < 4` are skipped but `i` still advances, desynchronizing atoms and coordinates), and `Chem.Conformer()` is created without the atom count. + +**Fix:** use `Chem.RWMol()`, track a separate atom counter, size the conformer, and call `.GetMol()` at the end. + +### H5. XYZ parser silently drops the first atom for standard files with a blank or `#` title line — ✅ Fixed +`molecule.parse_xyz_input()` strips blank lines and full-line comments **before** detecting the `count / title / atoms…` header layout (lines 364-405). For a completely standard file such as: + +``` +3 + ← blank title line (very common) +O 0.0 0.0 0.0 +H ... +H ... +``` + +the blank line is removed, then `start_idx = 2` skips the count line **and the `O` atom line**. The count mismatch is only a `logger.warning` (line 528), so the user gets a 2-atom water with no visible error. The same happens when the title line starts with `#` or `!`. A pasted `N\natom…` block without any title line also loses its first atom. + +**Fix:** detect the header on the raw line list (or only strip comments after header detection), and consider promoting the count-mismatch warning to a hard error. + +--- + +## Medium — incorrect behavior on specific paths / edge cases + +### M1. Element support is H–Kr, but error text and the resolver chain suggest otherwise — ✅ Fixed +`config.VALID_ATOMS` stops at Kr (Z=36), yet: + +- The parse error help text explicitly lists **`I` (iodine) as a valid symbol** (`molecule.py:485`). +- PubChem / CACTUS / SMILES resolution happily returns structures containing I, Sn, Ag, … which then fail `Molecule` validation with a misleading message (e.g. searching *thyroxine* or *amiodarone* resolves, then can't be loaded). +- `Molecule.get_electron_count()` returns 0 for unknown elements (`atomic_numbers.get(atom, 0)`), which would silently corrupt electron counts if the validation gate were ever relaxed. + +### M2. MP2 / CCSD / CCSD(T) silently use an RHF/ROHF reference for open-shell molecules — ⚠️ Re-scoped, fixed a different real bug +**Status:** empirical testing showed this finding's premise was wrong — PySCF's `scf.RHF()` factory already dispatches to ROHF correctly for multiplicity > 1, and MP2/CCSD/CCSD(T) run correctly against it. Investigating this surfaced the actual bug: `optimize_geometry`, `run_freq_calc`, `run_pes_scan`, `run_tddft_calc`, and `run_nmr_calc` had no guard against post-HF methods at all and crashed uninformatively for any calc type other than Single Point. Added an explicit `POST_HF_METHODS` guard with a clear `ValueError` to all five entry points instead. +`session_calc.py:331-337` builds `scf.RHF(mol)` for these methods regardless of spin. With multiplicity > 1, PySCF's `scf.RHF` factory returns ROHF; the post-HF kernels may fail or produce results inconsistent with the closed-shell "RHF reference" the docstrings and result fields describe. No guard or warning exists. + +### M3. UHF results silently lack dipole moment and Mulliken charges — ✅ Fixed +`session_calc.py:490` skips both for `method_upper == "UHF"` only — while UKS (open-shell DFT) flows through the same extraction successfully. Both properties are well-defined for UHF; the skip appears stale rather than principled, and the result card just shows nothing. + +### M4. NMR silently falls back to B3LYP/6-31G* TMS reference constants — ✅ Fixed +`nmr_calc.py:319-322`: any method/basis combination absent from `NMR_REFERENCE_SHIELDINGS` (e.g. M06-2X, wB97X-D, def2-TZVP, cc-pVTZ) silently uses the B3LYP/6-31G* constants, and `NMRResult` does not record which reference was applied. Chemical shifts can be systematically offset by several ppm (¹³C) with no indication to the student. The lookup is also case-sensitive. + +### M5. IR-intensity worker rebuilds RHF/UHF from spin alone; failure mode contradicts the docs — ✅ Fixed +`freq_ir_workers.run_displaced_scf` (line 149) and the serial `_displaced_scf_dipole` (freq_calc.py:390) pick `RHF if spin == 0 else UHF`. A UHF-singlet parent run becomes RHF in the displaced SCFs, and the parent's `(2, nao, nao)` UHF `dm0` initial guess then has the wrong shape. The resulting exception aborts the whole IR-intensity step — and contrary to `freq_ir_workers`' module docstring ("the driver … falls back to the serial loop so the user's calc still completes"), `freq_calc.py:555` does **not** fall back to serial; it skips IR intensities entirely. + +### M6. PES scan records the wrong geometry and NaN-poisoned statistics on failed points — ✅ Fixed +`pes_scan.py:329-333`: a failed scan point appends the **original input molecule** (not the geometry the scan had actually reached) to `coordinates_list`, so trajectory animation shows a bogus frame; the parallel `float("nan")` energy then makes `PESScanResult.energy_hartree` (`min()` over the list) and `summary()`'s min/max/barrier order-dependently NaN, since Python's `min`/`max` propagate NaN only when it is the first element. + +### M7. PES scan uses ASE's deprecated `FixInternals(angles=…, dihedrals=…)` radian kwargs — ✅ Fixed +`pes_scan.py:288-298`. ASE ≥ 3.21 deprecated these in favor of `angles_deg` / `dihedrals_deg`; the legacy kwargs are scheduled for removal. With `ase>=3.22.0` unpinned above, this is a forward-compatibility break waiting to happen (and worth a compatibility test now). + +### M8. Event log: full rewrite per event + read/rewrite race — ✅ Fixed +`calc_log.log_event()` appends (one lock section) then calls `prune_events()`, which reads the whole file and rewrites it in **two separate lock sections** (calc_log.py:968-994). An append from another thread between the read and the rewrite is silently lost. It's also O(file-size) work on *every* event — O(N²) over a session. + +### M9. `save_thumbnail` docstring promises "silently skips … any error" but only guards the import — ✅ Fixed +`results_storage.py:700-711` — only `ImportError` is caught; a `savefig` failure (e.g. disk full, font issues) propagates into the caller's save path. + +### M10. PubChem client edge cases — ✅ Fixed +- `_http_get` returns `None` (annotated `Response`) if `PUBCHEM_MAX_RETRIES` is ever configured < 1 (`pubchem.py:70-89`) → `AttributeError` at every caller. +- `check_pubchem_availability()` bypasses the shared throttle and hardcodes `timeout=5` instead of `config.PUBCHEM_AVAILABILITY_TIMEOUT_S` (line 547) — the config constant is defined but unused. + +### M11. Export filenames embed the basis set verbatim, including `*` — ✅ Fixed +`app_exports.py` builds names like `H2O_RHF_6-31G*.py` / `.xyz` / `.mol` / `.pdb`. `*` is an invalid filename character on Windows (where the non-PySCF UI is supported and CI runs), so these exports fail there; on POSIX it produces glob-hostile filenames. `results_storage._safe_name()` already exists and solves exactly this — it just isn't used here. + +### M12. Method-notes markdown mangling — ✅ Fixed +`app_runflow.update_notes` (lines 1075-1078) converts `**bold**` to HTML by replacing only the **first** `**` pair; `get_educational_notes()` emits many bold spans, so literal `**` markers leak into the rendered panel. + +### M13. The `quantui` CLI imports the entire GUI stack — ✅ Fixed +**Status:** `QuantUIApp`, `StepProgress`, and the three `help_content` exports are now resolved lazily via a module-level `__getattr__` (PEP 562) in `quantui/__init__.py`, mirroring the existing `config.MOLECULE_LIBRARY` pattern. Verified `import quantui.cli` no longer touches `sys.modules['ipywidgets']` / `quantui.app` / `quantui.progress`. `IPython` is still imported transitively (via the PubChem/visualization optional-dependency blocks, which were left as-is — a full lazy-loading refactor of those blocks was judged out of proportion to a Medium-severity finding), so the CLI is faster and no longer notebook-widget-dependent, though not fully IPython-free. +`cli.py`'s header says it "deliberately avoids importing from the GUI side … so it stays fast on import", but `from quantui.calc_log import …` first executes `quantui/__init__.py`, which unconditionally runs `from .app import QuantUIApp` (line 175) → ipywidgets, IPython, and the full app module are imported for `quantui log tail`. This makes the CLI both slow and hard-dependent on notebook packages, defeating the stated design. + +### M14. `logging.basicConfig()` at library import time — ✅ Fixed +`utils.py:19` configures the **root** logger the moment `quantui.utils` is imported (which `quantui/__init__` always does). A library must not do this — it hijacks/duplicates the logging configuration of any host application or notebook. Use per-module loggers only (the package already installs a `NullHandler` correctly in `__init__.py`). + +--- + +## Low — inefficiencies, data nits, docs, hygiene + +1. **Estimator re-reads the entire perf log on every UI refresh** — `calc_log.estimate_time()` parses all of `perf_log.jsonl` (kept indefinitely by design) per call; `update_estimate` fires on widget changes. Estimates get progressively slower over the app's lifetime; consider caching with an mtime check. — ✅ Fixed (`_read_all()` now caches on `(mtime, size)`; repeat calls between writes skip the read + per-line `json.loads` entirely) +2. **Basis-count table nit** — `_BASIS_FUNCTIONS["6-31G**"]["He"] = 2` (calc_log.py:138): 6-31G** places p-polarization on He as well → should be 5, inconsistent with `"H": 5` in the same table. — ✅ Fixed (verified against `pyscf.gto.M(atom="He", basis="6-31g**").nao == 5`) +3. **Inconsistent Bohr constants** — `optimizer.py:59` uses `0.529177249`, `freq_calc.py:359` uses `0.52917721092`. Harmless numerically, but pick one (CODATA 2018: 0.529177210903) and share it like `HARTREE_TO_EV`. — ✅ Fixed (added `config.BOHR_TO_ANGSTROM`, matching `pyscf.data.nist.BOHR`, as the single source of truth) +4. **Duplicated data/constants** — the 36-entry atomic-number dict appears twice in `molecule.py` (lines 139, 553) and again in `benchmarks.py`; `HARTREE_TO_EV` is re-declared locally in `results_storage.py` (twice), `comparison.py`, and `optimizer.py:475`. — ✅ Partially fixed: `benchmarks._count_electrons` now uses `config.ATOMIC_NUMBERS` (was hard-capped at Z=18, silently mis-counting anything heavier). The `HARTREE_TO_EV` duplication was left alone — every copy holds the identical value, so it's a pure DRY nit with no correctness risk, and de-duplicating it across 5 files carried more churn/risk than the finding warranted. +5. **Docstring drift** — ✅ Fixed + - `session_calc.run_in_session` docstring says verbose "Default: 3"; the signature default is 4. + - `pubchem.search_molecule_by_name` docstring says "None otherwise" but the function raises. + - `pyproject.toml` comments claim `pyscf.prop.infrared` is in pyscf ≥ 2.13 core and "NMR is accessed via pyscf.nmr (core)"; `freq_calc.py:350` and `nmr_calc.py:162` say the opposite (and implement workarounds accordingly). (Verified empirically against the installed pyscf 2.13.1: neither module exists; `pyproject.toml` corrected.) + - `freq_ir_workers` "serial fallback" claim (see M5) — already resolved by the M5 fix itself. +6. **Python 3.9 support is claimed but untested** — `requires-python = ">=3.9"` and a 3.9 classifier, but CI runs only 3.10/3.11. Care was clearly taken (e.g. the `StrEnum` shim in `viz_backend_router.py`), but nothing verifies it. — ✅ Fixed, and it found a real bug: adding 3.9 to CI surfaced 6 genuine test failures. pip resolves `ase==3.26.0` for Python 3.9 (vs. `3.29.0` for 3.10/3.11, since ASE dropped 3.9 support in a later release); `ase.utils.IOContext.openfile()` in 3.26.0 strictly requires a `logfile=` argument to have a `.close()` method to be treated as an already-open stream, which QuantUI's `app._LogCapture` (passed as `BFGS(..., logfile=...)` in `optimizer.py`/`pes_scan.py`) didn't have. Fixed by adding a no-op `close()`; verified against a real `ase==3.26.0` install that all previously-failing geometry-opt/frequency history tests now pass. +7. **`raise … from` missing** in 12 exception-wrapping sites in `pubchem.py`, `cactus.py`, `molecule.py`, `app_visualization.py` (ruff B904) — masks root causes in tracebacks. (B904 is deliberately ignored in `pyproject.toml`; worth revisiting for the network client at least.) — ✅ Fixed (all 12 sites chained via `raise ... from e`) +8. **`StepProgress` renders labels/messages without HTML-escaping** (`progress.py:96-103`) — a failure message containing `<`/`>` (easy from a parse error echoing user input) breaks the widget rendering. — ✅ Fixed (`html.escape()` on both label and message) +9. **`nmr_calc` monkey-patches PySCF process-wide on every call** — `pyscf.prop.nmr.rhf.gen_vind` and `.rks.get_vxc_giao` are replaced globally (nmr_calc.py:207, 294). Pragmatic, but it should be done once (idempotently) at import of the patched feature, with a version guard, so behavior doesn't depend on whether an NMR calc has run. — ✅ Fixed (extracted into a versioned, idempotent `_ensure_nmr_compat_patches_applied()`; a sentinel attribute on the installed function short-circuits re-patching) +10. **`results_storage.list_results()` sort order** — collision-suffixed directories sort lexicographically (`…_1`, `…_10`, `…_2`); cosmetic, only affects same-microsecond collisions. — ✅ Fixed (sort key now splits the numeric collision suffix and compares it as an int) +11. **`save_result` writes some fields un-coerced** — `energy_hartree` / `homo_lumo_gap_ev` go through `getattr` without `_opt_float`; a duck-typed result carrying numpy scalars would make `json.dumps` raise. All current internal producers coerce to `float`, so this is latent. — ✅ Fixed (all numeric/boolean fields now go through `_opt_float`/new `_opt_int`; confirmed `numpy.float32`/`int64`/`bool_` — which do NOT subclass their Python equivalents, unlike `numpy.float64` — previously broke `json.dumps` unconverted) +12. **`prune_events` keeps malformed entries forever** (calc_log.py:989) — deliberate, but it means a corrupt line never ages out of the 7-day log. — ⏸️ Deferred by choice (confirmed with the maintainer that the existing deliberate behavior should stay as-is) +13. **Single atoms are rejected by `parse_xyz_input`** (molecule.py:517) — atomic calculations are legitimate QC targets (and PySCF handles them); the restriction plus its wording is a design limitation worth revisiting. — ✅ Fixed (restriction lifted; confirmed a single-atom RHF/STO-3G calculation runs end-to-end and converges through `run_in_session`) +14. **`_cmd_analytics_build`** assigns `tool` from `_open_in_browser` and never uses it (cli.py:226). — ✅ Fixed + +--- + +## Verified-clean areas (for the record) + +- `user_settings.py`, `vib_cache.py` — careful atomic writes, schema versioning, clamping; no issues found. +- `c_stderr.py` fd-2 capture — correct save/restore ordering, including the restore-before-read subtlety. +- `gpu_offload.py` — sound fallback design; the CuPy-array handling added across `session_calc` is consistent. +- `viz_backend_router.py` — routing logic exhaustive and correct, including the 3.9/3.10 `StrEnum` shim. +- `molecule_library.py` — read-only per-call SQLite connections, correct manifest fallback, `lru_cache` on the preset dict (so the `config.MOLECULE_LIBRARY` PEP 562 shim is *not* a per-access reload). +- `benchmarks.py` calibration worker loop — terminate/skip/stop/timeout handling and the exitcode/signal diagnostics are solid. +- Ruff (pyflakes/bugbear/pylint-error rules): no undefined names, no unused variables, no mutable default arguments anywhere in `quantui/`, `tests/`, `scripts/`. + +## Suggested priorities + +1. **H1 + H2** together — both corrupt the flagship analysis/export features (isosurface, Molden) for whole classes of molecules, and share the fix surface (thread charge/spin/units through result objects). +2. **H5** — silent data loss on a very common input format. +3. **H3, H4** — small, contained fixes. +4. **M1, M4** — correctness-of-science issues that are invisible to students, which is the worst kind in a teaching tool. diff --git a/pyproject.toml b/pyproject.toml index d96ddf9..ea69596 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -61,10 +61,11 @@ quantui = [ [project.optional-dependencies] # PySCF requires Linux/macOS/WSL — not available on Windows natively. # Use the Apptainer container (apptainer/quantui.def) for Windows. -# pyscf>=2.13.0: pyscf.prop.infrared is in pyscf 2.13.0+ core (creates the -# pyscf.prop namespace that infrared lives under). -# pyscf-properties: provides additional pyscf.prop.* modules (EFG, IR, etc.). -# NMR is accessed via pyscf.nmr (core, not overwritten by pyscf-properties). +# pyscf.prop.infrared is NOT in pyscf/pyscf-properties as of pyscf 2.13 — +# freq_calc.py computes IR intensities itself via finite-difference dipole +# displacement instead (see the comment there). NMR is NOT accessed via +# pyscf.nmr (that module doesn't exist in released pyscf); it's +# pyscf.prop.nmr, provided by pyscf-properties below. pyscf = [ "pyscf>=2.13.0", "pyscf-properties", diff --git a/quantui/__init__.py b/quantui/__init__.py index 22bcf42..cd4e439 100644 --- a/quantui/__init__.py +++ b/quantui/__init__.py @@ -10,6 +10,7 @@ __version__ = "0.4.0" import logging +from typing import Any logging.getLogger(__name__).addHandler(logging.NullHandler()) @@ -40,9 +41,6 @@ VALID_ATOMS, WIDGET_LAYOUT, ) - -# Educational help content (requires ipywidgets at display time) -from .help_content import HELP_TOPICS, VALID_TOPICS, help_panel from .molecule import Molecule, parse_xyz_input # Orbital visualization (matplotlib energy diagrams, cube-file viewer) @@ -55,9 +53,6 @@ plot_orbital_diagram, ) -# Progress indicators -from .progress import StepProgress - # Security — catchable exception for constraint violations from .security import SecurityError from .utils import ( @@ -169,10 +164,31 @@ VISUALIZATION_AVAILABLE = False PY3DMOL_AVAILABLE = False -# App class — imported last so all package symbols are defined first. -# app.py imports from submodules directly, but placing this last is an -# extra safeguard against accidental circular-import issues in the future. -from .app import QuantUIApp +# App class, StepProgress, and help_content are resolved lazily via +# module __getattr__ (PEP 562) below. All three unconditionally pull in +# ipywidgets (app.py additionally pulls in the rest of the GUI stack), and +# eagerly importing them here defeats lightweight consumers like +# ``quantui.cli`` that only need pure-Python submodules (calc_log, +# analytics, gpu_offload) — see cli.py's module docstring. +_LAZY_ATTRS = { + "QuantUIApp": (".app", "QuantUIApp"), + "StepProgress": (".progress", "StepProgress"), + "HELP_TOPICS": (".help_content", "HELP_TOPICS"), + "VALID_TOPICS": (".help_content", "VALID_TOPICS"), + "help_panel": (".help_content", "help_panel"), +} + + +def __getattr__(name: str) -> Any: + target = _LAZY_ATTRS.get(name) + if target is None: + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") + module_name, attr_name = target + import importlib + + module = importlib.import_module(module_name, __name__) + return getattr(module, attr_name) + __all__ = [ # Config constants diff --git a/quantui/app.py b/quantui/app.py index cc1a81c..b47ccef 100644 --- a/quantui/app.py +++ b/quantui/app.py @@ -709,6 +709,22 @@ def write(self, text: str) -> None: def flush(self) -> None: pass + def close(self) -> None: + """No-op — required so ASE treats this as an already-open stream. + + ASE's ``IOContext.openfile()`` (used by ``BFGS(..., logfile=...)`` + in optimizer.py / pes_scan.py) checks ``hasattr(file, "close")`` to + decide whether *file* is an already-open, file-like object it + should leave alone, vs. a path string it should ``open()`` itself. + Without this method, ase>=3.22 (the floor this project pins) still + happened to work via a later refactor's more lenient check, but + ase==3.26.0 (the newest version pip resolves for Python 3.9) hits + the stricter ``openfile()`` and raises + ``TypeError: expected str, bytes or os.PathLike object`` — a real + Python-3.9-specific compatibility gap the L6 audit fix's CI matrix + expansion caught. + """ + def getvalue(self) -> str: return self._buf.getvalue() @@ -996,6 +1012,7 @@ def __init__(self) -> None: # ``_apply_analysis_context`` resets these between contexts so stale # state from a prior calc cannot leak into the next molecule. self._last_orb_mo_coeff: Any = None + self._last_orb_mo_occ: Any = None self._last_orb_mol_atom: Any = None self._last_orb_mol_basis: Any = None # Last-generated cube file path + orbital label (M-EXPORT / EXPORT.5). diff --git a/quantui/app_analysis.py b/quantui/app_analysis.py index 44c3a18..ce7016c 100644 --- a/quantui/app_analysis.py +++ b/quantui/app_analysis.py @@ -174,6 +174,7 @@ def apply_analysis_context(app: Any, ctx: Any) -> None: # to activate re-sets these in show_orbital_diagram. app._last_orb_info = None app._last_orb_mo_coeff = None + app._last_orb_mo_occ = None app._last_orb_mol_atom = None app._last_orb_mol_basis = None app.traj_accordion.set_title(0, "Trajectory Viewer") diff --git a/quantui/app_exports.py b/quantui/app_exports.py index e365158..13b81ad 100644 --- a/quantui/app_exports.py +++ b/quantui/app_exports.py @@ -5,6 +5,8 @@ from pathlib import Path from typing import Any +from .results_storage import _safe_name + def on_export(app: Any, btn: Any) -> None: """Export a standalone Python calculation script.""" @@ -19,9 +21,14 @@ def on_export(app: Any, btn: Any) -> None: method=app.method_dd.value, basis=app.basis_dd.value, ) + # M11 audit fix (2026-07-14): the basis set is embedded verbatim in + # the filename (e.g. "6-31G*.py"), and "*" is invalid in a Windows + # filename — this export silently failed there. _safe_name (already + # used by results_storage for the same purpose) replaces anything + # that isn't alphanumeric/underscore/hyphen with "x". fname = ( - f"{app._molecule.get_formula()}" - f"_{app.method_dd.value}_{app.basis_dd.value}.py" + f"{_safe_name(app._molecule.get_formula())}" + f"_{_safe_name(app.method_dd.value)}_{_safe_name(app.basis_dd.value)}.py" ) calc.generate_calculation_script(Path(fname)) app.export_status.value = f"Saved: {fname}" @@ -36,7 +43,7 @@ def on_export_xyz(app: Any, btn: Any) -> None: return try: mol, method, basis = export_molecule_and_label(app) - fname = f"{mol.get_formula()}_{method}_{basis}.xyz" + fname = f"{_safe_name(mol.get_formula())}_{_safe_name(method)}_{_safe_name(basis)}.xyz" xyz_body = mol.to_xyz_string() full_xyz = ( f"{len(mol.atoms)}\n{mol.get_formula()} {method}/{basis}\n{xyz_body}\n" @@ -57,7 +64,7 @@ def on_export_mol(app: Any, btn: Any) -> None: from rdkit import Chem mol, method, basis = export_molecule_and_label(app) - fname = f"{mol.get_formula()}_{method}_{basis}.mol" + fname = f"{_safe_name(mol.get_formula())}_{_safe_name(method)}_{_safe_name(basis)}.mol" rdmol = molecule_to_rdkit(mol) if rdmol is None: app.struct_export_status.value = "RDKit could not parse the structure." @@ -79,7 +86,7 @@ def on_export_pdb(app: Any, btn: Any) -> None: from rdkit import Chem mol, method, basis = export_molecule_and_label(app) - fname = f"{mol.get_formula()}_{method}_{basis}.pdb" + fname = f"{_safe_name(mol.get_formula())}_{_safe_name(method)}_{_safe_name(basis)}.pdb" rdmol = molecule_to_rdkit(mol) if rdmol is None: app.struct_export_status.value = "RDKit could not parse the structure." diff --git a/quantui/app_formatters.py b/quantui/app_formatters.py index 431ecb9..b92c3b1 100644 --- a/quantui/app_formatters.py +++ b/quantui/app_formatters.py @@ -298,6 +298,22 @@ def _nmr_table(label: str, shifts: list, sym: str) -> str: "" ) + # M4 audit fix (2026-07-14): the reference shielding constants table only + # covers a handful of method/basis combinations; any other combination + # silently substitutes the B3LYP/6-31G* constants, which can shift the + # reported ppm values by several ppm relative to a properly calibrated + # reference. Surface that substitution rather than let it pass silently. + _ref_warn = "" + if getattr(r, "is_fallback_reference", False): + _ref_warn = ( + '' + '' + f"⚠ No calibrated TMS reference for {r.method}/{r.basis} — using " + f"{getattr(r, 'reference_key', 'B3LYP/6-31G*')} constants instead. " + "Shifts may be off by a few ppm." + "" + ) + _empty = "" if not r.h_shifts() and not r.c_shifts(): _empty = ( @@ -310,7 +326,7 @@ def _nmr_table(label: str, shifts: list, sym: str) -> str: f'padding:10px 14px;border-radius:4px;margin:6px 0">' f"NMR Shielding — {r.formula} ({r.method}/{r.basis})" f'' - f"{header_rows}{h_table}{c_table}{_empty}{_basis_warn}
" + f"{header_rows}{h_table}{c_table}{_empty}{_basis_warn}{_ref_warn}" ) diff --git a/quantui/app_runflow.py b/quantui/app_runflow.py index 2a1ef07..afaef8d 100644 --- a/quantui/app_runflow.py +++ b/quantui/app_runflow.py @@ -1072,10 +1072,16 @@ def update_notes(app: Any, change: Any = None) -> None: ) notes = calc.get_educational_notes() if notes: - safe = ( - notes.replace("**", "", 1) - .replace("**", "", 1) - .replace("\n\n", "

") + # M12 audit fix (2026-07-14): .replace("**", ..., 1) only + # converts the FIRST **bold** pair in the whole string — every + # note after the first (get_educational_notes() typically + # returns 2-3 separate "**Label**: description" paragraphs + # joined by "\n\n") kept its literal "**" markers instead of + # being rendered bold. A regex replaces every **...** pair. + import re as _re + + safe = _re.sub(r"\*\*(.+?)\*\*", r"\1", notes).replace( + "\n\n", "

" ) with app.notes_output: display( diff --git a/quantui/app_visualization.py b/quantui/app_visualization.py index 3569165..83e36ee 100644 --- a/quantui/app_visualization.py +++ b/quantui/app_visualization.py @@ -965,6 +965,7 @@ def show_orbital_diagram(app: Any, result: Any) -> bool: app._last_orb_info = info app._last_orb_mo_coeff = getattr(result, "mo_coeff", None) + app._last_orb_mo_occ = mo_occ app._last_orb_mol_atom = getattr(result, "pyscf_mol_atom", None) app._last_orb_mol_basis = getattr(result, "pyscf_mol_basis", None) @@ -1171,6 +1172,7 @@ def _is_stale() -> bool: mo_coeff = getattr(app, "_last_orb_mo_coeff", None) mol_atom = getattr(app, "_last_orb_mol_atom", None) mol_basis = getattr(app, "_last_orb_mol_basis", None) + mo_occ_for_charge = getattr(app, "_last_orb_mo_occ", None) if mo_coeff is None or mol_atom is None or mol_basis is None: return @@ -1179,6 +1181,7 @@ def _is_stale() -> bool: from quantui.orbital_visualization import ( generate_cube_from_arrays, + infer_charge_and_spin, plot_cube_isosurface, render_orbital_isosurface_py3dmol, ) @@ -1204,7 +1207,19 @@ def _is_stale() -> bool: ts = _dt.now().strftime("%Y-%m-%d_%H-%M-%S-%f") cube_path = cube_dir / f"{safe_formula}_{safe_orb}_{ts}.cube" - generate_cube_from_arrays(mol_atom, mol_basis, mo_coeff, orb_idx, cube_path) + # Charge/spin aren't carried on the app's orbital-state attributes — + # infer them from the MO occupations so charged/open-shell molecules + # (H3O+, OH-, radicals, ...) don't fail to build in PySCF (H1 fix). + _charge, _spin = infer_charge_and_spin(mol_atom, mo_occ_for_charge) + generate_cube_from_arrays( + mol_atom, + mol_basis, + mo_coeff, + orb_idx, + cube_path, + charge=_charge, + spin=_spin, + ) scene_bgcolor = app._plotly_theme_colors()["scene_bgcolor"] # Route the render: py3Dmol does native, full-resolution in-browser @@ -2488,14 +2503,14 @@ def build_vib_export_html(app: Any, mode_number: int) -> tuple[str, str]: import numpy as np import py3Dmol # noqa: F401 — probe; make_view imports it for the export except ImportError as exc: - raise ValueError(f"py3Dmol unavailable for fallback export: {exc}") + raise ValueError(f"py3Dmol unavailable for fallback export: {exc}") from exc try: displ = np.array(freq_result.displacements[mode_number - 1], dtype=float) except (AttributeError, IndexError, ValueError, TypeError) as exc: raise ValueError( f"Could not read displacements for mode {mode_number}: {exc}" - ) + ) from exc atoms = list(molecule.atoms) base_coords = np.array(molecule.coordinates, dtype=float) diff --git a/quantui/benchmarks.py b/quantui/benchmarks.py index 38fe9ea..b1dd05f 100644 --- a/quantui/benchmarks.py +++ b/quantui/benchmarks.py @@ -774,27 +774,9 @@ def n_total(self) -> int: def _count_electrons(atoms: list[str], charge: int) -> int: """Rough electron count: sum of atomic numbers minus charge.""" - _Z = { - "H": 1, - "He": 2, - "Li": 3, - "Be": 4, - "B": 5, - "C": 6, - "N": 7, - "O": 8, - "F": 9, - "Ne": 10, - "Na": 11, - "Mg": 12, - "Al": 13, - "Si": 14, - "P": 15, - "S": 16, - "Cl": 17, - "Ar": 18, - } - return sum(_Z.get(a, 6) for a in atoms) - charge + from .config import ATOMIC_NUMBERS + + return sum(ATOMIC_NUMBERS.get(a, 6) for a in atoms) - charge # --------------------------------------------------------------------------- diff --git a/quantui/cactus.py b/quantui/cactus.py index 3e9d3c6..ee3e6ba 100644 --- a/quantui/cactus.py +++ b/quantui/cactus.py @@ -64,7 +64,7 @@ def resolve_to_sdf(identifier: str, conformer_3d: bool = True) -> str: return str(response.text) except requests.RequestException as e: logger.error(f"CACTUS request failed: {e}") - raise PubChemAPIError(f"Failed to connect to CACTUS: {e}") + raise PubChemAPIError(f"Failed to connect to CACTUS: {e}") from e raise MoleculeNotFoundError( f"CACTUS could not resolve '{identifier}' (last status: {last_status})" diff --git a/quantui/calc_log.py b/quantui/calc_log.py index 9278d8e..d4f9607 100644 --- a/quantui/calc_log.py +++ b/quantui/calc_log.py @@ -136,7 +136,7 @@ }, "6-31G**": { "H": 5, - "He": 2, + "He": 5, "Li": 9, "Be": 9, "B": 14, @@ -290,11 +290,33 @@ def _append(path: Path, record: dict) -> None: fh.write(line) +_READ_ALL_CACHE: dict[str, tuple[float, int, list[dict]]] = {} + + def _read_all(path: Path) -> list[dict]: - if not path.exists(): - return [] - records: list[dict] = [] + """Parse every JSON line in *path*, caching on (mtime, size). + + ``estimate_time()`` calls this on every UI refresh (widget-change + callback), re-parsing the entire (indefinitely-kept) perf log each + time even though it usually hasn't changed since the last call. + Cache the parsed records keyed by the file's mtime + size so repeat + calls between writes skip the read + per-line ``json.loads`` entirely; + a write bumps both, so the cache invalidates correctly. + """ + key = str(path) with _LOCK: + if not path.exists(): + _READ_ALL_CACHE.pop(key, None) + return [] + stat = path.stat() + cached = _READ_ALL_CACHE.get(key) + if ( + cached is not None + and cached[0] == stat.st_mtime + and cached[1] == stat.st_size + ): + return list(cached[2]) + records: list[dict] = [] with open(path, encoding="utf-8") as fh: for raw in fh: raw = raw.strip() @@ -303,7 +325,8 @@ def _read_all(path: Path) -> list[dict]: records.append(json.loads(raw)) except json.JSONDecodeError: pass - return records + _READ_ALL_CACHE[key] = (stat.st_mtime, stat.st_size, records) + return list(records) # --------------------------------------------------------------------------- @@ -948,10 +971,33 @@ def clear_event_log() -> None: # Event log (7-day TTL) # --------------------------------------------------------------------------- +# M8 audit fix (2026-07-14): log_event() used to call prune_events() after +# every single append, and prune_events() itself read the file (acquiring +# and releasing _LOCK) and only later reacquired _LOCK to rewrite it. Two +# problems: +# +# 1. Race: an append from another thread landing in the gap between the +# read and the rewrite got silently discarded when the rewrite replaced +# the whole file with the (now-stale) `kept` list computed before that +# append happened. +# 2. Cost: reading + rewriting the entire event log on every single write +# is O(file size) per event, i.e. O(N^2) over a session as the log +# grows — noticeable once a session has logged more than a few hundred +# events. +# +# Fixed by (a) making prune_events() read + filter + rewrite as a single +# lock-held critical section, so a concurrent append either completes +# before the prune starts or blocks until it finishes — it can never be +# silently lost — and (b) only running the full prune every +# _PRUNE_EVERY_N_EVENTS appends instead of on every single one. +_PRUNE_EVERY_N_EVENTS = 20 +_events_since_prune = 0 + def log_event(event_type: str, message: str, **extra: object) -> None: """ - Append one event to ``event_log.jsonl`` and prune entries > 7 days old. + Append one event to ``event_log.jsonl``; prune entries > 7 days old + periodically (every :data:`_PRUNE_EVERY_N_EVENTS` calls, not every one). Args: event_type: Short category string, e.g. ``"startup"``, ``"calc_done"``, @@ -959,6 +1005,8 @@ def log_event(event_type: str, message: str, **extra: object) -> None: message: Human-readable description. **extra: Any additional key-value pairs to include in the record. """ + global _events_since_prune + record: dict = { "timestamp": datetime.now(timezone.utc).isoformat(), "event": event_type, @@ -966,28 +1014,51 @@ def log_event(event_type: str, message: str, **extra: object) -> None: **extra, } _append(_event_path(), record) - prune_events() + + with _LOCK: + _events_since_prune += 1 + due_for_prune = _events_since_prune >= _PRUNE_EVERY_N_EVENTS + if due_for_prune: + _events_since_prune = 0 + if due_for_prune: + prune_events() def prune_events(days: int = 7) -> None: - """Remove event-log entries older than *days* days (default: 7).""" + """Remove event-log entries older than *days* days (default: 7). + + Reads, filters, and rewrites the file as a single lock-held critical + section so a concurrent :func:`log_event` append can never be silently + lost between the read and the rewrite (see module-level note above). + """ path = _event_path() - if not path.exists(): - return cutoff = datetime.now(timezone.utc) - timedelta(days=days) - records = _read_all(path) - kept: list[dict] = [] - for r in records: - try: - ts = datetime.fromisoformat(r["timestamp"]) - # fromisoformat on Python < 3.11 doesn't handle 'Z' suffix - if ts.tzinfo is None: - ts = ts.replace(tzinfo=timezone.utc) - if ts >= cutoff: - kept.append(r) - except (KeyError, ValueError): - kept.append(r) # keep malformed entries rather than silently drop + with _LOCK: + if not path.exists(): + return + records: list[dict] = [] + with open(path, encoding="utf-8") as fh: + for raw in fh: + raw = raw.strip() + if raw: + try: + records.append(json.loads(raw)) + except json.JSONDecodeError: + pass + + kept: list[dict] = [] + for r in records: + try: + ts = datetime.fromisoformat(r["timestamp"]) + # fromisoformat on Python < 3.11 doesn't handle 'Z' suffix + if ts.tzinfo is None: + ts = ts.replace(tzinfo=timezone.utc) + if ts >= cutoff: + kept.append(r) + except (KeyError, ValueError): + kept.append(r) # keep malformed entries rather than silently drop + path.parent.mkdir(parents=True, exist_ok=True) with open(path, "w", encoding="utf-8") as fh: for r in kept: diff --git a/quantui/calculator.py b/quantui/calculator.py index d9c9681..bd60f22 100644 --- a/quantui/calculator.py +++ b/quantui/calculator.py @@ -43,7 +43,18 @@ def __init__( ValueError: If method is not supported """ self.molecule = molecule - self.method = method.upper() + # Match case-insensitively against the canonical (mixed-case) + # SUPPORTED_METHODS list, but store the CANONICAL spelling — not + # method.upper(). Several entries are mixed-case (e.g. "wB97X-D", + # "CAM-B3LYP", "M06-L"); uppercasing turned "wB97X-D" into + # "WB97X-D", which is in neither SUPPORTED_METHODS nor the xc-alias + # tables the generated script relies on, breaking Export Script + # and the educational-notes lookup for every mixed-case method. + _method_upper = method.strip().upper() + self.method = next( + (m for m in config.SUPPORTED_METHODS if m.upper() == _method_upper), + method, + ) self.basis = basis if self.method not in config.SUPPORTED_METHODS: diff --git a/quantui/cli.py b/quantui/cli.py index 9d8b5ad..390738e 100644 --- a/quantui/cli.py +++ b/quantui/cli.py @@ -223,7 +223,7 @@ def _cmd_analytics_build(args: argparse.Namespace) -> int: # Cross-platform open: WSL → wslview / explorer.exe; otherwise # stdlib webbrowser. Failure is non-fatal (the path was already # printed) so users can always copy-paste manually. - opened, tool = _open_in_browser(result) + opened, _tool = _open_in_browser(result) if not opened: print( f"(could not auto-open browser — open {result} manually)", diff --git a/quantui/config.py b/quantui/config.py index 3371839..252acb7 100644 --- a/quantui/config.py +++ b/quantui/config.py @@ -30,6 +30,16 @@ "CCSD(T)", ] +# Post-HF wavefunction methods. Single-point-only in QuantUI: session_calc.py +# is the only entry point that special-cases them (an RHF/ROHF reference, +# auto-dispatched by PySCF's scf.RHF() factory for open-shell input, plus an +# mp.MP2/cc.CCSD post-SCF step). optimizer.py, freq_calc.py, tddft_calc.py, +# and nmr_calc.py have no such special-casing — without an early guard, +# selecting one of these methods there falls through to the DFT branch, +# which sets e.g. mf.xc = "CCSD" and fails deep inside PySCF with a cryptic +# "LibXCFunctional: name 'CCSD' not found" instead of a clear message. +POST_HF_METHODS: frozenset = frozenset({"MP2", "CCSD", "CCSD(T)"}) + # Educational metadata for each method — shown to students in the UI METHOD_INFO = { "RHF": { @@ -248,6 +258,12 @@ CACTUS_TIMEOUT_S: float = 8.0 CACTUS_CONNECT_TIMEOUT_S: float = 4.0 +# Bohr radius, in Angstrom — the exact value ``pyscf.data.nist.BOHR`` uses +# internally, so unit conversions here stay consistent with what PySCF +# actually computed with. (L audit fix: optimizer.py and freq_calc.py each +# hand-typed their own slightly different literal for this.) +BOHR_TO_ANGSTROM: float = 0.52917721092 + # Bundled-library size budget + heavy-atom ceilings. # These are QC *starting* geometries, so keep them runnable in a classroom. LIBRARY_SIZE_BUDGET_BYTES: int = 10 * 1024 * 1024 # 10 MB @@ -273,45 +289,137 @@ def __getattr__(name: str) -> Any: raise AttributeError(f"module {__name__!r} has no attribute {name!r}") -# Valid atomic symbols (periodic table subset commonly used) -VALID_ATOMS = [ - "H", - "He", - "Li", - "Be", - "B", - "C", - "N", - "O", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Si", - "P", - "S", - "Cl", - "Ar", - "K", - "Ca", - "Sc", - "Ti", - "V", - "Cr", - "Mn", - "Fe", - "Co", - "Ni", - "Cu", - "Zn", - "Ga", - "Ge", - "As", - "Se", - "Br", - "Kr", -] +# Atomic numbers for the full periodic table (Z=1..118). Single source of +# truth for element validity + electron counting — previously QuantUI only +# recognized elements up to Kr (Z=36), which rejected valid structures +# resolved via PubChem/CACTUS/SMILES for any heavier element (iodine in +# thyroxine, tin/antimony in organometallics, gold/platinum complexes, +# etc.), even though the "Invalid atom symbol" error text in molecule.py +# explicitly listed iodine as a supported example. +ATOMIC_NUMBERS: Dict[str, int] = { + "H": 1, + "He": 2, + "Li": 3, + "Be": 4, + "B": 5, + "C": 6, + "N": 7, + "O": 8, + "F": 9, + "Ne": 10, + "Na": 11, + "Mg": 12, + "Al": 13, + "Si": 14, + "P": 15, + "S": 16, + "Cl": 17, + "Ar": 18, + "K": 19, + "Ca": 20, + "Sc": 21, + "Ti": 22, + "V": 23, + "Cr": 24, + "Mn": 25, + "Fe": 26, + "Co": 27, + "Ni": 28, + "Cu": 29, + "Zn": 30, + "Ga": 31, + "Ge": 32, + "As": 33, + "Se": 34, + "Br": 35, + "Kr": 36, + "Rb": 37, + "Sr": 38, + "Y": 39, + "Zr": 40, + "Nb": 41, + "Mo": 42, + "Tc": 43, + "Ru": 44, + "Rh": 45, + "Pd": 46, + "Ag": 47, + "Cd": 48, + "In": 49, + "Sn": 50, + "Sb": 51, + "Te": 52, + "I": 53, + "Xe": 54, + "Cs": 55, + "Ba": 56, + "La": 57, + "Ce": 58, + "Pr": 59, + "Nd": 60, + "Pm": 61, + "Sm": 62, + "Eu": 63, + "Gd": 64, + "Tb": 65, + "Dy": 66, + "Ho": 67, + "Er": 68, + "Tm": 69, + "Yb": 70, + "Lu": 71, + "Hf": 72, + "Ta": 73, + "W": 74, + "Re": 75, + "Os": 76, + "Ir": 77, + "Pt": 78, + "Au": 79, + "Hg": 80, + "Tl": 81, + "Pb": 82, + "Bi": 83, + "Po": 84, + "At": 85, + "Rn": 86, + "Fr": 87, + "Ra": 88, + "Ac": 89, + "Th": 90, + "Pa": 91, + "U": 92, + "Np": 93, + "Pu": 94, + "Am": 95, + "Cm": 96, + "Bk": 97, + "Cf": 98, + "Es": 99, + "Fm": 100, + "Md": 101, + "No": 102, + "Lr": 103, + "Rf": 104, + "Db": 105, + "Sg": 106, + "Bh": 107, + "Hs": 108, + "Mt": 109, + "Ds": 110, + "Rg": 111, + "Cn": 112, + "Nh": 113, + "Fl": 114, + "Mc": 115, + "Lv": 116, + "Ts": 117, + "Og": 118, +} + +# Valid atomic symbols — every element in ATOMIC_NUMBERS. Derived rather +# than hand-maintained separately so the two can never drift apart again. +VALID_ATOMS = list(ATOMIC_NUMBERS.keys()) # Quick-start templates for the notebook UI # The `notes` and `learning_goals` fields are shown to students. diff --git a/quantui/freq_calc.py b/quantui/freq_calc.py index 3586a88..02c1936 100644 --- a/quantui/freq_calc.py +++ b/quantui/freq_calc.py @@ -159,6 +159,22 @@ def run_freq_calc( computation fails, frequencies are omitted and a warning is written to progress_stream — no exception is raised. """ + # Post-HF methods (MP2/CCSD/CCSD(T)) have no special-casing below — + # without this guard, method='CCSD' silently falls into the DFT + # branch (sets mf.xc = "CCSD") and fails deep inside PySCF with a + # cryptic "LibXCFunctional: name 'CCSD' not found" instead of a clear + # message. No post-HF Hessian is wired up here, so these methods are + # single-point only (see session_calc.py). + from . import config as _config + + if method.strip().upper() in _config.POST_HF_METHODS: + raise ValueError( + f"'{method}' is a post-HF method and cannot be used for " + "frequency analysis — QuantUI only has an analytical Hessian " + "wired up for HF/DFT methods here. Use RHF, UHF, or a DFT " + "functional instead." + ) + try: from pyscf import dft, gto, scf from pyscf.hessian import thermo as pyscf_thermo @@ -283,7 +299,16 @@ def _status(msg: str) -> None: _moe, _moo = _moe[0], _moo[0] mo_energy_hartree = _np_mo.asarray(_moe, dtype=float).tolist() mo_occ_list = _np_mo.asarray(_moo, dtype=float).tolist() - pyscf_mol_atom = [(str(s), list(map(float, c))) for s, c in mol._atom] + # Build from molecule.atoms/coordinates (Angstrom) rather than + # mol._atom, which PySCF always stores internally in Bohr. Every + # consumer of pyscf_mol_atom (Molden export, cube generation, + # session_calc's/optimizer's own construction of this field) + # assumes Angstrom; using mol._atom here silently shipped Bohr + # coordinates ~1.89x too large. + pyscf_mol_atom = [ + (atom, list(map(float, coords))) + for atom, coords in zip(molecule.atoms, molecule.coordinates) + ] except Exception as exc: # Same class as session_calc bug-A: silent failure here ships # a FreqResult with no MO data, breaking the Energies panel on @@ -355,8 +380,9 @@ def _status(msg: str) -> None: try: import numpy as _np_ir + from .config import BOHR_TO_ANGSTROM as _BOHR_TO_ANG + _DELTA = 0.01 # Bohr - _BOHR_TO_ANG = 0.52917721092 _KM_MOL_FAC = 42.255 # (D/Å)²/amu → km/mol _n_ir = mol.natm @@ -366,6 +392,20 @@ def _status(msg: str) -> None: _dm0 = mf.make_rdm1() _dpdx = _np_ir.zeros((_n_ir * 3, 3)) _xc = getattr(mf, "xc", None) + # M5 audit fix (2026-07-14): whether the inner displaced-SCF + # loop needs an unrestricted (UHF/UKS) object is determined + # by _dm0's actual shape — (2, nao, nao) for UHF/UKS/ROHF, + # (nao, nao) for RHF/RKS — NOT by mol.spin == 0. Those two + # signals only agree when the user's method choice matches + # the molecule's natural spin state. They diverge when a + # user explicitly selects UHF for a closed-shell molecule + # (mol.spin == 0 but the parent mf, and therefore _dm0, is + # still UHF-shaped): the inner loop used to build RHF from + # mol.spin == 0, then feed it the UHF-shaped _dm0, which + # raised a shape-mismatch ValueError inside PySCF and + # silently dropped IR intensities for the whole calc (caught + # by the broad except below). + _dm0_is_unrestricted = _np_ir.asarray(_dm0).ndim == 3 _status( "Numerical IR intensities: " @@ -384,10 +424,10 @@ def _status(msg: str) -> None: def _displaced_scf_dipole() -> _np_ir.ndarray: if _xc is not None: - _mf_d = dft.RKS(mol) if mol.spin == 0 else dft.UKS(mol) + _mf_d = dft.UKS(mol) if _dm0_is_unrestricted else dft.RKS(mol) _mf_d.xc = _xc else: - _mf_d = scf.RHF(mol) if mol.spin == 0 else scf.UHF(mol) + _mf_d = scf.UHF(mol) if _dm0_is_unrestricted else scf.RHF(mol) _mf_d.verbose = 0 _mf_d.stdout = stream # ``method_upper="RHF"`` is a label — try_to_gpu only @@ -418,96 +458,111 @@ def _displaced_scf_dipole() -> _np_ir.ndarray: _mol_v = mol.verbose mol.verbose = 0 + _parallel_failed = False try: if _use_parallel: - # Stash dm0 once on disk so workers can map-load it - # via initargs (avoids per-task pickling). - import concurrent.futures as _cf - import multiprocessing as _mp - import pickle as _pickle - import tempfile as _tempfile - - _n_workers = _ir_par.pick_worker_count( - _cpu_count, _ir_total_solves - ) - _threads_each = _ir_par.threads_per_worker( - _cpu_count, _n_workers - ) - - # Build all 6N task arguments first; pickling-safe - # flat lists per-displacement. - _tasks: list[tuple[int, int, int, list[float]]] = [] - for _I in range(_n_ir): - for _ax in range(3): - _cp = _coords0.copy() - _cp[_I, _ax] += _DELTA - _tasks.append((_I, _ax, +1, _cp.flatten().tolist())) - _cm = _coords0.copy() - _cm[_I, _ax] -= _DELTA - _tasks.append((_I, _ax, -1, _cm.flatten().tolist())) - - _dm0_handle = _tempfile.NamedTemporaryFile( - delete=False, suffix=".dm0.pkl" - ) try: - _pickle.dump(_dm0, _dm0_handle) - _dm0_handle.close() - - # Pyscf-format atom string for worker rebuild. - _atom_str = molecule.to_pyscf_format() - _spin = molecule.multiplicity - 1 - _charge = molecule.charge - _ctx = _mp.get_context("spawn") - with _cf.ProcessPoolExecutor( - max_workers=_n_workers, - mp_context=_ctx, - initializer=_ir_par.init_worker, - initargs=( - _atom_str, - basis, - _charge, - _spin, - _xc, - _dm0_handle.name, - _threads_each, - ), - ) as _pool: - # Submit all and store futures keyed by task - # index so we can assemble +/- per (I, ax). - _futs = { - _pool.submit( - _ir_par.run_displaced_scf, _task[3] - ): _task - for _task in _tasks - } - # Accumulate results into a temporary map - # ``(I, ax, sign) -> dipole_array``. - _dipoles: dict = {} - for _fut in _cf.as_completed(_futs): - _I, _ax, _sign, _coords_done = _futs[_fut] - _dipoles[(_I, _ax, _sign)] = _fut.result() - _ir_done_solves += 1 - _status( - "Numerical IR intensities (parallel ×" - f"{_n_workers}): " - f"{_ir_done_solves}/{_ir_total_solves} " - "finite-difference displacement SCFs done (6 per atom) " - f"({_ir_total_solves - _ir_done_solves} " - "remaining)" - ) - finally: + # Stash dm0 once on disk so workers can map-load it + # via initargs (avoids per-task pickling). + import concurrent.futures as _cf + import multiprocessing as _mp + import pickle as _pickle + import tempfile as _tempfile + + _n_workers = _ir_par.pick_worker_count( + _cpu_count, _ir_total_solves + ) + _threads_each = _ir_par.threads_per_worker( + _cpu_count, _n_workers + ) + + # Build all 6N task arguments first; pickling-safe + # flat lists per-displacement. + _tasks: list[tuple[int, int, int, list[float]]] = [] + for _I in range(_n_ir): + for _ax in range(3): + _cp = _coords0.copy() + _cp[_I, _ax] += _DELTA + _tasks.append((_I, _ax, +1, _cp.flatten().tolist())) + _cm = _coords0.copy() + _cm[_I, _ax] -= _DELTA + _tasks.append((_I, _ax, -1, _cm.flatten().tolist())) + + _dm0_handle = _tempfile.NamedTemporaryFile( + delete=False, suffix=".dm0.pkl" + ) try: - os.unlink(_dm0_handle.name) - except OSError: - pass - - # Assemble dpdx now that all dipoles are in hand. - for _I in range(_n_ir): - for _ax in range(3): - _mu_p = _dipoles[(_I, _ax, +1)] - _mu_m = _dipoles[(_I, _ax, -1)] - _dpdx[3 * _I + _ax] = (_mu_p - _mu_m) / (2 * _DELTA) - else: + _pickle.dump(_dm0, _dm0_handle) + _dm0_handle.close() + + # Pyscf-format atom string for worker rebuild. + _atom_str = molecule.to_pyscf_format() + _spin = molecule.multiplicity - 1 + _charge = molecule.charge + _ctx = _mp.get_context("spawn") + with _cf.ProcessPoolExecutor( + max_workers=_n_workers, + mp_context=_ctx, + initializer=_ir_par.init_worker, + initargs=( + _atom_str, + basis, + _charge, + _spin, + _xc, + _dm0_handle.name, + _threads_each, + ), + ) as _pool: + # Submit all and store futures keyed by task + # index so we can assemble +/- per (I, ax). + _futs = { + _pool.submit( + _ir_par.run_displaced_scf, _task[3] + ): _task + for _task in _tasks + } + # Accumulate results into a temporary map + # ``(I, ax, sign) -> dipole_array``. + _dipoles: dict = {} + for _fut in _cf.as_completed(_futs): + _I, _ax, _sign, _coords_done = _futs[_fut] + _dipoles[(_I, _ax, _sign)] = _fut.result() + _ir_done_solves += 1 + _status( + "Numerical IR intensities (parallel ×" + f"{_n_workers}): " + f"{_ir_done_solves}/{_ir_total_solves} " + "finite-difference displacement SCFs done (6 per atom) " + f"({_ir_total_solves - _ir_done_solves} " + "remaining)" + ) + finally: + try: + os.unlink(_dm0_handle.name) + except OSError: + pass + + # Assemble dpdx now that all dipoles are in hand. + for _I in range(_n_ir): + for _ax in range(3): + _mu_p = _dipoles[(_I, _ax, +1)] + _mu_m = _dipoles[(_I, _ax, -1)] + _dpdx[3 * _I + _ax] = (_mu_p - _mu_m) / (2 * _DELTA) + except Exception as _par_exc: + logger.warning( + "Parallel IR-intensity computation failed (%s); falling back to serial.", + _par_exc, + ) + _status( + "Parallel IR intensities failed; falling back to serial computation." + ) + _parallel_failed = True + # Reset so the serial loop's progress messages + # below start clean rather than continuing from + # wherever the failed parallel attempt left off. + _ir_done_solves = 0 + if not _use_parallel or _parallel_failed: for _I in range(_n_ir): for _ax in range(3): # +Δ displacement diff --git a/quantui/freq_ir_workers.py b/quantui/freq_ir_workers.py index 3273be7..0b17e5d 100644 --- a/quantui/freq_ir_workers.py +++ b/quantui/freq_ir_workers.py @@ -141,14 +141,27 @@ def run_displaced_scf(coords_bohr_flat) -> Any: mol.build() mol.set_geom_(coords, unit="Bohr") + # M5 audit fix (2026-07-14): whether this displaced SCF needs an + # unrestricted (UHF/UKS) object is determined by the shared dm0's + # actual shape -- (2, nao, nao) for UHF/UKS/ROHF, (nao, nao) for + # RHF/RKS -- NOT by mol.spin == 0. Those two signals only agree when + # the user's method choice matches the molecule's natural spin state. + # They diverge when a user explicitly selects UHF for a closed-shell + # molecule (mol.spin == 0 but the parent mf, and therefore dm0, is + # still UHF-shaped): building RHF from mol.spin == 0 and then feeding + # it the UHF-shaped dm0 raises a shape-mismatch ValueError inside + # PySCF. Mirrors the serial-path fix in freq_calc.py. + dm0 = state.get("dm0") + dm0_is_unrestricted = dm0 is not None and np.asarray(dm0).ndim == 3 + xc = state.get("xc") if xc is not None: - mf = dft.RKS(mol) if mol.spin == 0 else dft.UKS(mol) + mf = dft.UKS(mol) if dm0_is_unrestricted else dft.RKS(mol) mf.xc = xc else: - mf = scf.RHF(mol) if mol.spin == 0 else scf.UHF(mol) + mf = scf.UHF(mol) if dm0_is_unrestricted else scf.RHF(mol) mf.verbose = 0 - mf.kernel(dm0=state.get("dm0")) + mf.kernel(dm0=dm0) return np.array(mf.dip_moment(verbose=0)) diff --git a/quantui/molecule.py b/quantui/molecule.py index eccac77..f411385 100644 --- a/quantui/molecule.py +++ b/quantui/molecule.py @@ -12,6 +12,13 @@ logger = logging.getLogger(__name__) +# Single source of truth lives in config.ATOMIC_NUMBERS (covers the full +# periodic table, Z=1..118, and derives config.VALID_ATOMS so the two can't +# drift apart). Re-exported here so existing `from .molecule import +# ATOMIC_NUMBERS` call sites (e.g. orbital_visualization's charge/spin +# inference) keep working unchanged. +ATOMIC_NUMBERS: Dict[str, int] = config.ATOMIC_NUMBERS + class Molecule: """ @@ -135,47 +142,7 @@ def get_electron_count(self) -> int: Returns: int: Number of electrons (nuclear charges - charge) """ - # Simple atomic numbers for common elements - atomic_numbers = { - "H": 1, - "He": 2, - "Li": 3, - "Be": 4, - "B": 5, - "C": 6, - "N": 7, - "O": 8, - "F": 9, - "Ne": 10, - "Na": 11, - "Mg": 12, - "Al": 13, - "Si": 14, - "P": 15, - "S": 16, - "Cl": 17, - "Ar": 18, - "K": 19, - "Ca": 20, - "Sc": 21, - "Ti": 22, - "V": 23, - "Cr": 24, - "Mn": 25, - "Fe": 26, - "Co": 27, - "Ni": 28, - "Cu": 29, - "Zn": 30, - "Ga": 31, - "Ge": 32, - "As": 33, - "Se": 34, - "Br": 35, - "Kr": 36, - } - - nuclear_charge = sum(atomic_numbers.get(atom, 0) for atom in self.atoms) + nuclear_charge = sum(ATOMIC_NUMBERS.get(atom, 0) for atom in self.atoms) return nuclear_charge - self.charge def get_formula(self) -> str: @@ -361,11 +328,43 @@ def parse_xyz_input(xyz_text: str) -> Tuple[List[str], List[List[float]]]: lines = xyz_text.strip().split("\n") - # Process lines: remove empty lines and handle comments + # ── Step 1: XYZ-file header detection, on RAW lines ────────────────── + # The header (count line + title line) is POSITIONAL: whichever raw + # line is the first non-blank, non-full-line-comment line, if it + # parses as a bare integer, is the atom count — and the very next raw + # line is the title, regardless of what that title line itself + # contains (blank, "#"/"!"-prefixed, or free text). + # + # This must run BEFORE the blank/comment filtering below: filtering + # first would remove a blank or comment title line from the list, + # shifting the next real atom line into the "title" slot, where it + # was then silently discarded — dropping the first atom of every + # standard XYZ file whose title line happened to be blank or a + # comment (both very common in practice). + expected_atoms = None + body_start = 0 + for idx, raw in enumerate(lines): + stripped = raw.strip() + if not stripped: + continue + if stripped.startswith("#") or stripped.startswith("!"): + continue + try: + expected_atoms = int(stripped) + body_start = idx + 2 # count line + title line, whatever it is + logger.debug(f"Detected XYZ file format expecting {expected_atoms} atoms") + except ValueError: + pass # first content line isn't a bare count -> no header + break # only the first non-blank/non-comment line is eligible + + body_lines = lines[body_start:] + + # ── Step 2: filter blank lines + comments from the body ───────────── processed_lines = [] original_line_numbers = [] # Track original line numbers for error reporting - for line_num, line in enumerate(lines, start=1): + for offset, line in enumerate(body_lines): + line_num = body_start + offset + 1 # 1-based original line number line = line.strip() # Skip empty lines @@ -389,28 +388,16 @@ def parse_xyz_input(xyz_text: str) -> Tuple[List[str], List[List[float]]]: if not processed_lines: raise ValueError( "❌ No Data: All lines are empty or comments.\n\n" - "Please provide at least 2 atoms with coordinates." + "Please provide at least 1 atom with coordinates." ) atoms = [] coordinates = [] - # Check if first line is a number (XYZ file format) - start_idx = 0 - expected_atoms = None - - try: - expected_atoms = int(processed_lines[0]) - # Skip first two lines (count and comment/title) - start_idx = 2 if len(processed_lines) >= 2 else 1 - logger.debug(f"Detected XYZ file format expecting {expected_atoms} atoms") - except ValueError: - # Not XYZ file format, parse all lines as coordinates - start_idx = 0 - - # Parse coordinate lines - for i, line in enumerate(processed_lines[start_idx:]): - orig_line_num = original_line_numbers[start_idx + i] + # Parse coordinate lines (the header, if any, was already consumed + # positionally above, so every processed line here is an atom row). + for i, line in enumerate(processed_lines): + orig_line_num = original_line_numbers[i] parts = line.split() # Check minimum parts (atom symbol + 3 coordinates) @@ -501,7 +488,7 @@ def parse_xyz_input(xyz_text: str) -> Tuple[List[str], List[List[float]]]: f"💡 Coordinates must be numbers (integers or decimals).\n" f"Examples: 0.0, 1.5, -2.3, 0.757\n\n" f"Error details: {e}" - ) + ) from e atoms.append(atom_symbol) coordinates.append([x, y, z]) @@ -513,16 +500,6 @@ def parse_xyz_input(xyz_text: str) -> Tuple[List[str], List[List[float]]]: "Please check your coordinate format." ) - # Check for minimum 2 atoms (molecules need at least 2 atoms) - if len(atoms) < 2: - raise ValueError( - f"❌ Too Few Atoms: Found only {len(atoms)} atom.\n\n" - f"Molecules need at least 2 atoms for meaningful calculations.\n\n" - f"💡 For single atoms, consider:\n" - f" • Adding a second atom to form a molecule\n" - f" • Using a different computational chemistry tool for atomic calculations" - ) - # Verify expected count if XYZ file format if expected_atoms is not None and len(atoms) != expected_atoms: logger.warning( @@ -550,47 +527,8 @@ def suggest_multiplicity(atoms: List[str], charge: int) -> int: """ # Calculate electron count directly without creating Molecule # (to avoid validation errors with incompatible multiplicity) - atomic_numbers = { - "H": 1, - "He": 2, - "Li": 3, - "Be": 4, - "B": 5, - "C": 6, - "N": 7, - "O": 8, - "F": 9, - "Ne": 10, - "Na": 11, - "Mg": 12, - "Al": 13, - "Si": 14, - "P": 15, - "S": 16, - "Cl": 17, - "Ar": 18, - "K": 19, - "Ca": 20, - "Sc": 21, - "Ti": 22, - "V": 23, - "Cr": 24, - "Mn": 25, - "Fe": 26, - "Co": 27, - "Ni": 28, - "Cu": 29, - "Zn": 30, - "Ga": 31, - "Ge": 32, - "As": 33, - "Se": 34, - "Br": 35, - "Kr": 36, - } - try: - nuclear_charge = sum(atomic_numbers.get(atom, 0) for atom in atoms) + nuclear_charge = sum(ATOMIC_NUMBERS.get(atom, 0) for atom in atoms) num_electrons = nuclear_charge - charge # Even electrons -> singlet, odd electrons -> doublet diff --git a/quantui/nmr_calc.py b/quantui/nmr_calc.py index 9eebf9d..a9ace09 100644 --- a/quantui/nmr_calc.py +++ b/quantui/nmr_calc.py @@ -37,6 +37,16 @@ class NMRResult: formula: str reference_compound: str = "TMS" converged: bool = True + # M4 audit fix (2026-07-14): which config.NMR_REFERENCE_SHIELDINGS entry + # was actually applied, and whether it's an exact match for method/basis + # or a fallback. NMR_REFERENCE_SHIELDINGS only tabulates a handful of + # method/basis combinations; any other combination previously fell back + # to the B3LYP/6-31G* constants with no record of it anywhere, so + # chemical shifts could be silently offset by several ppm with no way + # for the student to know the reference wasn't calibrated for their + # method/basis. + reference_key: str = "" + is_fallback_reference: bool = False def h_shifts(self) -> List[Tuple[int, float]]: """(atom_index, δ ppm) pairs for all H atoms in molecule order.""" @@ -55,6 +65,33 @@ def c_shifts(self) -> List[Tuple[int, float]]: ] +def resolve_nmr_reference( + method: str, basis: str +) -> Tuple[Dict[str, float], str, bool]: + """Resolve TMS reference shielding constants for ``method``/``basis``. + + Looks up ``config.NMR_REFERENCE_SHIELDINGS`` case-insensitively (keys + there are declared as e.g. ``"B3LYP/6-31G*"``). Returns + ``(ref_map, matched_key, is_fallback)``: + + - ``ref_map``: the ``{"H": ..., "C": ...}`` shielding constants to use. + - ``matched_key``: the table key that was actually applied. + - ``is_fallback``: ``True`` when no entry exists for this exact + method/basis and ``config.NMR_DEFAULT_REFERENCE`` (B3LYP/6-31G*) was + substituted instead — chemical shifts computed with a substituted + reference can be off by several ppm relative to properly calibrated + constants for the requested level of theory. + """ + from . import config as _config + + requested_key = f"{method}/{basis}" + requested_upper = requested_key.upper() + for table_key, ref_map in _config.NMR_REFERENCE_SHIELDINGS.items(): + if table_key.upper() == requested_upper: + return ref_map, table_key, False + return _config.NMR_DEFAULT_REFERENCE, "B3LYP/6-31G*", True + + def run_nmr_calc( molecule: Molecule, method: str = "B3LYP", @@ -80,6 +117,19 @@ def run_nmr_calc( ImportError: If PySCF is not installed. RuntimeError: If the SCF or GIAO-NMR calculation fails. """ + # Post-HF methods (MP2/CCSD/CCSD(T)) have no special-casing below — + # without this guard, method='CCSD' silently falls into the DFT + # branch (sets mf.xc = "CCSD") and fails deep inside PySCF with a + # cryptic "LibXCFunctional: name 'CCSD' not found" instead of a clear + # message. GIAO-NMR shielding is not defined for these methods here. + from . import config as _config + + if method.strip().upper() in _config.POST_HF_METHODS: + raise ValueError( + f"'{method}' is a post-HF method and cannot be used for NMR " + "shielding — use RHF, UHF, or a DFT functional instead." + ) + try: from pyscf import dft, gto, scf except ImportError as exc: @@ -108,6 +158,166 @@ def run_nmr_calc( ) +# L audit fix (2026-07-14): bump this whenever the patch bodies below +# change — it doubles as the idempotency sentinel's value, so a bumped +# version forces re-patching instead of silently keeping stale closures +# from an older QuantUI version installed earlier in the process. +_NMR_COMPAT_PATCH_VERSION = 1 +_NMR_PATCH_VERSION_ATTR = "_quantui_nmr_compat_patch_version" + + +def _ensure_nmr_compat_patches_applied() -> None: + """Idempotently patch pyscf.prop.nmr for QuantUI-specific compatibility fixes. + + Both patches below used to be applied unconditionally on every NMR + calculation, re-defining the same closures and reassigning the same + module attributes on every call even though nothing about the + installed pyscf/pyscf-properties changes between calls. Each patch + is now a one-time, idempotent operation per process: a sentinel + attribute (versioned via ``_NMR_COMPAT_PATCH_VERSION``) on the + currently-installed function is checked first, so repeated NMR runs + are no-ops here. + """ + # pyscf-properties 0.1.0 gen_vind hardcodes reshape(3, nmo, nocc). + # pyscf 2.x krylov reduces the batch below 3 via linear-dependency masking, + # causing "cannot reshape array of size N into shape (3,nmo,nocc)". + # Patch gen_vind to use reshape(-1, nmo, nocc) so any batch size works. + try: + from functools import reduce as _reduce_nmr + + import numpy as _np + import pyscf.prop.nmr.rhf as _prop_nmr_rhf + from pyscf import lib as _pyscf_lib_nmr + + if ( + getattr(_prop_nmr_rhf.gen_vind, _NMR_PATCH_VERSION_ATTR, 0) + < _NMR_COMPAT_PATCH_VERSION + ): + + def _fixed_gen_vind(mf_arg, mo_coeff, mo_occ): + vresp = mf_arg.gen_response(singlet=True, hermi=2) + occidx = mo_occ > 0 + orbo = mo_coeff[:, occidx] + nocc = orbo.shape[1] + _nao, nmo = mo_coeff.shape + + def vind(mo1): + _mo1 = _np.asarray(mo1).reshape(-1, nmo, nocc) + dm1 = _np.asarray( + [ + _reduce_nmr(_np.dot, (mo_coeff, x * 2, orbo.T.conj())) + for x in _mo1 + ] + ) + dm1 = dm1 - dm1.transpose(0, 2, 1).conj() + v1mo = _pyscf_lib_nmr.einsum( + "xpq,pi,qj->xij", vresp(dm1), mo_coeff.conj(), orbo + ) + return v1mo.ravel() + + return vind + + setattr(_fixed_gen_vind, _NMR_PATCH_VERSION_ATTR, _NMR_COMPAT_PATCH_VERSION) + _prop_nmr_rhf.gen_vind = _fixed_gen_vind + except (ImportError, AttributeError) as exc: # noqa: BLE001 — optional probe + logger.debug("pyscf.prop.nmr.rhf.gen_vind patch not applied: %s", exc) + + # pyscf-properties 0.1.0 get_vxc_giao computes + # blksize = min(int(X*BLKSIZE)*BLKSIZE, ngrids) + # which equals ngrids when ngrids < X*BLKSIZE, and ngrids may not be + # divisible by BLKSIZE. pyscf 2.x block_loop asserts blksize%BLKSIZE==0. + # Patch get_vxc_giao to round blksize down to the nearest BLKSIZE multiple. + try: + import numpy as _np_rks + import pyscf.prop.nmr.rks as _prop_nmr_rks + from pyscf.dft import numint as _numint_rks + + if ( + getattr(_prop_nmr_rks.get_vxc_giao, _NMR_PATCH_VERSION_ATTR, 0) + < _NMR_COMPAT_PATCH_VERSION + ): + + def _fixed_get_vxc_giao( + ni, mol, grids, xc_code, dms, max_memory=2000, verbose=None + ): + xctype = ni._xc_type(xc_code) + make_rho, nset, nao = ni._gen_rho_evaluator(mol, dms, hermi=1) + ngrids = len(grids.weights) + _BLKSIZE = _numint_rks.BLKSIZE + _raw_blk = int(max_memory / 12 * 1e6 / 8 / nao / _BLKSIZE) * _BLKSIZE + blksize = max(_BLKSIZE, (min(_raw_blk, ngrids) // _BLKSIZE) * _BLKSIZE) + shls_slice = (0, mol.nbas) + ao_loc = mol.ao_loc_nr() + + vmat = _np_rks.zeros((3, nao, nao)) + if xctype == "LDA": + buf = _np_rks.empty((4, blksize, nao)) + ao_deriv = 0 + for ao, mask, weight, coords in ni.block_loop( + mol, grids, nao, ao_deriv, max_memory, blksize=blksize, buf=buf + ): + rho = make_rho(0, ao, mask, "LDA") + vxc = ni.eval_xc(xc_code, rho, 0, deriv=1)[1] + vrho = vxc[0] + aow = _np_rks.einsum("pi,p->pi", ao, weight * vrho) + giao = mol.eval_gto( + "GTOval_ig", coords, comp=3, non0tab=mask, out=buf[1:] + ) + vmat[0] += _numint_rks._dot_ao_ao( + mol, aow, giao[0], mask, shls_slice, ao_loc + ) + vmat[1] += _numint_rks._dot_ao_ao( + mol, aow, giao[1], mask, shls_slice, ao_loc + ) + vmat[2] += _numint_rks._dot_ao_ao( + mol, aow, giao[2], mask, shls_slice, ao_loc + ) + rho = vxc = vrho = aow = None + elif xctype == "GGA": + buf = _np_rks.empty((10, blksize, nao)) + ao_deriv = 1 + for ao, mask, weight, coords in ni.block_loop( + mol, grids, nao, ao_deriv, max_memory, blksize=blksize, buf=buf + ): + rho = make_rho(0, ao, mask, "GGA") + vxc = ni.eval_xc(xc_code, rho, 0, deriv=1)[1] + vrho, vsigma = vxc[:2] + wv = _np_rks.empty_like(rho) + wv[0] = weight * vrho + wv[1:] = rho[1:] * (weight * vsigma * 2) + aow = _np_rks.einsum("npi,np->pi", ao[:4], wv) + giao = mol.eval_gto( + "GTOval_ig", coords, 3, non0tab=mask, out=buf[4:] + ) + vmat[0] += _numint_rks._dot_ao_ao( + mol, aow, giao[0], mask, shls_slice, ao_loc + ) + vmat[1] += _numint_rks._dot_ao_ao( + mol, aow, giao[1], mask, shls_slice, ao_loc + ) + vmat[2] += _numint_rks._dot_ao_ao( + mol, aow, giao[2], mask, shls_slice, ao_loc + ) + giao = mol.eval_gto( + "GTOval_ipig", coords, 9, non0tab=mask, out=buf[1:] + ) + _prop_nmr_rks._gga_sum_( + vmat, mol, ao, giao, wv, mask, shls_slice, ao_loc + ) + rho = vxc = vrho = vsigma = wv = aow = None + elif xctype == "MGGA": + raise NotImplementedError("meta-GGA") + + return vmat - vmat.transpose(0, 2, 1) + + setattr( + _fixed_get_vxc_giao, _NMR_PATCH_VERSION_ATTR, _NMR_COMPAT_PATCH_VERSION + ) + _prop_nmr_rks.get_vxc_giao = _fixed_get_vxc_giao + except (ImportError, AttributeError) as exc: # noqa: BLE001 — optional probe + logger.debug("pyscf.prop.nmr.rks.get_vxc_giao patch not applied: %s", exc) + + def _run_nmr_calc_body( *, molecule: Molecule, @@ -171,129 +381,7 @@ def _run_nmr_calc_body( "Install pyscf-properties: pip install pyscf-properties" ) from exc - # pyscf-properties 0.1.0 gen_vind hardcodes reshape(3, nmo, nocc). - # pyscf 2.x krylov reduces the batch below 3 via linear-dependency masking, - # causing "cannot reshape array of size N into shape (3,nmo,nocc)". - # Patch gen_vind to use reshape(-1, nmo, nocc) so any batch size works. - try: - from functools import reduce as _reduce_nmr - - import pyscf.prop.nmr.rhf as _prop_nmr_rhf - from pyscf import lib as _pyscf_lib_nmr - - def _fixed_gen_vind(mf_arg, mo_coeff, mo_occ): - vresp = mf_arg.gen_response(singlet=True, hermi=2) - occidx = mo_occ > 0 - orbo = mo_coeff[:, occidx] - nocc = orbo.shape[1] - _nao, nmo = mo_coeff.shape - - def vind(mo1): - _mo1 = _np.asarray(mo1).reshape(-1, nmo, nocc) - dm1 = _np.asarray( - [ - _reduce_nmr(_np.dot, (mo_coeff, x * 2, orbo.T.conj())) - for x in _mo1 - ] - ) - dm1 = dm1 - dm1.transpose(0, 2, 1).conj() - v1mo = _pyscf_lib_nmr.einsum( - "xpq,pi,qj->xij", vresp(dm1), mo_coeff.conj(), orbo - ) - return v1mo.ravel() - - return vind - - _prop_nmr_rhf.gen_vind = _fixed_gen_vind - except (ImportError, AttributeError) as exc: # noqa: BLE001 — optional probe - logger.debug("pyscf.prop.nmr.rhf.gen_vind patch not applied: %s", exc) - - # pyscf-properties 0.1.0 get_vxc_giao computes - # blksize = min(int(X*BLKSIZE)*BLKSIZE, ngrids) - # which equals ngrids when ngrids < X*BLKSIZE, and ngrids may not be - # divisible by BLKSIZE. pyscf 2.x block_loop asserts blksize%BLKSIZE==0. - # Patch get_vxc_giao to round blksize down to the nearest BLKSIZE multiple. - try: - import numpy as _np_rks - import pyscf.prop.nmr.rks as _prop_nmr_rks - from pyscf.dft import numint as _numint_rks - - def _fixed_get_vxc_giao( - ni, mol, grids, xc_code, dms, max_memory=2000, verbose=None - ): - xctype = ni._xc_type(xc_code) - make_rho, nset, nao = ni._gen_rho_evaluator(mol, dms, hermi=1) - ngrids = len(grids.weights) - _BLKSIZE = _numint_rks.BLKSIZE - _raw_blk = int(max_memory / 12 * 1e6 / 8 / nao / _BLKSIZE) * _BLKSIZE - blksize = max(_BLKSIZE, (min(_raw_blk, ngrids) // _BLKSIZE) * _BLKSIZE) - shls_slice = (0, mol.nbas) - ao_loc = mol.ao_loc_nr() - - vmat = _np_rks.zeros((3, nao, nao)) - if xctype == "LDA": - buf = _np_rks.empty((4, blksize, nao)) - ao_deriv = 0 - for ao, mask, weight, coords in ni.block_loop( - mol, grids, nao, ao_deriv, max_memory, blksize=blksize, buf=buf - ): - rho = make_rho(0, ao, mask, "LDA") - vxc = ni.eval_xc(xc_code, rho, 0, deriv=1)[1] - vrho = vxc[0] - aow = _np_rks.einsum("pi,p->pi", ao, weight * vrho) - giao = mol.eval_gto( - "GTOval_ig", coords, comp=3, non0tab=mask, out=buf[1:] - ) - vmat[0] += _numint_rks._dot_ao_ao( - mol, aow, giao[0], mask, shls_slice, ao_loc - ) - vmat[1] += _numint_rks._dot_ao_ao( - mol, aow, giao[1], mask, shls_slice, ao_loc - ) - vmat[2] += _numint_rks._dot_ao_ao( - mol, aow, giao[2], mask, shls_slice, ao_loc - ) - rho = vxc = vrho = aow = None - elif xctype == "GGA": - buf = _np_rks.empty((10, blksize, nao)) - ao_deriv = 1 - for ao, mask, weight, coords in ni.block_loop( - mol, grids, nao, ao_deriv, max_memory, blksize=blksize, buf=buf - ): - rho = make_rho(0, ao, mask, "GGA") - vxc = ni.eval_xc(xc_code, rho, 0, deriv=1)[1] - vrho, vsigma = vxc[:2] - wv = _np_rks.empty_like(rho) - wv[0] = weight * vrho - wv[1:] = rho[1:] * (weight * vsigma * 2) - aow = _np_rks.einsum("npi,np->pi", ao[:4], wv) - giao = mol.eval_gto( - "GTOval_ig", coords, 3, non0tab=mask, out=buf[4:] - ) - vmat[0] += _numint_rks._dot_ao_ao( - mol, aow, giao[0], mask, shls_slice, ao_loc - ) - vmat[1] += _numint_rks._dot_ao_ao( - mol, aow, giao[1], mask, shls_slice, ao_loc - ) - vmat[2] += _numint_rks._dot_ao_ao( - mol, aow, giao[2], mask, shls_slice, ao_loc - ) - giao = mol.eval_gto( - "GTOval_ipig", coords, 9, non0tab=mask, out=buf[1:] - ) - _prop_nmr_rks._gga_sum_( - vmat, mol, ao, giao, wv, mask, shls_slice, ao_loc - ) - rho = vxc = vrho = vsigma = wv = aow = None - elif xctype == "MGGA": - raise NotImplementedError("meta-GGA") - - return vmat - vmat.transpose(0, 2, 1) - - _prop_nmr_rks.get_vxc_giao = _fixed_get_vxc_giao - except (ImportError, AttributeError) as exc: # noqa: BLE001 — optional probe - logger.debug("pyscf.prop.nmr.rks.get_vxc_giao patch not applied: %s", exc) + _ensure_nmr_compat_patches_applied() try: if method_upper == "RHF": @@ -316,8 +404,7 @@ def _fixed_get_vxc_giao( else: shielding_iso.append(float(arr)) - key = f"{method}/{basis}" - ref_map = _config.NMR_REFERENCE_SHIELDINGS.get(key, _config.NMR_DEFAULT_REFERENCE) + ref_map, matched_ref_key, is_fallback_ref = resolve_nmr_reference(method, basis) ref_H = float(ref_map.get("H", _config.NMR_DEFAULT_REFERENCE["H"])) ref_C = float(ref_map.get("C", _config.NMR_DEFAULT_REFERENCE["C"])) @@ -337,4 +424,6 @@ def _fixed_get_vxc_giao( basis=basis, formula=molecule.get_formula(), converged=converged, + reference_key=matched_ref_key, + is_fallback_reference=is_fallback_ref, ) diff --git a/quantui/optimizer.py b/quantui/optimizer.py index 3a69924..3e048fe 100644 --- a/quantui/optimizer.py +++ b/quantui/optimizer.py @@ -50,14 +50,12 @@ from typing import IO, Any, List, Optional from .ase_bridge import ASE_AVAILABLE, atoms_to_molecule, molecule_to_atoms +from .config import BOHR_TO_ANGSTROM as _BOHR_TO_ANG from .molecule import Molecule from .session_calc import HARTREE_TO_EV logger = logging.getLogger(__name__) -# Unit conversion: 1 Bohr = 0.529177249 Angstrom (NIST 2018 CODATA) -_BOHR_TO_ANG: float = 0.529177249 - # Defaults also exposed in config.py for the notebook UI DEFAULT_FMAX: float = 0.05 # eV/Å — tight enough for educational use DEFAULT_OPT_STEPS: int = 200 # generous upper limit for small molecules @@ -349,6 +347,23 @@ def optimize_geometry( " # or: conda install -c conda-forge ase" ) + # Post-HF methods (MP2/CCSD/CCSD(T)) have no special-casing in + # _QuantUIPySCFCalc — without this guard, method='CCSD' silently falls + # into the DFT branch (sets mf.xc = "CCSD") and fails deep inside PySCF + # with a cryptic "LibXCFunctional: name 'CCSD' not found" instead of a + # clear message. No analytical post-HF nuclear gradients are wired up + # here, so these methods are single-point only (see session_calc.py). + from . import config as _config + + if method.strip().upper() in _config.POST_HF_METHODS: + raise ValueError( + f"'{method}' is a post-HF method and cannot be used for geometry " + "optimization — QuantUI only has analytical gradients wired up " + "for HF/DFT methods here. Optimize with RHF, UHF, or a DFT " + f"functional, then run a Single Point calculation with '{method}' " + "on the optimized geometry." + ) + try: import pyscf as _pyscf # noqa: F401 — presence check except ImportError as exc: diff --git a/quantui/orbital_visualization.py b/quantui/orbital_visualization.py index c247147..132cb3d 100644 --- a/quantui/orbital_visualization.py +++ b/quantui/orbital_visualization.py @@ -531,6 +531,47 @@ def orbital_summary_html(info: OrbitalInfo) -> str: # ============================================================================ +def infer_charge_and_spin(mol_atom: list, mo_occ: np.ndarray | list) -> Tuple[int, int]: + """Infer ``(charge, spin)`` for a ``gto.Mole`` from atoms + MO occupations. + + Cube/isosurface generation from saved MO data does not have direct access + to the original ``Molecule.charge`` / ``Molecule.multiplicity`` — only the + atom list and the MO coefficients/occupations that came out of the SCF. + PySCF's ``Mole.build()`` requires ``charge``/``spin`` consistent with the + actual electron count, so passing the default ``charge=0, spin=0`` fails + to build for any charged or open-shell (odd-electron) molecule. + + This reconstructs both from data that's always available: + + - ``spin`` (PySCF's ``2S = n_alpha - n_beta``) is 0 when ``mo_occ`` is + 1-D (closed-shell RHF/RKS — including the MP2/CCSD/CCSD(T) paths, which + always run on an RHF reference), or ``n_alpha - n_beta`` when ``mo_occ`` + is 2-D (UHF/UKS, shape ``(2, n_mo)``). + - ``charge`` is the nuclear charge (sum of atomic numbers in ``mol_atom``) + minus the total electron count (``sum(mo_occ)`` over all spin channels). + + Returns ``(0, 0)`` if ``mol_atom`` or ``mo_occ`` is falsy/``None`` so + callers can pass through directly without a separate None-check. + """ + if not mol_atom or mo_occ is None: + return 0, 0 + from .molecule import ATOMIC_NUMBERS + + occ = np.asarray(mo_occ, dtype=float) + if occ.ndim == 2: + n_alpha = float(occ[0].sum()) + n_beta = float(occ[1].sum()) + spin = int(round(n_alpha - n_beta)) + n_electrons = n_alpha + n_beta + else: + spin = 0 + n_electrons = float(occ.sum()) + + nuclear_charge = sum(ATOMIC_NUMBERS.get(sym, 0) for sym, _ in mol_atom) + charge = int(round(nuclear_charge - n_electrons)) + return charge, spin + + def generate_cube_file( results_path: Path, orbital_index: int, @@ -582,6 +623,7 @@ def generate_cube_file( data = np.load(results_path, allow_pickle=True) mo_coeff = data["mo_coeff"] + mo_occ = data["mo_occ"] if "mo_occ" in data else None if mo_coeff.ndim == 3: mo_coeff = mo_coeff[0] @@ -594,7 +636,22 @@ def generate_cube_file( "Re-run the calculation with the updated script template." ) - mol = gto.M(atom=atom_str, basis=basis_str, unit="Angstrom") + # Charge/spin aren't stored in results.npz — infer them from the MO + # occupations (when present) so charged/open-shell molecules don't fail + # to build. mol_atom here is a PySCF-format string, not the (symbol, + # coords) tuple list infer_charge_and_spin expects, so parse it first. + charge, spin = 0, 0 + if mo_occ is not None: + parsed_atoms = [ + (tok.split()[0], [0.0, 0.0, 0.0]) + for tok in atom_str.replace(";", "\n").splitlines() + if tok.strip() + ] + charge, spin = infer_charge_and_spin(parsed_atoms, mo_occ) + + mol = gto.M( + atom=atom_str, basis=basis_str, unit="Angstrom", charge=charge, spin=spin + ) output_path = Path(output_path) output_path.parent.mkdir(parents=True, exist_ok=True) @@ -623,6 +680,7 @@ def generate_cube_from_arrays( ny: int = 60, nz: int = 60, margin: float = 5.0, + charge: int = 0, spin: int = 0, ) -> Path: """ @@ -650,6 +708,13 @@ def generate_cube_from_arrays( Grid resolution along each axis. margin : float Extra space (Bohr) beyond atomic extents. + charge : int + Total molecular charge. Required for charged species (e.g. H3O+, + NH4+, OH-) — without it PySCF's electron-count check fails at + ``mol.build()``. Default 0 (neutral). + spin : int + PySCF's ``2S = n_alpha - n_beta``. Required for open-shell + (odd-electron) molecules — default 0 assumes closed-shell. Returns ------- @@ -670,7 +735,9 @@ def generate_cube_from_arrays( " conda install -c conda-forge pyscf" ) from exc - mol = gto.M(atom=mol_atom, basis=mol_basis, unit="Angstrom", spin=spin) + mol = gto.M( + atom=mol_atom, basis=mol_basis, unit="Angstrom", charge=charge, spin=spin + ) coeff = np.asarray(mo_coeff) if coeff.ndim == 3: diff --git a/quantui/pes_scan.py b/quantui/pes_scan.py index 4bf7393..d6db109 100644 --- a/quantui/pes_scan.py +++ b/quantui/pes_scan.py @@ -76,10 +76,24 @@ class PESScanResult: # ── Convenience properties ────────────────────────────────────────────── + def _finite_energies(self) -> List[float]: + """``energies_hartree`` with failed-point NaN placeholders dropped. + + M6 audit fix (2026-07-14): a failed scan point appends + ``float("nan")`` to ``energies_hartree`` (see :func:`run_pes_scan`). + Python's ``min``/``max`` are order-dependent with NaN present — a + NaN as the first element "wins" (everything compares False against + it) and poisons the result; a NaN later in the list is correctly + ignored. Filtering NaN out before any min/max call makes the result + deterministic regardless of *which* scan point failed. + """ + return [e for e in self.energies_hartree if math.isfinite(e)] + @property def energy_hartree(self) -> float: - """Minimum SCF energy across all scan points (Hartrees).""" - return min(self.energies_hartree) if self.energies_hartree else float("nan") + """Minimum SCF energy across all *successful* scan points (Hartrees).""" + finite = self._finite_energies() + return min(finite) if finite else float("nan") @property def energy_ev(self) -> float: @@ -98,10 +112,15 @@ def n_steps(self) -> int: @property def energies_relative_kcal(self) -> List[float]: - """Energy relative to the lowest scan point, in kcal/mol.""" - if not self.energies_hartree: + """Energy relative to the lowest successful scan point, in kcal/mol. + + Failed points (NaN in ``energies_hartree``) stay NaN here too — + callers plotting this list should skip non-finite entries. + """ + finite = self._finite_energies() + if not finite: return [] - e_min = min(self.energies_hartree) + e_min = min(finite) return [(e - e_min) * _HARTREE_TO_KCAL for e in self.energies_hartree] @property @@ -121,10 +140,11 @@ def scan_coordinate_label(self) -> str: def summary(self) -> str: """Return a multi-line human-readable result summary.""" - if not self.energies_hartree: + finite = self._finite_energies() + if not finite: return "No scan points computed." - e_min = min(self.energies_hartree) - e_max = max(self.energies_hartree) + e_min = min(finite) + e_max = max(finite) barrier = (e_max - e_min) * _HARTREE_TO_KCAL min_idx = self.energies_hartree.index(e_min) lines = [ @@ -201,6 +221,22 @@ def run_pes_scan( "ASE is not installed — cannot run PES scan.\n" " pip install 'ase>=3.22.0'" ) + + # Post-HF methods (MP2/CCSD/CCSD(T)) have no special-casing in + # _QuantUIPySCFCalc (shared with optimizer.py) — without this guard, + # method='CCSD' silently falls into the DFT branch (sets mf.xc = + # "CCSD") and fails deep inside PySCF with a cryptic "LibXCFunctional: + # name 'CCSD' not found" instead of a clear message. + from . import config as _config + + if method.strip().upper() in _config.POST_HF_METHODS: + raise ValueError( + f"'{method}' is a post-HF method and cannot be used for a PES " + "scan — QuantUI only has analytical gradients wired up for " + "HF/DFT methods here. Scan with RHF, UHF, or a DFT functional " + "instead." + ) + try: import pyscf as _pyscf # noqa: F401 except ImportError as exc: @@ -261,6 +297,13 @@ def run_pes_scan( energies_hartree: List[float] = [] coordinates_list: List[Molecule] = [] converged_all = True + # M6 audit fix (2026-07-14): on a failed scan point, fall back to the + # last successfully-computed geometry rather than the original input + # molecule. Snapping every failed frame back to the starting geometry + # produced a bogus discontinuity in the trajectory animation/plot — + # the last-good geometry is a far more sensible placeholder for "we + # don't know where this point landed, but it wasn't back at the start." + _last_good_molecule = molecule i1, i2 = atom_indices[0], atom_indices[1] i3 = atom_indices[2] if len(atom_indices) >= 3 else 0 @@ -285,17 +328,25 @@ def run_pes_scan( ok = True else: if scan_type == "bond": - constraint = FixInternals(bonds=[(val, [i1, i2])]) + constraint = FixInternals(bonds=[[val, [i1, i2]]]) elif scan_type == "angle": atoms.set_angle(i1, i2, i3, val) - constraint = FixInternals( - angles=[(math.radians(val), [i1, i2, i3])] - ) + # M7 audit fix (2026-07-14): ASE's radian-based `angles=` + # kwarg is not just deprecated, it's flat-out broken with + # the currently-targeted ASE (>=3.22, verified against + # 3.29.0) — internally it does + # ``np.asarray(angles); angles[:, 0] = ...`` to convert + # to degrees, which raises "setting an array element + # with a sequence" for any real angle constraint (the + # per-entry [value, [3 indices]] shape isn't + # rectangular). Every angle/dihedral PES scan silently + # failed at 100% of its points as a result. `angles_deg` + # takes the value directly in degrees and skips that + # broken reshape entirely. + constraint = FixInternals(angles_deg=[[val, [i1, i2, i3]]]) else: # dihedral atoms.set_dihedral(i1, i2, i3, i4, val) - constraint = FixInternals( - dihedrals=[(math.radians(val), [i1, i2, i3, i4])] - ) + constraint = FixInternals(dihedrals_deg=[[val, [i1, i2, i3, i4]]]) atoms.set_constraint(constraint) @@ -321,6 +372,7 @@ def run_pes_scan( atoms, charge=molecule.charge, multiplicity=molecule.multiplicity ) coordinates_list.append(mol_at_point) + _last_good_molecule = mol_at_point _stream.write( f" E = {e_ha:.8f} Ha ({'converged' if ok else 'not converged'})\n" @@ -329,7 +381,7 @@ def run_pes_scan( except Exception as exc: _stream.write(f" ⚠ Scan point {step_num} failed: {exc}\n") energies_hartree.append(float("nan")) - coordinates_list.append(molecule) + coordinates_list.append(_last_good_molecule) converged_all = False finally: diff --git a/quantui/progress.py b/quantui/progress.py index 4abba1b..97f2f99 100644 --- a/quantui/progress.py +++ b/quantui/progress.py @@ -24,6 +24,7 @@ from __future__ import annotations +import html from typing import List, Optional import ipywidgets as widgets @@ -96,10 +97,10 @@ def _render(self) -> None: line = ( f'
' - f"{icon} Step {i + 1}: {label}" + f"{icon} Step {i + 1}: {html.escape(label)}" ) if self._messages[i]: - line += f" — {self._messages[i]}" + line += f" — {html.escape(self._messages[i])}" line += "
" lines.append(line) diff --git a/quantui/pubchem.py b/quantui/pubchem.py index b7e5e33..4db54b1 100644 --- a/quantui/pubchem.py +++ b/quantui/pubchem.py @@ -66,14 +66,21 @@ def _http_get( to :class:`PubChemAPIError` exactly as before. """ timeout = timeout if timeout is not None else config.PUBCHEM_TIMEOUT_S + # M10 audit fix (2026-07-14): if config.PUBCHEM_MAX_RETRIES were ever 0 + # (or negative), `range(config.PUBCHEM_MAX_RETRIES)` would iterate zero + # times, leaving `response` at its None initializer and returning None + # from a function typed to return requests.Response — every caller + # then hits AttributeError on response.status_code. Always attempt at + # least once regardless of the configured retry count. + max_attempts = max(1, config.PUBCHEM_MAX_RETRIES) response = None - for attempt in range(config.PUBCHEM_MAX_RETRIES): + for attempt in range(max_attempts): _throttle() response = requests.get(url, params=params, timeout=timeout) if response.status_code != 503: return response # Throttled — back off (capped) and retry, unless this was the last try. - if attempt < config.PUBCHEM_MAX_RETRIES - 1: + if attempt < max_attempts - 1: backoff = min( config.PUBCHEM_BACKOFF_BASE_S * (2**attempt), config.PUBCHEM_BACKOFF_MAX_S, @@ -82,7 +89,7 @@ def _http_get( "PubChem throttled (503); retrying in %.1fs (attempt %d/%d)", backoff, attempt + 1, - config.PUBCHEM_MAX_RETRIES, + max_attempts, ) time.sleep(backoff) # Exhausted retries — hand the last 503 back; caller raises via raise_for_status. @@ -115,7 +122,7 @@ def search_molecule_by_name(name: str) -> int: name: Common name or IUPAC name of the molecule Returns: - int: PubChem Compound ID (CID) if found, None otherwise + int: PubChem Compound ID (CID) Raises: PubChemAPIError: If API request fails @@ -143,7 +150,7 @@ def search_molecule_by_name(name: str) -> int: except requests.RequestException as e: logger.error(f"PubChem API request failed: {e}") - raise PubChemAPIError(f"Failed to connect to PubChem: {e}") + raise PubChemAPIError(f"Failed to connect to PubChem: {e}") from e def search_cid_by_inchikey(inchikey: str) -> int: @@ -164,7 +171,7 @@ def search_cid_by_inchikey(inchikey: str) -> int: return int(cids[0]) except requests.RequestException as e: logger.error(f"PubChem InChIKey request failed: {e}") - raise PubChemAPIError(f"Failed to connect to PubChem: {e}") + raise PubChemAPIError(f"Failed to connect to PubChem: {e}") from e def search_cids_by_name(name: str) -> list: @@ -184,7 +191,7 @@ def search_cids_by_name(name: str) -> list: ] except requests.RequestException as e: logger.error(f"PubChem CID-list request failed: {e}") - raise PubChemAPIError(f"Failed to connect to PubChem: {e}") + raise PubChemAPIError(f"Failed to connect to PubChem: {e}") from e def search_pubchem_candidates(query: str, max_results: int = 10) -> list: @@ -208,7 +215,7 @@ def search_pubchem_candidates(query: str, max_results: int = 10) -> list: props = response.json().get("PropertyTable", {}).get("Properties", []) except requests.RequestException as e: logger.error(f"PubChem property request failed: {e}") - raise PubChemAPIError(f"Failed to connect to PubChem: {e}") + raise PubChemAPIError(f"Failed to connect to PubChem: {e}") from e # Preserve the CID search order (the property endpoint may reorder). by_cid = {int(p.get("CID")): p for p in props if p.get("CID") is not None} @@ -274,7 +281,7 @@ def get_molecule_sdf(cid: int, conformer_3d: bool = True) -> str: except requests.RequestException as e: logger.error(f"PubChem SDF request failed: {e}") - raise PubChemAPIError(f"Failed to retrieve molecule: {e}") + raise PubChemAPIError(f"Failed to retrieve molecule: {e}") from e def _separate_fragments(mol: Any, min_gap: float = 3.0) -> None: @@ -403,7 +410,7 @@ def sdf_to_xyz(sdf_content: str) -> Tuple[str, Dict[str, Any]]: except Exception as e: logger.error(f"SDF to XYZ conversion failed: {e}") - raise ValueError(f"Failed to convert SDF to XYZ: {e}") + raise ValueError(f"Failed to convert SDF to XYZ: {e}") from e def fetch_molecule( @@ -542,9 +549,21 @@ def check_pubchem_availability() -> bool: Returns: bool: True if PubChem is accessible, False otherwise """ + # M10 audit fix (2026-07-14): this used to call requests.get() directly, + # bypassing the shared client-side rate limiter (_throttle()) that every + # other function in this module goes through — a burst of concurrent + # availability checks (e.g. several students in a classroom clicking + # "check connection" around the same time) could exceed PubChem's + # server-side throttle. It also hardcoded timeout=5 instead of the + # config.PUBCHEM_AVAILABILITY_TIMEOUT_S constant defined specifically + # for this probe, so changing that constant silently had no effect + # here. Calls _throttle() directly (not the full _http_get retry loop — + # this is meant to be a quick, no-retry reachability probe) and uses + # the configured timeout. try: url = f"{PUBCHEM_BASE_URL}/compound/cid/962/property/MolecularFormula/JSON" - response = requests.get(url, timeout=5) + _throttle() + response = requests.get(url, timeout=config.PUBCHEM_AVAILABILITY_TIMEOUT_S) return bool(response.status_code == 200) except Exception: return False @@ -635,7 +654,7 @@ def smiles_to_xyz(smiles: str, optimize_3d: bool = True) -> Tuple[str, Dict[str, except Exception as e: logger.error(f"SMILES to XYZ conversion failed: {e}") - raise ValueError(f"Failed to convert SMILES to XYZ: {e}") + raise ValueError(f"Failed to convert SMILES to XYZ: {e}") from e def inchi_to_xyz(inchi: str, optimize_3d: bool = True) -> Tuple[str, Dict[str, Any]]: @@ -691,7 +710,7 @@ def inchi_to_xyz(inchi: str, optimize_3d: bool = True) -> Tuple[str, Dict[str, A except Exception as e: logger.error(f"InChI to XYZ conversion failed: {e}") - raise ValueError(f"Failed to convert InChI to XYZ: {e}") + raise ValueError(f"Failed to convert InChI to XYZ: {e}") from e def student_friendly_smiles_to_xyz(smiles: str) -> Tuple[Optional[str], str]: @@ -805,21 +824,38 @@ def generate_2d_structure_svg( # Build mol from XYZ from rdkit.Chem import rdDetermineBonds - rdkit_mol = Chem.Mol() - conf = Chem.Conformer() - - for i, line in enumerate(lines[2:]): # Skip first 2 lines + # Collect valid atom lines first so the conformer can be sized + # up front and the atom index used for SetAtomPosition always + # matches the just-added atom (skipped malformed lines used to + # desync the two). + atom_lines = [] + for line in lines[2:]: # Skip first 2 lines (count + comment) parts = line.split() if len(parts) < 4: continue + atom_lines.append(parts) + + if not atom_lines: + raise ValueError("No atom lines found in XYZ string") + + # Chem.Mol() is immutable — AddAtom only exists on RWMol. + # Build on RWMol, then convert to an immutable Mol via + # GetMol() before DetermineBonds (mirrors the working + # Chem.MolFromXYZBlock() + DetermineBonds() pattern used + # elsewhere in QuantUI, e.g. preopt.py). + rw_mol = Chem.RWMol() + conf = Chem.Conformer(len(atom_lines)) + + for i, parts in enumerate(atom_lines): symbol = parts[0] x, y, z = float(parts[1]), float(parts[2]), float(parts[3]) atom = Chem.Atom(symbol) - rdkit_mol.AddAtom(atom) # type: ignore[attr-defined] + rw_mol.AddAtom(atom) conf.SetAtomPosition(i, (x, y, z)) - rdkit_mol.AddConformer(conf) # type: ignore[attr-defined] + rw_mol.AddConformer(conf) + rdkit_mol = rw_mol.GetMol() # Determine bonds rdDetermineBonds.DetermineBonds(rdkit_mol) diff --git a/quantui/results_storage.py b/quantui/results_storage.py index 289274b..48eb819 100644 --- a/quantui/results_storage.py +++ b/quantui/results_storage.py @@ -28,13 +28,25 @@ import re from datetime import datetime from pathlib import Path -from typing import TYPE_CHECKING, Optional +from typing import TYPE_CHECKING, Any, Optional + +from .config import BOHR_TO_ANGSTROM as _BOHR_TO_ANGSTROM if TYPE_CHECKING: pass # result types accepted via duck typing; no hard import needed _SCHEMA_VERSION = 2 +# Molden's [FR-COORD] block is defined (theochem.ru.nl/molden/molden_format.html) +# to always be in Bohr, regardless of the unit tag on [Atoms] — a Molden-format +# quirk. pyscf_mol_atom (the source for both blocks) is Angstrom throughout +# QuantUI, so [FR-COORD] needs an explicit conversion; [Atoms]/[GTO]/[MO] (via +# molden.from_mo / molden.header, built from a mol with implicit unit="Angstrom") +# do not. Derived from config.BOHR_TO_ANGSTROM (pyscf.data.nist.BOHR) rather +# than a separately hand-typed literal, so this stays consistent with the +# other Bohr<->Angstrom conversions in the codebase. +_ANGSTROM_TO_BOHR = 1.0 / _BOHR_TO_ANGSTROM + def _default_results_dir() -> Path: env = os.environ.get("QUANTUI_RESULTS_DIR") @@ -56,6 +68,22 @@ def _opt_float(x: object) -> Optional[float]: return None +def _opt_int(x: object) -> Optional[int]: + """Coerce an optional (possibly numpy) scalar to a JSON-safe int or None. + + ``json.dumps`` accepts ``numpy.float64``/``numpy.float32`` transparently + (they subclass ``float``), but ``numpy.int64``/``numpy.bool_`` do not + subclass ``int``/``bool`` and raise ``TypeError`` unconverted — this + normalizes any duck-typed result's numpy scalar to a plain ``int``. + """ + if x is None: + return None + try: + return int(x) # type: ignore[arg-type] + except (TypeError, ValueError): + return None + + def _opt_float_list(x: object) -> Optional[list]: """Coerce an optional iterable of numbers to a JSON-safe list of floats.""" if x is None: @@ -143,10 +171,17 @@ def save_result( _collision += 1 dest.mkdir(parents=True) - _e_ha = getattr(result, "energy_hartree", float("nan")) + _e_ha_raw = getattr(result, "energy_hartree", float("nan")) + _e_ha = _opt_float(_e_ha_raw) + if _e_ha is None: + _e_ha = float("nan") # energy_ev may be a property (SessionResult) or absent (OptimizationResult # and new types also define it as a property, so getattr works for all). - _e_ev = getattr(result, "energy_ev", _e_ha * _HARTREE_TO_EV) + _e_ev = _opt_float(getattr(result, "energy_ev", _e_ha * _HARTREE_TO_EV)) + if _e_ev is None: + _e_ev = _e_ha * _HARTREE_TO_EV + + _converged = getattr(result, "converged", None) data: dict = { "_schema_version": _SCHEMA_VERSION, @@ -157,14 +192,20 @@ def save_result( "basis": getattr(result, "basis", "?"), "energy_hartree": _e_ha, "energy_ev": _e_ev, - "homo_lumo_gap_ev": getattr(result, "homo_lumo_gap_ev", None), - "converged": getattr(result, "converged", None), - "n_iterations": getattr(result, "n_iterations", -1), + "homo_lumo_gap_ev": _opt_float(getattr(result, "homo_lumo_gap_ev", None)), + "converged": None if _converged is None else bool(_converged), + "n_iterations": _opt_int(getattr(result, "n_iterations", -1)), # Post-HF correlation breakdown — None for HF/DFT. Persisted so the # saved-result card can show the HF reference + correlation rows. - "mp2_correlation_hartree": getattr(result, "mp2_correlation_hartree", None), - "ccsd_correlation_hartree": getattr(result, "ccsd_correlation_hartree", None), - "ccsd_t_correction_hartree": getattr(result, "ccsd_t_correction_hartree", None), + "mp2_correlation_hartree": _opt_float( + getattr(result, "mp2_correlation_hartree", None) + ), + "ccsd_correlation_hartree": _opt_float( + getattr(result, "ccsd_correlation_hartree", None) + ), + "ccsd_t_correction_hartree": _opt_float( + getattr(result, "ccsd_t_correction_hartree", None) + ), # Persisted so the saved-result card matches the live card (M-CLEAN # formatter-parity fix). Additive — absent on older results, where the # history card falls back exactly as before (CPU / no dipole / no @@ -187,6 +228,23 @@ def save_result( return dest +_COLLISION_SUFFIX_RE = re.compile(r"^(.*)_(\d+)$") + + +def _result_dir_sort_key(d: Path) -> tuple: + """Sort key that orders same-timestamp collision suffixes numerically. + + Directory names are ``___``, with a + ``_`` counter appended on same-microsecond collisions (N=1, 2, ...). + A plain string sort put ``..._10`` before ``..._2`` (lexicographic, not + numeric); split the trailing counter and sort on it as an int instead. + """ + m = _COLLISION_SUFFIX_RE.match(d.name) + if m: + return (m.group(1), int(m.group(2))) + return (d.name, -1) + + def list_results(results_dir: Optional[Path] = None) -> list: """Return result directories sorted newest-first. @@ -197,6 +255,7 @@ def list_results(results_dir: Optional[Path] = None) -> list: return [] return sorted( (d for d in base.iterdir() if d.is_dir() and (d / "result.json").exists()), + key=_result_dir_sort_key, reverse=True, ) @@ -357,7 +416,9 @@ def _append_molden_vibrations( ``normal_modes``). ``normal_modes`` is a list of length-N entries, each a list of per-atom (x, y, z) displacement triples. The ``[FR-COORD]`` block repeats the equilibrium geometry from - ``pyscf_mol_atom`` so the file is self-contained. + ``pyscf_mol_atom`` (converted Angstrom -> Bohr; the Molden spec + requires ``[FR-COORD]`` in Bohr regardless of ``[Atoms]``'s unit tag) + so the file is self-contained. """ with open(path, "a", encoding="utf-8") as fh: fh.write("\n[FREQ]\n") @@ -367,8 +428,9 @@ def _append_molden_vibrations( fh.write("\n[FR-COORD]\n") for sym, coords in pyscf_mol_atom: fh.write( - f"{sym} {float(coords[0]):.6f} {float(coords[1]):.6f} " - f"{float(coords[2]):.6f}\n" + f"{sym} {float(coords[0]) * _ANGSTROM_TO_BOHR:.6f} " + f"{float(coords[1]) * _ANGSTROM_TO_BOHR:.6f} " + f"{float(coords[2]) * _ANGSTROM_TO_BOHR:.6f}\n" ) fh.write("\n[FR-NORM-COORD]\n") @@ -710,6 +772,32 @@ def save_thumbnail(result_dir: Path, data: dict) -> None: except ImportError: return + # M9 audit fix (2026-07-14): only the matplotlib import itself was + # guarded — figure construction, text rendering, and fig.savefig() (a + # real filesystem write, so it can hit disk-full / permission errors) + # could all raise past this function despite the docstring's promise + # to silently skip "any error". Wrap the whole body so that promise + # actually holds; fig.close() still runs via finally regardless of + # where in the body a failure happened. + fig = None + try: + fig = _build_thumbnail_figure(plt, data) + fig.savefig( + str(result_dir / "thumbnail.png"), + dpi=144, + bbox_inches="tight", + facecolor=fig.get_facecolor(), + pad_inches=0.05, + ) + except Exception: + pass + finally: + if fig is not None: + plt.close(fig) + + +def _build_thumbnail_figure(plt: Any, data: dict) -> Any: + """Build (but don't save) the thumbnail matplotlib Figure for :func:`save_thumbnail`.""" _colors: dict = { "single_point": ("#2563eb", "#dbeafe"), "geometry_opt": ("#7c3aed", "#ede9fe"), @@ -809,13 +897,4 @@ def save_thumbnail(result_dir: Path, data: dict) -> None: transform=ax.transAxes, ) - try: - fig.savefig( - str(result_dir / "thumbnail.png"), - dpi=144, - bbox_inches="tight", - facecolor=bg, - pad_inches=0.05, - ) - finally: - plt.close(fig) + return fig diff --git a/quantui/session_calc.py b/quantui/session_calc.py index 8a8825a..d3fb065 100644 --- a/quantui/session_calc.py +++ b/quantui/session_calc.py @@ -230,7 +230,7 @@ def run_in_session( ``'6-31G'``, ``'6-31G*'``, ``'cc-pVDZ'``). Default: ``'6-31G'``. verbose: PySCF verbosity level (0 = silent … 9 = very detailed). Level 3 prints per-iteration SCF energies; level 4 adds - convergence diagnostics. Default: 3. + convergence diagnostics. Default: 4. progress_stream: Optional writable text stream. All PySCF output during the calculation is written here. Pass a widget-backed stream (e.g. ``_WidgetStream``) in the notebook for live display; @@ -329,11 +329,21 @@ def _run_session_calc_body( elif method_upper == "UHF": mf = scf.UHF(mol) elif method_upper == "MP2": - mf = scf.RHF(mol) # MP2 runs on top of RHF + # ``scf.RHF(mol)`` is a factory: for a closed-shell molecule + # (mol.spin == 0) it returns a true RHF object; for an open-shell + # molecule it auto-dispatches to ROHF instead (verified against + # PySCF's own factory behavior — this is not a QuantUI branch). + # ``mp.MP2(mf)`` below then further auto-dispatches: RMP2 on an + # RHF reference, UMP2 (ROHF-based) on an ROHF reference. Both are + # standard, well-defined methods; MP2 is not restricted to + # closed-shell input here. + mf = scf.RHF(mol) elif method_upper in ("CCSD", "CCSD(T)"): - # Coupled cluster builds on an RHF reference (M8.1). The correlation - # energy (and optional perturbative-triples correction) is added - # post-SCF below. + # Same auto-dispatch as MP2 above: scf.RHF(mol) yields RHF for + # closed-shell input and ROHF for open-shell input, and + # cc.CCSD(mf) below correspondingly dispatches to RCCSD or + # ROHF-based UCCSD. The correlation energy (and optional + # perturbative-triples correction) is added post-SCF below. mf = scf.RHF(mol) else: # DFT: resolve alias then auto-select RKS / UKS. ``resolve_xc`` @@ -487,26 +497,32 @@ def _to_numpy_array(arr: Any) -> Any: mulliken_charges: Optional[List[float]] = None dipole_moment_debye: Optional[float] = None - if method_upper != "UHF": - try: - # gpu4pyscf doesn't implement population analysis on the GPU object - # (``mf.mulliken_pop`` is NotImplemented), so on a GPU-offloaded run - # fall back to the host (CPU) object via ``to_cpu()``. ``chg`` is - # then host NumPy; _to_numpy_array also covers the CuPy case. - mf_pop = mf - if not callable(getattr(mf, "mulliken_pop", None)) and callable( - getattr(mf, "to_cpu", None) - ): - mf_pop = mf.to_cpu() - _, chg = mf_pop.mulliken_pop(verbose=0) - mulliken_charges = [float(c) for c in _to_numpy_array(chg)] - except Exception as exc: - logger.debug("Mulliken population extraction failed: %s", exc) - try: - dip = _to_numpy_array(mf.dip_moment(verbose=0)) - dipole_moment_debye = float(_np.linalg.norm(dip)) - except Exception as exc: - logger.debug("Dipole moment extraction failed: %s", exc) + # M3 audit fix (2026-07-14): both mf.mulliken_pop() and mf.dip_moment() + # are well-defined and work correctly for a genuine UHF object (verified + # empirically against PySCF) — the previous ``method_upper != "UHF"`` + # guard around this whole block was an unnecessary restriction that + # left the result card blank for both properties on every UHF run, + # while UKS (open-shell DFT) went through the identical extraction + # successfully. + try: + # gpu4pyscf doesn't implement population analysis on the GPU object + # (``mf.mulliken_pop`` is NotImplemented), so on a GPU-offloaded run + # fall back to the host (CPU) object via ``to_cpu()``. ``chg`` is + # then host NumPy; _to_numpy_array also covers the CuPy case. + mf_pop = mf + if not callable(getattr(mf, "mulliken_pop", None)) and callable( + getattr(mf, "to_cpu", None) + ): + mf_pop = mf.to_cpu() + _, chg = mf_pop.mulliken_pop(verbose=0) + mulliken_charges = [float(c) for c in _to_numpy_array(chg)] + except Exception as exc: + logger.debug("Mulliken population extraction failed: %s", exc) + try: + dip = _to_numpy_array(mf.dip_moment(verbose=0)) + dipole_moment_debye = float(_np.linalg.norm(dip)) + except Exception as exc: + logger.debug("Dipole moment extraction failed: %s", exc) # MO arrays for orbital visualization (non-fatal if extraction fails). # Uses the same ``_to_numpy_array`` CuPy→host helper defined above diff --git a/quantui/tddft_calc.py b/quantui/tddft_calc.py index 1660652..c344604 100644 --- a/quantui/tddft_calc.py +++ b/quantui/tddft_calc.py @@ -132,6 +132,20 @@ def run_tddft_calc( TD calculation fails, excitation lists are empty and a warning is written to progress_stream — no exception is raised. """ + # Post-HF methods (MP2/CCSD/CCSD(T)) have no special-casing below — + # without this guard, method='CCSD' silently falls into the DFT + # branch (sets mf.xc = "CCSD") and fails deep inside PySCF with a + # cryptic "LibXCFunctional: name 'CCSD' not found" instead of a clear + # message. TD-DFT/TDHF is not defined for these methods here. + from . import config as _config + + if method.strip().upper() in _config.POST_HF_METHODS: + raise ValueError( + f"'{method}' is a post-HF method and cannot be used for " + "TD-DFT/UV-Vis — use RHF/UHF (TDHF) or a DFT functional " + "instead." + ) + try: from pyscf import dft, gto, scf except ImportError as exc: diff --git a/quantui/utils.py b/quantui/utils.py index d7f383c..d8fee58 100644 --- a/quantui/utils.py +++ b/quantui/utils.py @@ -15,8 +15,14 @@ from . import config -# Configure logging -logging.basicConfig(level=getattr(logging, config.LOG_LEVEL), format=config.LOG_FORMAT) +# M14 audit fix (2026-07-14): logging.basicConfig() configures the *root* +# logger process-wide the moment this module is imported (which +# quantui/__init__ always does) — a library must never do this, since it +# silently hijacks/duplicates whatever logging setup the host application +# or notebook already has. quantui/__init__.py already attaches a +# NullHandler to the package logger; host apps that want console output +# can call logging.basicConfig() themselves, optionally using +# config.LOG_LEVEL / config.LOG_FORMAT as defaults. logger = logging.getLogger(__name__) diff --git a/tests/test_app.py b/tests/test_app.py index 9dac227..c17e403 100644 --- a/tests/test_app.py +++ b/tests/test_app.py @@ -392,6 +392,30 @@ def test_empty_write_is_noop(self): cap.write("") assert cap.getvalue() == "" + def test_close_is_noop_and_does_not_raise(self): + """Regression (found via the L6 audit fix's Python 3.9 CI matrix): + ase.utils.IOContext.openfile() — used by BFGS(..., logfile=...) in + optimizer.py / pes_scan.py — checks hasattr(file, "close") to decide + whether *file* is an already-open stream it should leave alone vs. a + path string it should open() itself. ase==3.26.0 (the newest version + pip resolves for Python 3.9) enforces this strictly and raised + TypeError for a _LogCapture instance, which had no close() method; + ase==3.29.0 (resolved for 3.10/3.11) happened to tolerate it via a + later refactor, masking the gap until 3.9 was added to CI. + """ + cap, _ = self._make_capture() + cap.close() # Must not raise + assert hasattr(cap, "close") + + def test_satisfies_ase_openfile_already_open_contract(self): + """Directly exercises the exact duck-typing check ASE performs.""" + cap, _ = self._make_capture() + assert hasattr(cap, "close"), ( + "ase.utils.IOContext.openfile() treats any object without a " + "'close' attribute as a path to open() itself, which fails for " + "a non-path file-like object like _LogCapture" + ) + # --------------------------------------------------------------------------- # _do_run dispatch @@ -642,6 +666,89 @@ def test_xyz_no_molecule_shows_error(self): app._on_export_xyz(None) assert "molecule" in app.struct_export_status.value.lower() + def test_xyz_filename_sanitizes_basis_with_asterisk(self, tmp_path): + """M11 audit fix (2026-07-14): a basis like "6-31G*" embedded + verbatim in a filename is invalid on Windows ("*" is a reserved + character there) and glob-hostile on POSIX. The exported filename + must not contain "*". + """ + app = QuantUIApp() + app._set_molecule(_water()) + app._last_result_dir = tmp_path + app.basis_dd.value = "6-31G*" + + app._on_export_xyz(None) + + xyz_files = list(tmp_path.glob("*.xyz")) + assert len(xyz_files) == 1 + assert "*" not in xyz_files[0].name + assert "Error" not in app.struct_export_status.value + + +class TestExportScriptCallback: + """_on_export (standalone PySCF script) sanitizes its filename too.""" + + def test_script_filename_sanitizes_basis_with_asterisk(self, tmp_path, monkeypatch): + """M11 audit fix (2026-07-14): same filename-sanitization bug as + the XYZ export, for the "Export Script" button — the script is + written to a bare relative filename in the current directory. + """ + monkeypatch.chdir(tmp_path) + app = QuantUIApp() + app._set_molecule(_water()) + app.basis_dd.value = "6-31G*" + + app._on_export(None) + + py_files = list(tmp_path.glob("*.py")) + assert len(py_files) == 1 + assert "*" not in py_files[0].name + assert "Error" not in app.export_status.value + + +class TestUpdateNotesBoldRendering: + """_update_notes converts every **bold** span, not just the first. + + Regression (M12 audit fix, 2026-07-14): the old implementation was + ``notes.replace("**", "", 1).replace("**", "", 1)`` — string + .replace(..., 1) only touches the FIRST occurrence in the whole + string, so only the first "**bold**" pair converted; every later one + (get_educational_notes() typically returns 2-3 separate + "**Label**: description" paragraphs) kept its literal "**" markers + and leaked into the rendered panel instead of rendering bold. + """ + + def test_multiple_bold_spans_all_converted(self, monkeypatch): + # UHF + 6-31G* + multiplicity 2 -> 3 separate "**Label**" spans + # in get_educational_notes()'s output. + # + # `with app.notes_output: display(HTML(...))` only populates the + # widget's `.outputs` under a live IPython display hook, which + # isn't present in a plain pytest process — so this test captures + # what's passed to `display()` directly instead of inspecting + # `notes_output.outputs` (the pattern used by tests of the + # newer `_set_html_output` atomic-swap helper, which manipulates + # `.outputs` directly and doesn't have this limitation). + import quantui.app_runflow as app_runflow + + captured: list = [] + monkeypatch.setattr(app_runflow, "display", lambda obj: captured.append(obj)) + + app = QuantUIApp() + mol = Molecule(["O", "H"], [[0.0, 0.0, 0.0], [0.96, 0.0, 0.0]], multiplicity=2) + app._set_molecule(mol) + app.method_dd.value = "UHF" + app.basis_dd.value = "6-31G*" + + captured.clear() # drop any renders triggered by the value changes above + app._update_notes() + + assert len(captured) == 1 + html = captured[0].data + assert "**" not in html, f"literal ** markers leaked into rendered HTML: {html}" + assert html.count("") >= 2 + assert html.count("") == html.count("") + class TestExportMoleculeAndLabel: """_export_molecule_and_label returns correct molecule and labels.""" @@ -1536,6 +1643,29 @@ def test_export_button_and_status_exist(self): assert isinstance(app._vib_export_status, widgets.HTML) assert app._vib_export_status.value == "" + def test_export_bad_mode_index_chains_original_exception(self): + """L audit fix (ruff B904): build_vib_export_html's py3Dmol-fallback + path must chain the original IndexError via `raise ... from exc` + when displacements[mode_number - 1] is out of range, not swallow it. + """ + from types import SimpleNamespace + + from quantui.app_visualization import build_vib_export_html + from quantui.viz_backend_router import BackendAvailability + + if not BackendAvailability.from_environment().py3dmol: + pytest.skip("py3Dmol not available for export fallback test") + + freq_stub = SimpleNamespace(displacements=[[[0.1, 0.0, 0.0]]]) + app_stub = SimpleNamespace( + _last_vib_freq_result=freq_stub, + _last_vib_molecule=self._water(), + _viz_availability=BackendAvailability(py3dmol=True, plotlymol=False), + ) + with pytest.raises(ValueError) as exc_info: + build_vib_export_html(app_stub, mode_number=5) # out of range + assert isinstance(exc_info.value.__cause__, IndexError) + def test_export_without_vib_state_shows_error_status(self, tmp_path, monkeypatch): monkeypatch.setenv("QUANTUI_RESULTS_DIR", str(tmp_path)) app = QuantUIApp() @@ -2430,7 +2560,7 @@ def test_render_orbital_isosurface_saves_cube_to_disk(self, tmp_path): captured: dict[str, object] = {} - def _fake_generate(_atom, _basis, _coeff, _idx, out_path): + def _fake_generate(_atom, _basis, _coeff, _idx, out_path, **_kwargs): captured["path"] = out_path out_path.write_text("cube", encoding="utf-8") return out_path @@ -2475,7 +2605,7 @@ def test_render_orbital_isosurface_py3dmol_path(self, tmp_path): captured: dict[str, object] = {} - def _fake_generate(_atom, _basis, _coeff, _idx, out_path): + def _fake_generate(_atom, _basis, _coeff, _idx, out_path, **_kwargs): captured["path"] = out_path out_path.write_text("cube", encoding="utf-8") return out_path diff --git a/tests/test_app_formatters.py b/tests/test_app_formatters.py index 6df05b5..d16fc0a 100644 --- a/tests/test_app_formatters.py +++ b/tests/test_app_formatters.py @@ -201,6 +201,46 @@ def test_format_nmr_result_warns_on_small_basis(): assert "qualitative NMR only" in html +def test_format_nmr_result_warns_on_fallback_reference(): + """M4 audit fix (2026-07-14): fallback reference constants are surfaced. + + Regression: an untabulated method/basis silently substituted the + B3LYP/6-31G* reference constants with no indication anywhere in the + result card — the student had no way to know their chemical shifts + might be off by a few ppm. + """ + result = _NMRStub( + converged=True, + reference_compound="TMS", + atom_symbols=["O", "H", "H"], + chemical_shifts_ppm={1: 4.5, 2: 4.5}, + formula="H2O", + method="CAM-B3LYP", + basis="6-31G*", + reference_key="B3LYP/6-31G*", + is_fallback_reference=True, + ) + html = format_nmr_result(result) + assert "No calibrated TMS reference" in html + assert "CAM-B3LYP/6-31G*" in html + + +def test_format_nmr_result_no_warning_for_exact_reference_match(): + result = _NMRStub( + converged=True, + reference_compound="TMS", + atom_symbols=["O", "H", "H"], + chemical_shifts_ppm={1: 4.5, 2: 4.5}, + formula="H2O", + method="B3LYP", + basis="6-31G*", + reference_key="B3LYP/6-31G*", + is_fallback_reference=False, + ) + html = format_nmr_result(result) + assert "No calibrated TMS reference" not in html + + def test_format_pes_scan_result_reports_range_and_convergence(): result = SimpleNamespace( converged_all=True, diff --git a/tests/test_benchmarks.py b/tests/test_benchmarks.py index 32fa6cf..12b0eff 100644 --- a/tests/test_benchmarks.py +++ b/tests/test_benchmarks.py @@ -18,6 +18,7 @@ BENCHMARK_SUITE, BenchmarkStep, CalibrationResult, + _count_electrons, load_last_calibration, run_calibration, ) @@ -38,6 +39,30 @@ not _PYSCF_AVAILABLE, reason="PySCF not installed (Linux/macOS/WSL only)" ) +# --------------------------------------------------------------------------- +# _count_electrons +# --------------------------------------------------------------------------- + + +class TestCountElectrons: + """L audit fix: _count_electrons used to hand-carry its own truncated + (Z=1-18) atomic-number table, silently falling back to carbon's Z=6 + for every heavier element (e.g. iodine, gold) instead of using the + shared, full-periodic-table config.ATOMIC_NUMBERS. + """ + + def test_light_elements_match_known_z(self): + assert _count_electrons(["O", "H", "H"], charge=0) == 10 # water + + def test_charge_is_subtracted(self): + assert _count_electrons(["O", "H", "H"], charge=1) == 9 + + def test_heavy_element_uses_real_atomic_number_not_carbon_fallback(self): + # Iodine: Z=53. The old 18-element table fell back to 6 for any + # element past Ar, undercounting badly for anything heavier. + assert _count_electrons(["I"], charge=0) == 53 + + # --------------------------------------------------------------------------- # BENCHMARK_SUITE contents # --------------------------------------------------------------------------- diff --git a/tests/test_calc_log.py b/tests/test_calc_log.py index 324b144..d7b38a3 100644 --- a/tests/test_calc_log.py +++ b/tests/test_calc_log.py @@ -128,3 +128,215 @@ def test_estimate_time_non_single_point_ignores_legacy_untyped_records( f"for 4 atoms), got {est_freq['seconds']:.1f} s — suggests the cost " "model isn't firing on legacy SP records" ) + + +# ============================================================================ +# Event log: append/prune race + throttled pruning (M8 audit fix, 2026-07-14) +# ============================================================================ + + +def test_log_event_basic_roundtrip(isolated_log_dir): + import quantui.calc_log as clog + + clog.log_event("test_event", "hello world", extra_field=42) + events = clog.get_recent_events(10) + assert len(events) == 1 + assert events[0]["event"] == "test_event" + assert events[0]["message"] == "hello world" + assert events[0]["extra_field"] == 42 + + +def test_prune_events_holds_lock_across_read_and_rewrite(isolated_log_dir): + """Concurrent log_event() calls must never lose an append to a + concurrent prune_events() rewrite. + + Regression: prune_events() used to read the file (acquiring and + releasing the module lock) and only later reacquire the lock to + rewrite it. An append landing in that gap got silently discarded when + the rewrite replaced the whole file with the (stale) filtered list + computed before that append happened. Fixed by making the read + + filter + rewrite one lock-held critical section. + """ + import threading + + import quantui.calc_log as clog + + n_threads = 8 + n_events_per_thread = 30 + + def worker(tid: int) -> None: + for i in range(n_events_per_thread): + clog.log_event("stress_test", f"thread {tid} event {i}") + + threads = [threading.Thread(target=worker, args=(t,)) for t in range(n_threads)] + for t in threads: + t.start() + for t in threads: + t.join() + + events = clog.get_recent_events(10_000) + expected = n_threads * n_events_per_thread + assert len(events) == expected, ( + f"expected {expected} events, got {len(events)} — " + f"{expected - len(events)} were lost to the append/prune race" + ) + + +def test_prune_not_triggered_on_every_single_event(isolated_log_dir, monkeypatch): + """log_event() must not call the full prune on every single append. + + Regression: log_event() called prune_events() (a full read + rewrite + of the whole file) after every single append — O(file size) work per + event, O(N^2) over a session. It should only prune periodically. + """ + import quantui.calc_log as clog + + prune_calls = {"n": 0} + real_prune = clog.prune_events + + def _counting_prune(*args, **kwargs): + prune_calls["n"] += 1 + return real_prune(*args, **kwargs) + + monkeypatch.setattr(clog, "prune_events", _counting_prune) + + n_events = clog._PRUNE_EVERY_N_EVENTS * 3 - 1 + for i in range(n_events): + clog.log_event("test_event", f"event {i}") + + # Should prune roughly every _PRUNE_EVERY_N_EVENTS calls, not once per event. + assert prune_calls["n"] < n_events + assert prune_calls["n"] == n_events // clog._PRUNE_EVERY_N_EVENTS + + +def test_prune_events_still_removes_old_entries(isolated_log_dir): + """The periodic-prune change must not break the actual 7-day TTL.""" + import json + from datetime import datetime, timedelta, timezone + + import quantui.calc_log as clog + + path = clog._event_path() + path.parent.mkdir(parents=True, exist_ok=True) + old_ts = (datetime.now(timezone.utc) - timedelta(days=10)).isoformat() + new_ts = datetime.now(timezone.utc).isoformat() + with open(path, "w", encoding="utf-8") as fh: + fh.write( + json.dumps({"timestamp": old_ts, "event": "old", "message": "stale"}) + "\n" + ) + fh.write( + json.dumps({"timestamp": new_ts, "event": "new", "message": "fresh"}) + "\n" + ) + + clog.prune_events() + + remaining = clog.get_recent_events(10) + assert len(remaining) == 1 + assert remaining[0]["event"] == "new" + + +def test_read_all_caches_between_unchanged_calls(isolated_log_dir): + """L audit fix: _read_all must not re-parse the perf log when it hasn't + changed since the last call — estimate_time() re-read the entire + (indefinitely-kept) file on every UI refresh even when nothing new + had been written. + """ + import quantui.calc_log as clog + + clog.log_calculation( + formula="H2", + n_atoms=2, + n_electrons=2, + method="RHF", + basis="STO-3G", + n_iterations=5, + elapsed_s=1.0, + converged=True, + n_basis=2, + n_cores=1, + calc_type="single_point", + ) + + first = clog._read_all(clog._perf_path()) + assert len(first) == 1 + + calls = {"n": 0} + real_loads = clog.json.loads + + def _counting_loads(s): + calls["n"] += 1 + return real_loads(s) + + clog.json.loads = _counting_loads + try: + second = clog._read_all(clog._perf_path()) + finally: + clog.json.loads = real_loads + + assert second == first + assert calls["n"] == 0 # cache hit: no re-parsing of the unchanged file + + +def test_read_all_cache_invalidates_on_new_write(isolated_log_dir): + """The caching fix above must not go stale after a genuine new write.""" + import quantui.calc_log as clog + + clog.log_calculation( + formula="H2", + n_atoms=2, + n_electrons=2, + method="RHF", + basis="STO-3G", + n_iterations=5, + elapsed_s=1.0, + converged=True, + n_basis=2, + n_cores=1, + calc_type="single_point", + ) + first = clog._read_all(clog._perf_path()) + assert len(first) == 1 + + clog.log_calculation( + formula="H2", + n_atoms=2, + n_electrons=2, + method="RHF", + basis="STO-3G", + n_iterations=6, + elapsed_s=2.0, + converged=True, + n_basis=2, + n_cores=1, + calc_type="single_point", + ) + second = clog._read_all(clog._perf_path()) + assert len(second) == 2 + + +def test_6_31gss_he_basis_count_matches_pyscf(isolated_log_dir): + """L audit fix: He under 6-31G** must have 5 basis functions, not 2. + + 6-31G** adds a p-polarization shell on H/He on top of 6-31G*'s bare + s-only He (2 bf), the same way 6-31G* already adds p on the heavy + atoms — so He should follow H's pattern (2 -> 5), not stay at 2. + Verified against ``pyscf.gto.M(atom="He", basis="6-31g**").nao == 5``. + """ + import quantui.calc_log as clog + + assert clog.count_basis_functions(["He"], "6-31G**") == 5 + + +def test_basis_function_table_internally_consistent(isolated_log_dir): + """H and He must carry equal counts in every basis in the lookup table. + + Both are period-1 elements with the same shell structure in every + basis set this table covers, so a basis that gives H and He different + counts indicates a transcription error (this caught the 6-31G** nit). + """ + import quantui.calc_log as clog + + for basis, table in clog._BASIS_FUNCTIONS.items(): + assert ( + table["H"] == table["He"] + ), f"{basis}: H={table['H']} but He={table['He']}, expected equal" diff --git a/tests/test_calculator.py b/tests/test_calculator.py index 777f661..42c7d59 100644 --- a/tests/test_calculator.py +++ b/tests/test_calculator.py @@ -88,6 +88,23 @@ def test_unsupported_method(self, water_molecule): with pytest.raises(ValueError, match="not supported"): PySCFCalculation(water_molecule, method="NONEXISTENT", basis="6-31G") + def test_mixed_case_method_preserved_canonically(self, water_molecule): + """Mixed-case SUPPORTED_METHODS entries must not be upper-cased. + + Regression: method.upper() turned "wB97X-D" into "WB97X-D", which + matches neither SUPPORTED_METHODS nor the xc-alias tables the + generated script depends on — silently breaking Export Script and + the educational-notes lookup for every mixed-case method name + (wB97X-D, CAM-B3LYP, M06-L, M06-2X, HSE06, PBE-D3). + """ + calc = PySCFCalculation(water_molecule, method="wB97X-D", basis="6-31G") + assert calc.method == "wB97X-D" + + def test_mixed_case_method_case_insensitive_input(self, water_molecule): + """A differently-cased spelling still resolves to the canonical form.""" + calc = PySCFCalculation(water_molecule, method="wb97x-d", basis="6-31G") + assert calc.method == "wB97X-D" + def test_nonstandard_basis_warning(self, water_molecule, caplog): """Test warning for non-standard basis set.""" calc = PySCFCalculation(water_molecule, method="RHF", basis="custom-basis") diff --git a/tests/test_cli.py b/tests/test_cli.py index cad6083..91a0630 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -10,6 +10,7 @@ import io import json +import subprocess import sys import pytest @@ -407,3 +408,59 @@ def _fake_open(url, *_args, **_kwargs): assert ok is True assert tool == "webbrowser" assert opened[0].startswith("file:") + + +class TestCliAvoidsGuiStackImport: + """M13 audit fix: ``import quantui.cli`` must not pull in ipywidgets. + + Importing any submodule of ``quantui`` always runs ``quantui/__init__.py`` + first (Python import semantics), and that module used to eagerly import + ``QuantUIApp`` (and, transitively, ``ipywidgets``/``IPython`` widget + machinery) even for callers that only want the pure-Python CLI helpers. + This must run in a subprocess — checking ``sys.modules`` in-process is + unreliable once other test modules in the same pytest session have + already imported ipywidgets/quantui.app themselves. + """ + + def test_import_cli_does_not_load_ipywidgets_or_app(self): + script = ( + "import sys\n" + "import quantui.cli\n" + "assert 'ipywidgets' not in sys.modules, 'ipywidgets was imported'\n" + "assert 'quantui.app' not in sys.modules, 'quantui.app was imported'\n" + "assert 'quantui.progress' not in sys.modules, " + "'quantui.progress was imported'\n" + "print('OK')\n" + ) + result = subprocess.run( + [sys.executable, "-c", script], + capture_output=True, + text=True, + timeout=60, + ) + assert result.returncode == 0, result.stdout + result.stderr + assert "OK" in result.stdout + + def test_lazy_attrs_still_resolve_on_demand(self): + script = ( + "import quantui\n" + "assert quantui.StepProgress(['a']).widget is not None\n" + "assert len(quantui.HELP_TOPICS) > 0\n" + "assert 'charge' in quantui.VALID_TOPICS\n" + "quantui.help_panel('charge')\n" + "assert quantui.QuantUIApp.__name__ == 'QuantUIApp'\n" + "print('OK')\n" + ) + result = subprocess.run( + [sys.executable, "-c", script], + capture_output=True, + text=True, + timeout=60, + ) + assert result.returncode == 0, result.stdout + result.stderr + assert "OK" in result.stdout + + def test_unknown_attribute_still_raises(self): + import quantui + + assert hasattr(quantui, "THIS_ATTR_DOES_NOT_EXIST") is False diff --git a/tests/test_export_molden.py b/tests/test_export_molden.py index 646bde1..a94633f 100644 --- a/tests/test_export_molden.py +++ b/tests/test_export_molden.py @@ -216,3 +216,40 @@ def test_writes_freq_block_alongside_orbitals_when_both_present(self, tmp_path): assert "[MO]" in text assert "[FREQ]" in text assert "[FR-NORM-COORD]" in text + + +@pyscf_only +class TestFrCoordUnits: + """[FR-COORD] must be Bohr per the Molden spec, regardless of [Atoms]'s + unit tag (theochem.ru.nl/molden/molden_format.html) — a Molden quirk. + pyscf_mol_atom is Angstrom throughout QuantUI, so the writer must + convert. Regression for a bug where [FR-COORD] was written directly + from (Angstrom) pyscf_mol_atom with no conversion, making it + inconsistent with [Atoms]/[GTO]/[MO] in the same file. + """ + + def test_fr_coord_is_bohr_not_angstrom(self, tmp_path): + _ANGSTROM_TO_BOHR = 1.8897261254578281 + atom_list = _water_atom_list() + out = save_molden( + tmp_path, + pyscf_mol_atom=atom_list, + pyscf_mol_basis="sto-3g", + frequencies_cm1=[1500.0], + normal_modes=[[[0.1, 0.0, 0.0], [-0.05, 0.0, 0.0], [-0.05, 0.0, 0.0]]], + ) + text = out.read_text(encoding="utf-8") + fr_coord_block = text.split("[FR-COORD]")[1].split("[FR-NORM-COORD]")[0] + lines = [ln for ln in fr_coord_block.strip().splitlines() if ln.strip()] + assert len(lines) == len(atom_list) + for line, (sym, coords) in zip(lines, atom_list): + parts = line.split() + assert parts[0] == sym + got = [float(p) for p in parts[1:4]] + expected = [c * _ANGSTROM_TO_BOHR for c in coords] + for g, e in zip(got, expected): + assert g == pytest.approx(e, abs=1e-5) + # Sanity: the Bohr value should NOT equal the raw Angstrom + # value (except for the trivial all-zero O atom). + if any(abs(c) > 1e-9 for c in coords): + assert got != pytest.approx(coords, abs=1e-3) diff --git a/tests/test_freq_calc.py b/tests/test_freq_calc.py index 38ac54a..da83182 100644 --- a/tests/test_freq_calc.py +++ b/tests/test_freq_calc.py @@ -175,6 +175,59 @@ def test_thermo_g_less_than_h(self): assert result.thermo.G_hartree < result.thermo.H_hartree +# ============================================================================ +# pyscf_mol_atom unit convention (H2 audit fix, 2026-07-14) +# ============================================================================ + + +class TestPyscfMolAtomUnits: + """``pyscf_mol_atom`` must be Angstrom, matching session_calc/optimizer. + + Regression for a bug where freq_calc built this field from PySCF's + internal ``mol._atom`` (always Bohr), while every consumer (Molden + export, cube generation, orbital replay) assumes Angstrom — silently + inflating exported geometries ~1.89x for Frequency results only. + """ + + @pyscf_only + @pytest.mark.slow + def test_pyscf_mol_atom_matches_input_geometry_in_angstrom(self): + from quantui.freq_calc import run_freq_calc + + molecule = _water() + result = run_freq_calc(molecule, method="RHF", basis="STO-3G") + assert result.pyscf_mol_atom is not None + for (sym, coords), (orig_sym, orig_coords) in zip( + result.pyscf_mol_atom, zip(molecule.atoms, molecule.coordinates) + ): + assert sym == orig_sym + for c, orig_c in zip(coords, orig_coords): + assert c == pytest.approx(orig_c, abs=1e-9) + + +# ============================================================================ +# Post-HF method guard (M2 audit fix, 2026-07-14) +# ============================================================================ + + +class TestRunFreqCalcPostHfGuard: + """Post-HF methods raise a clear ValueError instead of a cryptic LibXC error. + + Regression: run_freq_calc() had no special-casing for MP2/CCSD/CCSD(T) + — the SCF-selection branch silently treated them as a DFT xc functional + (mf.xc = "CCSD"), failing deep inside PySCF with "LibXCFunctional: name + 'CCSD' not found" instead of a clear message. The guard fires before any + PySCF import, so it needs neither PySCF nor ASE. + """ + + @pytest.mark.parametrize("method", ["MP2", "CCSD", "CCSD(T)"]) + def test_post_hf_method_raises_value_error(self, method): + from quantui.freq_calc import run_freq_calc + + with pytest.raises(ValueError, match="post-HF"): + run_freq_calc(_water(), method=method, basis="STO-3G") + + # ============================================================================ # IR intensities — PySCF required # ============================================================================ @@ -224,5 +277,118 @@ def test_ir_intensities_physically_reasonable(self): assert max(result.ir_intensities) > 1.0 +# ============================================================================ +# IR-intensity inner-loop dm0/method dispatch (M5 audit fix, 2026-07-14) +# ============================================================================ + + +class TestIrIntensityUhfClosedShellDispatch: + """The inner displaced-SCF loop must dispatch on dm0's actual shape. + + Regression: the serial (_displaced_scf_dipole in freq_calc.py) and + parallel (run_displaced_scf in freq_ir_workers.py) inner loops both + picked RHF/UHF based on mol.spin == 0 alone. That agrees with the + parent SCF's actual type only when the user's method choice matches + the molecule's natural spin state. Explicitly selecting UHF for a + closed-shell molecule (mol.spin == 0, but the parent mf — and its + dm0 — is still UHF-shaped, a legitimate technique for probing + symmetry-broken solutions) used to raise a shape-mismatch ValueError + deep in PySCF, silently dropping IR intensities for the whole run + (caught by the broad except around the entire IR-intensity block). + """ + + @pyscf_only + @pytest.mark.slow + def test_uhf_on_closed_shell_molecule_still_gets_ir_intensities(self): + from quantui.freq_calc import run_freq_calc + + result = run_freq_calc(_water(), method="UHF", basis="STO-3G") + assert result.ir_intensities, ( + "UHF on a closed-shell molecule should still produce IR " + "intensities via the serial inner-loop dispatch fix" + ) + assert len(result.ir_intensities) == len(result.frequencies_cm1) + + def test_freq_ir_workers_dispatches_on_dm0_shape_not_spin(self): + """Unit-level check of the parallel worker's dispatch logic directly. + + Builds a UHF parent for a closed-shell (spin=0) molecule — the + exact scenario that used to crash — and confirms the worker + picks UHF (matching the (2, nao, nao) dm0) rather than RHF + (which would raise on that dm0 shape). + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto, scf + + from quantui.freq_ir_workers import init_worker, run_displaced_scf + + atom_str = "O 0 0 0.119; H 0 0.763 -0.477; H 0 -0.763 -0.477" + mol = gto.M(atom=atom_str, basis="sto-3g", spin=0, charge=0, verbose=0) + mf = scf.UHF(mol) + mf.kernel() + dm0 = mf.make_rdm1() + assert dm0.ndim == 3 # UHF-shaped despite mol.spin == 0 + + tmp = tempfile.NamedTemporaryFile(delete=False, suffix=".pkl") + try: + pickle.dump(dm0, tmp) + tmp.close() + init_worker(atom_str, "sto-3g", 0, 0, None, tmp.name, 1) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + dip = run_displaced_scf(coords) # must not raise + assert len(dip) == 3 + finally: + os.unlink(tmp.name) + + @pyscf_only + @pytest.mark.slow + def test_parallel_failure_falls_back_to_serial(self, monkeypatch): + """A parallel-path failure must fall back to serial, not give up. + + Regression: run_displaced_scf's own docstring claims "The + freq_calc driver catches such failures and falls back to the + serial loop so the user's calc still completes" — but no such + fallback existed; any single worker failure propagated out to + the broad except around the whole IR-intensity block, dropping + IR intensities entirely. Forces the parallel path to be + selected and to fail immediately, then checks IR intensities + still come out via the serial fallback. + """ + import concurrent.futures as cf + import io as _io + + import quantui.freq_ir_workers as ir_workers + from quantui.freq_calc import run_freq_calc + + monkeypatch.setattr(ir_workers, "parallel_enabled_for_run", lambda **kw: True) + + class _FailingExecutor: + def __init__(self, *a, **kw): + pass + + def __enter__(self): + raise RuntimeError("simulated worker pool failure") + + def __exit__(self, *a): + return False + + monkeypatch.setattr(cf, "ProcessPoolExecutor", _FailingExecutor) + + buf = _io.StringIO() + result = run_freq_calc( + _water(), method="RHF", basis="STO-3G", progress_stream=buf + ) + log = buf.getvalue() + assert "falling back to serial" in log + assert result.ir_intensities, ( + "IR intensities should still be populated via the serial " + "fallback after a simulated parallel-path failure" + ) + + if __name__ == "__main__": pytest.main([__file__, "-v", "--tb=short"]) diff --git a/tests/test_molecule.py b/tests/test_molecule.py index f4dfa8a..037497d 100644 --- a/tests/test_molecule.py +++ b/tests/test_molecule.py @@ -378,8 +378,11 @@ def test_parse_non_numeric_coordinates(self): xyz_text = """H 0.0 0.0 abc H 0.0 0.0 0.74""" - with pytest.raises(ValueError, match="Could not parse coordinates"): + with pytest.raises(ValueError, match="Could not parse coordinates") as exc_info: parse_xyz_input(xyz_text) + # L audit fix (ruff B904): the original ValueError from float() + # must be chained via `raise ... from e`. + assert isinstance(exc_info.value.__cause__, ValueError) def test_parse_negative_coordinates(self): """Test parsing negative coordinates.""" @@ -393,6 +396,57 @@ def test_parse_negative_coordinates(self): assert coords[2][0] == -0.757 +class TestXyzHeaderTitleLine: + """XYZ-file-format header (count + title) edge cases (H5 audit fix, 2026-07-14). + + Regression: blank/comment title lines were stripped by the parser's + blank/comment filtering BEFORE header detection ran, which shifted the + first real atom line into the "title" slot — silently dropping the + first atom for any standard XYZ file whose title line was blank or a + "#"/"!" comment (both very common: many tools write an empty title). + """ + + def test_blank_title_line_keeps_all_atoms(self): + xyz_text = "3\n\nO 0.0 0.0 0.0\nH 0.757 0.587 0.0\nH -0.757 0.587 0.0" + atoms, coords = parse_xyz_input(xyz_text) + assert atoms == ["O", "H", "H"] + assert len(coords) == 3 + assert coords[0] == [0.0, 0.0, 0.0] + + def test_comment_title_line_keeps_all_atoms(self): + xyz_text = ( + "3\n# Water molecule\nO 0.0 0.0 0.0\n" + "H 0.757 0.587 0.0\nH -0.757 0.587 0.0" + ) + atoms, coords = parse_xyz_input(xyz_text) + assert atoms == ["O", "H", "H"] + assert len(coords) == 3 + + def test_bang_comment_title_line_keeps_all_atoms(self): + xyz_text = ( + "3\n! Water molecule\nO 0.0 0.0 0.0\n" + "H 0.757 0.587 0.0\nH -0.757 0.587 0.0" + ) + atoms, coords = parse_xyz_input(xyz_text) + assert atoms == ["O", "H", "H"] + + def test_free_text_title_line_still_works(self): + xyz_text = ( + "3\nWater molecule\nO 0.0 0.0 0.0\n" + "H 0.757 0.587 0.0\nH -0.757 0.587 0.0" + ) + atoms, coords = parse_xyz_input(xyz_text) + assert atoms == ["O", "H", "H"] + + def test_leading_blank_line_before_count_still_works(self): + xyz_text = ( + "\n3\ntitle\nO 0.0 0.0 0.0\n" + "H 0.757 0.587 0.0\nH -0.757 0.587 0.0" + ) + atoms, coords = parse_xyz_input(xyz_text) + assert atoms == ["O", "H", "H"] + + class TestEnhancedXYZParser: """Test enhanced XYZ parser features.""" @@ -435,12 +489,29 @@ def test_parse_with_inline_comments(self): assert len(coords) == 3 assert coords[0] == [0.0, 0.0, 0.0] - def test_parse_minimum_atoms_error(self): - """Test error for single atom (need at least 2).""" - xyz_text = """H 0.0 0.0 0.0""" + def test_parse_single_atom_succeeds(self): + """L13 audit fix: single-atom input must parse, not be rejected. - with pytest.raises(ValueError, match="Too Few Atoms"): - parse_xyz_input(xyz_text) + Atomic calculations (e.g. a lone Ar atom) are legitimate PySCF + targets; parse_xyz_input used to hard-reject anything with fewer + than 2 atoms. + """ + xyz_text = """Ar 0.0 0.0 0.0""" + + atoms, coords = parse_xyz_input(xyz_text) + + assert atoms == ["Ar"] + assert coords == [[0.0, 0.0, 0.0]] + + def test_parse_single_atom_with_xyz_header_succeeds(self): + xyz_text = """1 +Single argon atom +Ar 0.0 0.0 0.0""" + + atoms, coords = parse_xyz_input(xyz_text) + + assert atoms == ["Ar"] + assert coords == [[0.0, 0.0, 0.0]] def test_parse_invalid_atom_symbol_with_suggestion(self): """Test that invalid atom symbol provides helpful suggestion.""" diff --git a/tests/test_nmr_calc.py b/tests/test_nmr_calc.py index 1578100..148ddfa 100644 --- a/tests/test_nmr_calc.py +++ b/tests/test_nmr_calc.py @@ -151,6 +151,44 @@ def test_c_reference_plausible(self): assert 150.0 < NMR_DEFAULT_REFERENCE["C"] < 220.0 +class TestResolveNmrReference: + """M4 audit fix (2026-07-14): record which reference was actually used. + + Regression: any method/basis combination not in + config.NMR_REFERENCE_SHIELDINGS silently used the B3LYP/6-31G* + constants with no record anywhere that a substitution happened — + chemical shifts could be off by a few ppm with no indication to the + student. resolve_nmr_reference() (and the NMRResult.reference_key / + is_fallback_reference fields it feeds) makes that substitution visible. + No PySCF needed — pure lookup logic. + """ + + def test_exact_match_is_not_fallback(self): + from quantui.nmr_calc import resolve_nmr_reference + + ref_map, key, is_fallback = resolve_nmr_reference("B3LYP", "6-31G*") + assert key == "B3LYP/6-31G*" + assert is_fallback is False + assert ref_map == {"H": 31.72, "C": 183.71} + + def test_unlisted_combination_is_fallback(self): + from quantui.nmr_calc import resolve_nmr_reference + + ref_map, key, is_fallback = resolve_nmr_reference("CAM-B3LYP", "6-31G*") + assert is_fallback is True + assert key == "B3LYP/6-31G*" + from quantui.config import NMR_DEFAULT_REFERENCE + + assert ref_map == NMR_DEFAULT_REFERENCE + + def test_lookup_is_case_insensitive(self): + from quantui.nmr_calc import resolve_nmr_reference + + ref_map, key, is_fallback = resolve_nmr_reference("rhf", "6-31g*") + assert is_fallback is False + assert key == "RHF/6-31G*" + + # --------------------------------------------------------------------------- # run_nmr_calc — PySCF-gated # --------------------------------------------------------------------------- @@ -211,6 +249,88 @@ def test_shielding_iso_length_matches_atoms(self): result = run_nmr_calc(mol, method="RHF", basis="STO-3G") assert len(result.shielding_iso_ppm) == len(list(mol.atoms)) + @pyscf_only + @pytest.mark.slow + def test_exact_match_reference_not_flagged_as_fallback(self): + """M4 audit fix: RHF/STO-3G is tabulated -> not a fallback.""" + result = run_nmr_calc(_water(), method="RHF", basis="STO-3G") + assert result.reference_key == "RHF/STO-3G" + assert result.is_fallback_reference is False + + @pyscf_only + @pytest.mark.slow + def test_untabulated_combination_flagged_as_fallback(self): + """M4 audit fix: an untabulated method/basis is flagged, not silent.""" + result = run_nmr_calc(_water(), method="CAM-B3LYP", basis="6-31G*") + assert result.is_fallback_reference is True + assert result.reference_key == "B3LYP/6-31G*" + + +class TestNmrCompatPatchIdempotency: + """L audit fix: the pyscf.prop.nmr compat patches must apply once per + process, not redefine + reassign the same closures on every NMR call. + """ + + @pyscf_only + def test_patch_functions_are_not_redefined_on_repeat_calls(self): + import pyscf.prop.nmr.rhf as _prop_nmr_rhf + import pyscf.prop.nmr.rks as _prop_nmr_rks + + from quantui.nmr_calc import _ensure_nmr_compat_patches_applied + + _ensure_nmr_compat_patches_applied() + gen_vind_first = _prop_nmr_rhf.gen_vind + get_vxc_giao_first = _prop_nmr_rks.get_vxc_giao + + _ensure_nmr_compat_patches_applied() + # Same function object — the second call must be a no-op, not a + # fresh `def` + reassignment (which would install a *different* + # (if behaviorally identical) closure each time). + assert _prop_nmr_rhf.gen_vind is gen_vind_first + assert _prop_nmr_rks.get_vxc_giao is get_vxc_giao_first + + @pyscf_only + def test_patched_functions_carry_version_sentinel(self): + import pyscf.prop.nmr.rhf as _prop_nmr_rhf + import pyscf.prop.nmr.rks as _prop_nmr_rks + + from quantui.nmr_calc import ( + _NMR_COMPAT_PATCH_VERSION, + _NMR_PATCH_VERSION_ATTR, + _ensure_nmr_compat_patches_applied, + ) + + _ensure_nmr_compat_patches_applied() + assert ( + getattr(_prop_nmr_rhf.gen_vind, _NMR_PATCH_VERSION_ATTR, None) + == _NMR_COMPAT_PATCH_VERSION + ) + assert ( + getattr(_prop_nmr_rks.get_vxc_giao, _NMR_PATCH_VERSION_ATTR, None) + == _NMR_COMPAT_PATCH_VERSION + ) + + +# ============================================================================ +# Post-HF method guard (M2 audit fix, 2026-07-14) +# ============================================================================ + + +class TestRunNmrCalcPostHfGuard: + """Post-HF methods raise a clear ValueError instead of a cryptic LibXC error. + + Regression: run_nmr_calc() had no special-casing for MP2/CCSD/CCSD(T) + — the SCF-selection branch silently treated them as a DFT xc functional + (mf.xc = "CCSD"), failing deep inside PySCF with "LibXCFunctional: name + 'CCSD' not found" instead of a clear message. The guard fires before any + PySCF import, so it needs no PySCF. + """ + + @pytest.mark.parametrize("method", ["MP2", "CCSD", "CCSD(T)"]) + def test_post_hf_method_raises_value_error(self, method): + with pytest.raises(ValueError, match="post-HF"): + run_nmr_calc(_water(), method=method, basis="STO-3G") + if __name__ == "__main__": pytest.main([__file__, "-v", "--tb=short"]) diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py index 8ce2a49..2c0429e 100644 --- a/tests/test_optimizer.py +++ b/tests/test_optimizer.py @@ -183,6 +183,25 @@ def test_config_exports_match(self): assert config.DEFAULT_FMAX == DEFAULT_FMAX assert config.DEFAULT_OPT_STEPS == DEFAULT_OPT_STEPS + def test_bohr_to_angstrom_matches_shared_config_constant(self): + """L audit fix: optimizer.py used to hand-type its own Bohr->Angstrom + literal (0.529177249, a plain CODATA-2018 rounding) that didn't + match freq_calc.py's separately hand-typed literal + (0.52917721092, pyscf.data.nist.BOHR). Both now import + config.BOHR_TO_ANGSTROM, so force/gradient unit conversions in + optimize_geometry are computed with the same constant PySCF itself + uses internally. + """ + import quantui.optimizer as opt_mod + from quantui import config + from quantui.results_storage import _ANGSTROM_TO_BOHR + + assert opt_mod._BOHR_TO_ANG == config.BOHR_TO_ANGSTROM + assert config.BOHR_TO_ANGSTROM == pytest.approx(0.52917721092, abs=1e-11) + assert _ANGSTROM_TO_BOHR == pytest.approx( + 1.0 / config.BOHR_TO_ANGSTROM, abs=1e-15 + ) + # ============================================================================ # Import-guard tests — all platforms @@ -190,37 +209,61 @@ def test_config_exports_match(self): class TestOptimizeGeometryImportGuards: + """Patch quantui.optimizer's OWN ``ASE_AVAILABLE`` name directly rather + than reloading the module after monkeypatching ase_bridge's copy. + + ``optimize_geometry()`` reads the module-local ``ASE_AVAILABLE`` binding + (imported via ``from .ase_bridge import ASE_AVAILABLE`` at module load + time) — monkeypatching ``ase_bridge.ASE_AVAILABLE`` doesn't affect that + binding at all, which is why the original version of this test used + ``importlib.reload(opt_mod)`` to re-execute the import with the patched + value in scope. But ``monkeypatch`` only reverts attributes it directly + set (``ase_bridge.ASE_AVAILABLE``) — it has no way to know the reload + also needs undoing, so ``quantui.optimizer.ASE_AVAILABLE`` stayed stuck + at ``False`` for the rest of the test session once this test ran, + silently breaking every later test that called ``optimize_geometry()`` + expecting ASE to be available. Patching the local name directly avoids + the reload (and the leak) entirely. + """ + def test_raises_when_ase_unavailable(self, monkeypatch): - import quantui.ase_bridge as bridge + import quantui.optimizer as opt_mod - monkeypatch.setattr(bridge, "ASE_AVAILABLE", False) + monkeypatch.setattr(opt_mod, "ASE_AVAILABLE", False) - import importlib + with pytest.raises(ImportError, match="pip install"): + opt_mod.optimize_geometry(_h2()) + def test_error_message_is_actionable(self, monkeypatch): import quantui.optimizer as opt_mod - importlib.reload(opt_mod) - from quantui.optimizer import optimize_geometry + monkeypatch.setattr(opt_mod, "ASE_AVAILABLE", False) - with pytest.raises(ImportError, match="pip install"): - optimize_geometry(_h2()) + with pytest.raises(ImportError) as exc_info: + opt_mod.optimize_geometry(_h2()) + msg = str(exc_info.value) + assert "pip install" in msg or "conda install" in msg - def test_error_message_is_actionable(self, monkeypatch): - import quantui.ase_bridge as bridge - monkeypatch.setattr(bridge, "ASE_AVAILABLE", False) +class TestOptimizeGeometryPostHfGuard: + """M2 audit fix (2026-07-14): post-HF methods raise a clear ValueError. - import importlib + Regression: optimize_geometry() had no special-casing for MP2/CCSD/ + CCSD(T) — _QuantUIPySCFCalc.calculate() silently treated them as a DFT + xc functional (mf.xc = "CCSD"), failing deep inside PySCF with a cryptic + "LibXCFunctional: name 'CCSD' not found" instead of a clear message. + The guard fires before any PySCF import, so it needs ASE but not PySCF. + """ - import quantui.optimizer as opt_mod + ase_only = pytest.mark.skipif(not ASE_AVAILABLE, reason="ase not installed") - importlib.reload(opt_mod) + @ase_only + @pytest.mark.parametrize("method", ["MP2", "CCSD", "CCSD(T)"]) + def test_post_hf_method_raises_value_error(self, method): from quantui.optimizer import optimize_geometry - with pytest.raises(ImportError) as exc_info: - optimize_geometry(_h2()) - msg = str(exc_info.value) - assert "pip install" in msg or "conda install" in msg + with pytest.raises(ValueError, match="post-HF"): + optimize_geometry(_h2(), method=method, basis="STO-3G") # ============================================================================ diff --git a/tests/test_pes_scan.py b/tests/test_pes_scan.py index 069544f..585ae55 100644 --- a/tests/test_pes_scan.py +++ b/tests/test_pes_scan.py @@ -123,6 +123,63 @@ def test_empty_energies_returns_nan(self): assert r.energies_relative_kcal == [] +class TestPESScanResultNanHandling: + """M6 audit fix (2026-07-14): NaN from failed scan points must not + poison min()/max() based on ordering. + + Regression: Python's min()/max() are order-dependent when NaN is + present — a NaN as the first element "wins" (every comparison + against it is False) and poisons the result; a NaN later in the list + is correctly skipped. energy_hartree, energies_relative_kcal, and + summary() all called min()/max() directly on energies_hartree + without filtering, so whether a scan came out usable depended on + *which* point happened to fail. + """ + + def _make(self, energies, scan_type="bond"): + from quantui.pes_scan import PESScanResult + + n = len(energies) + mol = _h2() + return PESScanResult( + formula="H2", + method="RHF", + basis="STO-3G", + scan_type=scan_type, + atom_indices=[0, 1], + scan_parameter_values=[0.5 + i * 0.3 for i in range(n)], + energies_hartree=list(energies), + coordinates_list=[mol] * n, + converged_all=False, + ) + + def test_energy_hartree_ignores_leading_nan(self): + r = self._make([float("nan"), -1.10, -1.117, -1.05]) + assert r.energy_hartree == pytest.approx(-1.117) + + def test_energy_hartree_ignores_middle_nan(self): + r = self._make([-1.10, -1.117, float("nan"), -1.05]) + assert r.energy_hartree == pytest.approx(-1.117) + + def test_energy_hartree_ignores_trailing_nan(self): + r = self._make([-1.10, -1.117, -1.05, float("nan")]) + assert r.energy_hartree == pytest.approx(-1.117) + + def test_energies_relative_kcal_not_poisoned_by_leading_nan(self): + r = self._make([float("nan"), -1.10, -1.117, -1.05]) + rel = r.energies_relative_kcal + import math + + assert math.isnan(rel[0]) + assert min(v for v in rel[1:]) == pytest.approx(0.0, abs=1e-9) + + def test_summary_has_no_nan_with_leading_failed_point(self): + r = self._make([float("nan"), -1.10, -1.117, -1.05]) + summary = r.summary() + assert "nan" not in summary.lower() + assert "-1.11700000 Ha" in summary + + # ── run_pes_scan validation (no PySCF needed) ───────────────────────────────── @@ -159,6 +216,20 @@ def test_steps_less_than_2_raises(self): with pytest.raises((ImportError, ValueError)): run_pes_scan(_h2(), scan_type="bond", atom_indices=[0, 1], steps=1) + @pytest.mark.parametrize("method", ["MP2", "CCSD", "CCSD(T)"]) + def test_post_hf_method_raises(self, method): + """M2 audit fix (2026-07-14): post-HF methods raise a clear error. + + Regression: run_pes_scan() had no special-casing for MP2/CCSD/ + CCSD(T) — _QuantUIPySCFCalc.calculate() (shared with optimizer.py) + silently treated them as a DFT xc functional, failing deep inside + PySCF with a cryptic "LibXCFunctional" error instead of a clear one. + """ + from quantui.pes_scan import run_pes_scan + + with pytest.raises((ImportError, ValueError), match="post-HF|ASE"): + run_pes_scan(_h2(), method=method, scan_type="bond", atom_indices=[0, 1]) + # ── App widget integration ──────────────────────────────────────────────────── @@ -327,3 +398,113 @@ def test_h2_bond_scan_minimum_near_equilibrium(self): min_idx = e_rel.index(min(e_rel)) min_val = result.scan_parameter_values[min_idx] assert 0.5 <= min_val <= 1.5 # broad tolerance for 6-step coarse scan + + def test_failed_scan_point_falls_back_to_last_good_geometry(self, monkeypatch): + """M6 audit fix (2026-07-14): a failed point's frame must be the + last successful geometry, not always the original input molecule. + + Forces the 2nd scan point's BFGS relaxation to raise, then checks + that point's coordinates match point 1's (the last good frame) — + not the raw input geometry, which is what the bug used to record + regardless of how far the scan had already progressed. + """ + import math + + import ase.optimize as ase_optimize + + from quantui.pes_scan import run_pes_scan + + real_bfgs = ase_optimize.BFGS + call_count = {"n": 0} + + class _FlakyBFGS: + def __init__(self, *args, **kwargs): + call_count["n"] += 1 + self._real = real_bfgs(*args, **kwargs) + + def run(self, *args, **kwargs): + if call_count["n"] == 2: + raise RuntimeError("simulated optimizer failure") + return self._real.run(*args, **kwargs) + + monkeypatch.setattr(ase_optimize, "BFGS", _FlakyBFGS) + + result = run_pes_scan( + _water(), + method="RHF", + basis="STO-3G", + scan_type="bond", + atom_indices=[0, 1], # O-H bond; 3 atoms so BFGS actually runs + start=0.85, + stop=1.15, + steps=3, + ) + + assert result.converged_all is False + assert math.isnan(result.energies_hartree[1]) + failed_frame = result.coordinates_list[1] + last_good_frame = result.coordinates_list[0] + input_frame = _water() + assert failed_frame.coordinates == last_good_frame.coordinates + assert failed_frame.coordinates != input_frame.coordinates + + +@_pyscf_available +@pytest.mark.slow +class TestRunPesScanAngleDihedral: + """M7 audit fix (2026-07-14): angle/dihedral scans must actually run. + + Regression: FixInternals(angles=[...]) / FixInternals(dihedrals=[...]) + (the radian-based, deprecated kwargs) don't just emit a FutureWarning + with the ASE version QuantUI targets (verified against 3.29.0) — they + raise "setting an array element with a sequence" from an internal + np.asarray reshape, unconditionally, for every angle/dihedral + constraint. Every angle and dihedral PES scan failed at 100% of its + points as a result (silently, caught by the per-point try/except). + Switching to angles_deg/dihedrals_deg (plain degrees, no radian + conversion) fixes this. + """ + + def test_angle_scan_converges_without_nan(self): + from quantui.pes_scan import run_pes_scan + + result = run_pes_scan( + _water(), + method="RHF", + basis="STO-3G", + scan_type="angle", + atom_indices=[1, 0, 2], # H-O-H angle + start=95.0, + stop=115.0, + steps=3, + ) + assert result.converged_all is True + assert all(e == e for e in result.energies_hartree) # no NaN + + def test_dihedral_scan_converges_without_nan(self): + from quantui.molecule import Molecule + from quantui.pes_scan import run_pes_scan + + # Hydrogen peroxide (H-O-O-H) — smallest molecule with a real + # dihedral degree of freedom. + h2o2 = Molecule( + atoms=["H", "O", "O", "H"], + coordinates=[ + [0.9, 0.0, 0.9], + [0.0, 0.0, 0.75], + [0.0, 0.0, -0.75], + [-0.9, 0.3, -0.9], + ], + ) + result = run_pes_scan( + h2o2, + method="RHF", + basis="STO-3G", + scan_type="dihedral", + atom_indices=[0, 1, 2, 3], + start=60.0, + stop=120.0, + steps=3, + ) + assert result.converged_all is True + assert all(e == e for e in result.energies_hartree) # no NaN diff --git a/tests/test_progress.py b/tests/test_progress.py index 1446362..dd65456 100644 --- a/tests/test_progress.py +++ b/tests/test_progress.py @@ -125,5 +125,18 @@ def test_empty_message_ok(self): def test_html_special_chars_in_message(self): sp = StepProgress(["A"]) sp.complete(0, "Found <3 atoms & 2 bonds") - # Should not crash — message included as-is - assert sp.widget.value + # L audit fix: labels/messages must be HTML-escaped, not embedded + # as-is — a message containing "<"/">" (e.g. echoing a parse error + # back from user input) used to break the widget's markup instead + # of rendering as literal text. + assert "Found <3 atoms & 2 bonds" in sp.widget.value + assert "<3 atoms" not in sp.widget.value + + def test_html_escapes_script_tag_in_label_and_message(self): + sp = StepProgress(["Parse "]) + sp.fail(0, "") + value = sp.widget.value + assert "