diff --git a/pyproject.toml b/pyproject.toml index 966c18f..1519738 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -82,8 +82,19 @@ Issues = "https://github.com/The-Schultz-Lab/QuantUI/issues" # the terminal (``quantui log tail -n 50``, etc.). See ``quantui/cli.py``. quantui = "quantui.cli:main" -[tool.setuptools] -packages = ["quantui", "quantui.backends"] +[tool.setuptools.packages.find] +# AUDIT F20 — this used to be an explicit `packages = ["quantui", +# "quantui.backends"]` list under [tool.setuptools]. quantui.engines is a +# real subpackage (quantui/engines/__init__.py + base.py/pyscf_engine.py/ +# pyfock_engine.py) that was never added to that list: the built wheel +# silently omitted it entirely (ModuleNotFoundError importing +# quantui.engines from an installed, non-editable wheel), while every CI +# run stayed green because CI installs editable (`pip install -e`), which +# doesn't go through this list at all. Automatic discovery means the next +# new subpackage ships by construction instead of needing a second edit +# here that's easy to forget. +include = ["quantui", "quantui.*"] +exclude = ["tests", "tests.*"] [tool.setuptools.package-data] # Bundled molecule library (M-STRUCT): the indexed SQLite store + the @@ -176,6 +187,13 @@ dev = [ "pytest-cov>=4.0.0", "pytest-mock>=3.10.0", "pytest-xdist>=3.0.0", # parallel test execution (-n=auto in addopts) + # AUDIT F20 — the PEP 517 build frontend, used by + # tests/test_packaging.py to build a real wheel and inspect its + # contents (import/package-data smoke test against a *built* wheel, + # not the editable install every other test runs against — the gap + # that let quantui.engines silently vanish from the shipped package + # while CI, always editable, stayed green). + "build>=1.0.0", "mypy>=1.10,<2.4", # pinned, not an open floor — same reasoning as # black/ruff below: it must agree with .pre-commit-config.yaml's rev, # which is what CI enforces (M-TYPECHECK TYPE.2). A newer mypy silently diff --git a/quantui/app.py b/quantui/app.py index 46d5603..311cf94 100644 --- a/quantui/app.py +++ b/quantui/app.py @@ -5966,6 +5966,15 @@ def _run_required_final_single_point(target_mol, reason: str): method=self.method_dd.value, basis=self.basis_dd.value, progress_stream=log, # type: ignore[arg-type] + # AUDIT F16 — the app opens _ckpt and resolves _resume + # above (used by Geometry Opt / PES Scan / Reorganization + # Energy), but never passed either into interactive + # Frequency runs, so expensive IR/Raman displacement + # checkpointing (which freq_calc.py supports and the + # batch route already uses) could never be banked or + # resumed here. + checkpoint=_ckpt, + resume=_resume, ) result_html = self._format_freq_result(result) _displacements_serialized = None @@ -5978,6 +5987,25 @@ def _run_required_final_single_point(target_mol, reason: str): ).tolist() except Exception: pass + _thermo = getattr(result, "thermo", None) + _thermo_serialized = ( + { + "zpve_hartree": _thermo.zpve_hartree, + "H_hartree": _thermo.H_hartree, + "S_jmol": _thermo.S_jmol, + "G_hartree": _thermo.G_hartree, + "temperature_k": _thermo.temperature_k, + # AUDIT F18 — the temperature was already tracked on + # ThermoData; pressure and the model itself were not, + # and neither survived a save. Both are fixed by the + # harmonic-oscillator/rigid-rotor/ideal-gas model at + # 1 atm used throughout freq_calc.py's thermo block. + "pressure_atm": 1.0, + "approximation": "ideal_gas_rigid_rotor_harmonic_oscillator", + } + if _thermo is not None + else None + ) save_spectra = { "ir": { "frequencies_cm1": result.frequencies_cm1, @@ -5985,6 +6013,9 @@ def _run_required_final_single_point(target_mol, reason: str): "raman_activities": result.raman_activities, "zpve_hartree": result.zpve_hartree, "displacements": _displacements_serialized, + # AUDIT F18 — thermo (H, S, G) was computed and shown + # live but never made it into the saved result.json. + "thermo": _thermo_serialized, }, "molecule": { "atoms": list(calc_mol.atoms), @@ -6111,6 +6142,19 @@ def _run_required_final_single_point(target_mol, reason: str): str(k): v for k, v in result.chemical_shifts_ppm.items() }, "reference_compound": result.reference_compound, + # AUDIT F17 — the backend records which reference + # shielding constants were actually used and + # whether that was a fallback substitution (e.g. + # a method/basis combo with no matching reference, + # falling back to a different level of theory's + # constants); the local save used to drop both, + # losing that calibration provenance on replay. + # The batch NMR serializer (nmr_result_payload) + # already includes them. + "reference_key": getattr(result, "reference_key", ""), + "is_fallback_reference": getattr( + result, "is_fallback_reference", False + ), } } save_type = "nmr" diff --git a/quantui/app_analysis.py b/quantui/app_analysis.py index 354e0ab..56ce282 100644 --- a/quantui/app_analysis.py +++ b/quantui/app_analysis.py @@ -233,6 +233,7 @@ def apply_analysis_context(app: Any, ctx: Any) -> None: app._last_orb_info = None app._last_orb_mo_coeff = None app._last_orb_mo_occ = None + app._last_orb_method = None # Mulliken state consumed by the Populations panel — reset so a context # without charges cannot leak the prior calc's chart into this one. app._last_mulliken_symbols = None diff --git a/quantui/app_formatters.py b/quantui/app_formatters.py index 97b1525..de4602b 100644 --- a/quantui/app_formatters.py +++ b/quantui/app_formatters.py @@ -31,6 +31,28 @@ def _method_basis_label(method: str, basis: str, scf_variant: str | None) -> str return label +# Open-shell SCF variants whose orbitals split into separate alpha/beta +# channels (UHF/UKS: two independent sets; ROHF: a single spatial-orbital +# set but singly-occupied orbitals that still only have a well-defined +# alpha-channel HOMO/LUMO in the usual sense). +_OPEN_SHELL_SCF_VARIANTS = frozenset({"UHF", "UKS", "ROHF", "ROKS"}) + + +def _homo_lumo_gap_label(scf_variant: str | None) -> str: + """ "HOMO-LUMO gap", qualified "(α)" for an open-shell reference. + + AUDIT additional-concerns — session_calc.py's gap extraction always + reads spin channel 0 (alpha) for a 2-D ``mo_energy`` array (its own + comment: "UHF: ... use alpha spin for the gap estimate"). The result + card never said so, presenting a single-channel number as if it were + an unqualified property — silently dropping the beta-channel gap, + which can differ meaningfully for an open-shell system. + """ + if scf_variant and scf_variant.upper() in _OPEN_SHELL_SCF_VARIANTS: + return "HOMO-LUMO gap (α)" + return "HOMO-LUMO gap" + + def _result_card_open(*, accent: str | None = None, extra_style: str = "") -> str: border = accent or _theme.css.ACCENT_SUCCESS_ALT style = ( @@ -120,6 +142,22 @@ def _num(label: str, value: str) -> str: "(approximate 2-electron integrals)" ) + # AUDIT additional-concerns — session_calc.py extracts both properties + # from ``mf`` (the HF/DFT reference) even for MP2/CCSD/CCSD(T), which + # never builds a correlated density here. That is a real, potentially + # intentional teaching simplification (a correlated dipole/population + # needs a relaxed/unrelaxed density from the post-HF method itself, + # which this app does not compute) — but the card must say so instead + # of presenting an HF-level dipole/population as if it came from the + # requested correlated method. + _is_post_hf = _mp2 is not None or _ccsd is not None + _post_hf_note = ( + f' ' + "(HF reference — not a correlated MP2/CCSD property)" + if _is_post_hf + else "" + ) + _dip = get("dipole_moment_debye") if _dip is not None: _vec = get("dipole_vector_debye") @@ -136,7 +174,7 @@ def _num(label: str, value: str) -> str: f' ' "(magnitude only — μ components not saved)" ) - rows += _num("Dipole moment", _dip_str) + rows += _num("Dipole moment", _dip_str + _post_hf_note) _chg = get("mulliken_charges") _syms = get("atom_symbols") @@ -146,13 +184,19 @@ def _num(label: str, value: str) -> str: f'' f"Mulliken charges" f'{_charge_str}' + f'word-break:break-all">{_charge_str}{_post_hf_note}' ) return rows def format_result(r: Any) -> str: """Format a single-point-style result card.""" + # AUDIT F07 — for CCSD/CCSD(T), r.converged also folds in the CC + # amplitude solve's own convergence, so a bare "SCF converged" label + # would misleadingly blame the reference SCF for a CC-only failure. + _conv_label = ( + "Converged" if getattr(r, "cc_converged", None) is not None else "SCF converged" + ) _conv = "Yes" if r.converged else "No (treat results with caution)" _cc = _converged_color(r.converged) _gap = f"{r.homo_lumo_gap_ev:.4f} eV" if r.homo_lumo_gap_ev is not None else "N/A" @@ -167,8 +211,12 @@ def format_result(r: Any) -> str: f"{r.energy_hartree:.8f} Ha  ({r.energy_ev:.4f} eV)", _theme.css.TEXT_HEADING, ), - ("HOMO-LUMO gap", _gap, _theme.css.TEXT_HEADING), - ("SCF converged", _conv, _cc), + ( + _homo_lumo_gap_label(getattr(r, "scf_variant", None)), + _gap, + _theme.css.TEXT_HEADING, + ), + (_conv_label, _conv, _cc), ( "SCF iterations", ( @@ -222,6 +270,9 @@ def format_opt_result(r: Any) -> str: def format_freq_result(r: Any) -> str: """Format a frequency-analysis result card.""" + # AUDIT F15 — r.converged now also requires the Hessian/harmonic- + # analysis step to have completed, not just the reference SCF, so the + # row is labeled/colored on overall status rather than "SCF converged". _conv = "Yes" if r.converged else "No (treat with caution)" _cc = _converged_color(r.converged) n_real = r.n_real_modes() @@ -239,7 +290,7 @@ def format_freq_result(r: Any) -> str: _rows = ( f'SCF energy' f'{r.energy_hartree:.8f} Ha' - f'SCF converged' + f'Converged' f'{_conv}' f'Real modes' f'{n_real}' @@ -282,15 +333,21 @@ def format_freq_result(r: Any) -> str: def format_tddft_result(r: Any) -> str: """Format a TD-DFT / UV-Vis result card.""" + # AUDIT F08 — r.converged now folds in per-root TD convergence, so the + # row is labeled/colored on overall status, not just the ground SCF. _conv = "Yes" if r.converged else "No (treat with caution)" _cc = _converged_color(r.converged) + _n_converged = getattr(r, "n_converged_states", None) + _states_detail = str(len(r.excitation_energies_ev)) + if _n_converged is not None and _n_converged != len(r.excitation_energies_ev): + _states_detail += f" ({_n_converged} converged)" header_rows = ( f'Ground-state energy' f'{r.energy_hartree:.8f} Ha' - f'SCF converged' + f'Converged' f'{_conv}' f'States computed' - f'{len(r.excitation_energies_ev)}' + f'{_states_detail}' ) exc_table = "" if r.excitation_energies_ev: @@ -342,6 +399,16 @@ def format_nmr_result(r: Any) -> str: f'Reference' f'{r.reference_compound} ({r.method}/{r.basis})' ) + # AUDIT F17 — surface a fallback-reference substitution explicitly + # rather than letting the row above imply an exact-match reference. + if getattr(r, "is_fallback_reference", False): + _ref_key = getattr(r, "reference_key", "") or "a different level of theory" + header_rows += ( + f'' + f'' + f"⚠ No reference at {r.method}/{r.basis} — shifts use {_ref_key} " + "constants instead." + ) def _nmr_table(label: str, shifts: list, sym: str) -> str: if not shifts: @@ -672,6 +739,23 @@ def _label(kind: str) -> str: f'' f"{body}
" ) + # AUDIT additional-concerns — "Ion state" above reports a multiplicity + # chosen by electron-count parity/minimal spin (see + # reorganization_energy._ion_multiplicity's own docstring), which is a + # convenient default, NOT a ground-state determination — most relevant + # for a transition-metal ion, where the true ground state can be + # higher-spin. Any PCM solvent selected applies only to the four + # single-point energies above, evaluated at gas-phase-optimized + # geometries (the relaxations themselves are not solvent-optimized). + # Both scope notes are stated here, once, rather than left implicit. + blocks.append( + f'
' + "Ion multiplicity is the minimal-spin default from electron-count " + "parity, not a ground-state determination (relevant for transition-" + "metal ions). Solvent (if selected) applies only to the single-" + "point energies above; geometries are optimized in the gas phase." + "
" + ) return "".join(blocks) @@ -788,7 +872,11 @@ def format_past_result(data: dict[str, Any], result_dir: Optional[Path] = None) f"{data['energy_hartree']:.8f} Ha  ({data['energy_ev']:.4f} eV)", _theme.css.TEXT_HEADING, ), - ("HOMO-LUMO gap", _gap, _theme.css.TEXT_HEADING), + ( + _homo_lumo_gap_label(data.get("scf_variant")), + _gap, + _theme.css.TEXT_HEADING, + ), ("SCF converged", _conv, _cc), ( "SCF iterations", @@ -819,6 +907,32 @@ def format_past_result(data: dict[str, Any], result_dir: Optional[Path] = None) f'border:1px solid {_theme.css.BORDER}" width="173" height="108" />' ) + # AUDIT F18 — restore persisted frequency thermochemistry into the + # History card. Older saved results (or a Hessian-only run with no + # thermo block) have no "thermo" key at all; that's a silent no-op, + # not an error. + _thermo_html = "" + if ct == "frequency": + _thermo = ((data.get("spectra") or {}).get("ir") or {}).get("thermo") + if _thermo: + _kj = 2625.5 # kJ/mol per Hartree + _thermo_html = ( + f'' + f"— Thermochemistry at {_thermo.get('temperature_k', 298.15):.0f} K" + f" / {_thermo.get('pressure_atm', 1.0):.0f} atm —" + f"" + f'ZPVE' + f'{_thermo["zpve_hartree"]:.6f} Ha' + f'H' + f'{_thermo["H_hartree"]:.6f} Ha' + f'S' + f'{_thermo["S_jmol"]:.2f} J/(mol·K)' + f'G' + f'{_thermo["G_hartree"]:.6f} Ha' + f" ({_thermo['G_hartree'] * _kj:.2f} kJ/mol)" + ) + # Reorganization-energy channels (REORG.1). This is the reported bug: the # card came back without the numbers the calculation exists to produce. # Keyed on the calc type AND the payload, so a reorg result saved before λ @@ -840,6 +954,6 @@ def format_past_result(data: dict[str, Any], result_dir: Optional[Path] = None) f'{_method_basis_label(data["method"], data["basis"], data.get("scf_variant"))}' f' {ts}' + _result_card_table_open() - + f"{_rows}{_extra}{_reorg_html}" + + f"{_rows}{_extra}{_thermo_html}{_reorg_html}" + _RESULT_CARD_CLOSE ) diff --git a/quantui/app_runflow.py b/quantui/app_runflow.py index eb66235..b86c824 100644 --- a/quantui/app_runflow.py +++ b/quantui/app_runflow.py @@ -404,10 +404,55 @@ def on_basis_fix(app: Any, btn: Any = None) -> None: pass +# AUDIT F11 — run_freq_calc, run_tddft_calc, run_nmr_calc, and run_pes_scan +# don't accept a solvent argument at all, so checking "Implicit solvent +# (PCM)" for Frequency/UV-Vis/NMR Shielding/PES Scan used to be a complete +# no-op: the backend received molecule/method/basis/progress_stream and +# silently ran gas-phase while the UI kept showing the box checked. +# +# "Single Point" has full PCM support (session_calc.run_in_session). +# "Geometry Opt" and "Reorganization Energy" get a real, but partial, +# solvent treatment: the geometry optimization itself runs gas-phase, then +# a single point WITH solvent is computed at the final geometry and its +# energy/orbitals replace the last trajectory frame's — a real published +# approximation, not silently ignored, but not a solvated optimization +# either (see _run_required_final_single_point in app.py and +# reorganization_energy.py's own docstring). +_SOLVENT_SUPPORTED_CALC_TYPES = frozenset( + {"Single Point", "Geometry Opt", "Reorganization Energy"} +) + + +def _update_solvent_control_for_calc_type(app: Any, ct: str) -> None: + """Disable the solvent checkbox for calc types that would silently + ignore it, so a checked box can never mean "no effect" (AUDIT F11).""" + try: + cb = app.solvent_cb + except AttributeError: + return + if ct in _SOLVENT_SUPPORTED_CALC_TYPES: + cb.disabled = False + if ct in ("Geometry Opt", "Reorganization Energy"): + # AUDIT F11 — label the approximation explicitly rather than + # letting a checked box imply a fully solvated optimization. + cb.description = ( + "Implicit solvent (PCM) — gas-phase optimization, " + "solvated final single point" + ) + else: + cb.description = "Implicit solvent (PCM)" + else: + cb.value = False # also hides solvent_dd via on_solvent_cb_changed + cb.disabled = True + cb.description = f"Implicit solvent (PCM) — not supported for {ct}" + + def on_calc_type_changed(app: Any, change: Any, *, layout_fn: Any) -> None: """Update extra options panel based on selected calculation type.""" ct = change["new"] + _update_solvent_control_for_calc_type(app, ct) + from quantui.freq_calc import is_freq_mode_seed if ( @@ -481,7 +526,7 @@ def on_calc_type_changed(app: Any, change: Any, *, layout_fn: Any) -> None: app._tddft_seed_note, widgets.HTML( f'⚠ Requires a DFT ' - "functional (e.g. B3LYP, PBE0). RHF/UHF will run TDHF (CIS) " + "functional (e.g. B3LYP, PBE0). RHF/UHF will run TDHF/RPA " "instead." ), ] @@ -2285,11 +2330,30 @@ def checkpoint_identity(app: Any) -> Any: molecule = getattr(app, "_molecule", None) if molecule is None: return None + _ct = calc_type_key(app) + extra: tuple = () + if _ct == "pes_scan": + # AUDIT F10 — calc_type="pes_scan" alone doesn't distinguish a + # bond scan from an angle/dihedral scan, or which atoms it + # scans, so two different scan configurations of the same + # molecule/method/basis used to collide on the same + # resume_key. + try: + extra = ( + str(app._scan_type_dd.value), + str(app._scan_atom1.value), + str(app._scan_atom2.value), + str(app._scan_atom3.value), + str(app._scan_atom4.value), + ) + except Exception: # noqa: BLE001 — checkpointing is never load-bearing + extra = () return CalcIdentity.from_molecule( molecule, - calc_type=calc_type_key(app), + calc_type=_ct, method=app.method_dd.value, basis=app.basis_dd.value, + extra=extra, ) except Exception: # noqa: BLE001 — checkpointing is never load-bearing return None diff --git a/quantui/app_visualization.py b/quantui/app_visualization.py index dc4b36f..dfaaf7d 100644 --- a/quantui/app_visualization.py +++ b/quantui/app_visualization.py @@ -1046,6 +1046,10 @@ def update_uv_vis_figure(app: Any, mode: str, fwhm: float) -> None: n_points = max(600, int((x_max - x_min) * 2.0)) x_grid = _np.linspace(x_min, x_max, n_points) y_grid = _np.zeros_like(x_grid) + # AUDIT additional-concerns — same height-normalized Lorentzian + # convention as ir_plot.py/raman_plot.py (see ir_plot.py's + # module docstring): peak height = supplied oscillator + # strength, area scales with FWHM. Deliberate, not a bug. for x0, amp in zip(wl, osc): y_grid += amp * (gamma**2 / ((x_grid - x0) ** 2 + gamma**2)) fig.add_trace( @@ -1293,6 +1297,15 @@ def show_orbital_diagram(app: Any, result: Any) -> bool: 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) + # AUDIT additional-concerns — snapshot the method that actually + # produced this mo_coeff, from the result object itself, rather than + # reading the live Method dropdown at cube-generation time (below). + # The dropdown can change (or a different History result can be + # loaded) between this call and Generate being pressed, at which point + # the dropdown no longer describes the orbitals actually being + # exported — the cube's own provenance comment used to silently name + # whatever method the dropdown showed at that later moment instead. + app._last_orb_method = str(getattr(result, "method", "") or "") plotly_rendered = False try: @@ -1595,8 +1608,11 @@ def _show_range_err() -> None: # 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. - _charge, _spin = infer_charge_and_spin(mol_atom, mo_occ_for_charge) + # (H3O+, OH-, radicals, ...) don't fail to build in PySCF. Passing + # mol_basis lets AUDIT F14's ECP-aware charge inference apply. + _charge, _spin = infer_charge_and_spin( + mol_atom, mo_occ_for_charge, basis=mol_basis + ) # ORBX.2: the user-chosen cubegen grid. Read at generate time rather # than cached, so changing the dropdown affects the next Generate @@ -1614,13 +1630,17 @@ def _show_range_err() -> None: _grid = ISO_RESOLUTION_PRESETS.get( _res_key, ISO_RESOLUTION_PRESETS[DEFAULT_ISO_RESOLUTION] ) - # M-EXPORT2 EXP2.4 / M-ORBEXPORT ORBX.4: best-effort provenance, not a - # re-verified guarantee — the live method dropdown, not necessarily - # what actually produced the stored mo_coeff (e.g. after a History - # replay of a differently-computed result). - _method_for_provenance = str( - getattr(getattr(app, "method_dd", None), "value", "") or "" - ) + # M-EXPORT2 EXP2.4 / M-ORBEXPORT ORBX.4 — AUDIT additional-concerns: + # this used to read the LIVE method dropdown, which is not + # necessarily what actually produced the stored mo_coeff (e.g. + # after a History replay of a differently-computed result, or if + # the dropdown is changed between loading the orbitals and + # pressing Generate). ``_last_orb_method`` is snapshotted from the + # result object itself at the moment its orbitals were loaded + # (show_orbital_diagram), so it stays correct regardless of what + # the dropdown shows later — immutable result provenance instead + # of a mutable, disconnectable UI control. + _method_for_provenance = str(getattr(app, "_last_orb_method", "") or "") generate_cube_from_arrays( mol_atom, mol_basis, diff --git a/quantui/backends/cluster_config.py b/quantui/backends/cluster_config.py index 19327c9..17ceaeb 100644 --- a/quantui/backends/cluster_config.py +++ b/quantui/backends/cluster_config.py @@ -139,8 +139,8 @@ def default_staging_root() -> Path: #SBATCH --ntasks={cores} #SBATCH --mem={memory}G #SBATCH --time={walltime} -#SBATCH --output={output_file} -#SBATCH --error={error_file}{optional_directives} +#SBATCH --output="{output_file}" +#SBATCH --error="{error_file}"{optional_directives} set -euo pipefail diff --git a/quantui/backends/slurm.py b/quantui/backends/slurm.py index 7f4491c..5333745 100644 --- a/quantui/backends/slurm.py +++ b/quantui/backends/slurm.py @@ -9,6 +9,7 @@ import json import logging +import shlex import subprocess import sys import time @@ -169,13 +170,23 @@ def dispatch( return request.request_id def _worker_command(self, request_path: Path, staging_dir: Path) -> str: - py = sys.executable - inner = f"{py} -m quantui.backends.worker --request {request_path}" + # AUDIT F21 — every path here is inserted into shell text (this + # string is embedded verbatim into the generated sbatch script), + # not passed as an argv list, so an unquoted path containing a + # space (or any other shell metacharacter) splits into multiple + # arguments. E.g. unquoted "/tmp/audit folder/request.json" became + # "--request /tmp/audit" plus a stray "folder/request.json" + # argument. shlex.quote makes every one of these shell-safe + # regardless of what it contains. + py = shlex.quote(sys.executable) + request_arg = shlex.quote(str(request_path)) + inner = f"{py} -m quantui.backends.worker --request {request_arg}" if self.use_apptainer: - image = self.apptainer_image + image = shlex.quote(self.apptainer_image) + staging_arg = shlex.quote(str(staging_dir)) return ( - f'apptainer exec --nv --bind "$HOME:$HOME" --pwd "{staging_dir}" ' - f'"{image}" {inner}' + f'apptainer exec --nv --bind "$HOME:$HOME" --pwd {staging_arg} ' + f"{image} {inner}" ) return inner diff --git a/quantui/backends/slurm_ingest.py b/quantui/backends/slurm_ingest.py index dc0103b..5d59106 100644 --- a/quantui/backends/slurm_ingest.py +++ b/quantui/backends/slurm_ingest.py @@ -43,6 +43,24 @@ def _finalize_history_entry(saved_dir: Path) -> None: def _basic_result(payload: dict[str, Any], record: JobRecord) -> SimpleNamespace: + """Reconstruct a result-alike object from a worker's staging JSON. + + AUDIT F12 — this used to reconstruct only 7 generic fields even though + ``session_result_payload()`` (and the other ``*_result_payload`` + builders) already serialize Mulliken charges, dipole, atom symbols, + SCF variant/rescue provenance, post-HF correlation breakdown, and + solvent/GPU/density-fit metadata. ``save_result()`` reads every one of + these via ``getattr(result, ..., default)``, so silently omitting them + here made ``save_result`` write them as null regardless of whether the + JSON actually had real values — a real water round trip lost the + dipole, charges, atom symbols, and RHF provenance. + + Every field below is read defensively (``.get`` with no required key) + because not every calc type's payload builder sets every field — + absent ones round-trip as the same ``None``/default ``save_result`` + already treats as "not applicable for this calc type", matching its + own documented contract. + """ return SimpleNamespace( energy_hartree=float(payload.get("energy_hartree", float("nan"))), homo_lumo_gap_ev=payload.get("homo_lumo_gap_ev"), @@ -51,6 +69,21 @@ def _basic_result(payload: dict[str, Any], record: JobRecord) -> SimpleNamespace method=str(payload.get("method", record.request_obj.method)), basis=str(payload.get("basis", record.request_obj.basis)), formula=str(payload.get("formula", "?")), + mulliken_charges=payload.get("mulliken_charges"), + dipole_moment_debye=payload.get("dipole_moment_debye"), + dipole_vector_debye=payload.get("dipole_vector_debye"), + atom_symbols=payload.get("atom_symbols"), + scf_rescue_stage=payload.get("scf_rescue_stage", "none"), + scf_variant=payload.get("scf_variant") or None, + mp2_correlation_hartree=payload.get("mp2_correlation_hartree"), + ccsd_correlation_hartree=payload.get("ccsd_correlation_hartree"), + ccsd_t_correction_hartree=payload.get("ccsd_t_correction_hartree"), + cc_converged=payload.get("cc_converged"), + dispersion_applied=payload.get("dispersion_applied"), + solvent=payload.get("solvent"), + gpu_used=bool(payload.get("gpu_used", False)), + gpu_name=payload.get("gpu_name"), + density_fit=bool(payload.get("density_fit", False)), ) diff --git a/quantui/backends/worker.py b/quantui/backends/worker.py index 5325ae2..e59aca2 100644 --- a/quantui/backends/worker.py +++ b/quantui/backends/worker.py @@ -49,6 +49,17 @@ _SUPPORTED_CALC_TYPES = frozenset(CALC_TYPES) +# AUDIT F11 — the science APIs behind these calc types (run_freq_calc, +# run_tddft_calc, run_nmr_calc, run_pes_scan, optimize_geometry) don't +# accept a solvent argument at all, so request.solvent used to be silently +# dropped rather than either applied or rejected. Only these two runners +# actually thread request.solvent through to a PCM-capable API +# (session_calc.run_in_session / reorganization_energy.run_reorganization_ +# energy, both real gas+PCM or PCM-single-point implementations — see each +# runner below and reorganization_energy's own docstring for the +# gas-phase-optimization + PCM-single-point approximation it documents). +_SOLVENT_SUPPORTED_CALC_TYPES = frozenset({"single_point", "reorganization_energy"}) + def _write_progress( staging_dir: Path, stage: str, message: str, percent: float @@ -91,6 +102,7 @@ def _begin_worker_checkpoint( basis: str, staging_dir: Path, log_stream, + extra: tuple = (), ): """Open a checkpoint for this job, scoped to its own staging directory (M-CLUSTER2 CL2.8). @@ -112,7 +124,7 @@ def _begin_worker_checkpoint( from quantui.checkpoint import CalcIdentity, Checkpoint identity = CalcIdentity.from_molecule( - molecule, calc_type=calc_type, method=method, basis=basis + molecule, calc_type=calc_type, method=method, basis=basis, extra=extra ) ckpt = Checkpoint(identity, root=staging_dir / ".checkpoint") ckpt.attach_log(log_stream) @@ -487,6 +499,11 @@ def _run_pes_scan(request: CalculationRequest, staging_dir: Path, log_stream) -> basis=request.basis, staging_dir=staging_dir, log_stream=log_stream, + # AUDIT F10 — calc_type="pes_scan" alone can't tell a bond scan + # apart from an angle scan of atoms 2,3,4, so two different scan + # configurations of the same molecule/method/basis used to collide + # on the same resume_key. + extra=(scan_type, *(str(i) for i in atom_indices)), ) if resumable: n_points = len(ckpt.completed_points()) @@ -625,6 +642,27 @@ def run_worker_request(request_path: Path) -> CalculationResult: save_type=calc_type, ) + # AUDIT F11 — a solvent set for a calc_type whose science API doesn't + # accept one used to be silently dropped (gas-phase result, solvent + # label nowhere). Fail the request instead of running a calculation + # the user asked for solvated and got gas-phase. + if request.solvent and calc_type not in _SOLVENT_SUPPORTED_CALC_TYPES: + msg = ( + f"solvent={request.solvent!r} was requested for " + f"calc_type={calc_type!r}, which does not support PCM solvation " + f"(only {', '.join(sorted(_SOLVENT_SUPPORTED_CALC_TYPES))} do). " + "Resubmit without a solvent, or use a supported calc_type." + ) + _append_log(staging_dir, msg) + return _error_result( + request, + staging_dir, + code="UNSUPPORTED_CAPABILITY", + message=msg, + retryable=False, + save_type=calc_type, + ) + runners: dict[str, Callable[..., Any]] = { "single_point": _run_single_point, "geometry_opt": _run_geometry_opt, diff --git a/quantui/backends/worker_payload.py b/quantui/backends/worker_payload.py index 75635ea..855c906 100644 --- a/quantui/backends/worker_payload.py +++ b/quantui/backends/worker_payload.py @@ -84,6 +84,20 @@ def session_result_payload(result) -> Dict[str, Any]: # M-UX2 UXP2.10 — see the other *_result_payload functions' matching # field. "scf_variant": getattr(result, "scf_variant", "") or None, + # AUDIT F12 — these were computed onto SessionResult but never + # reached staging JSON at all (a serialization-layer gap distinct + # from _basic_result's ingest-layer one): post-HF correlation + # breakdown, solvent/GPU/density-fit provenance, and the AUDIT + # F04/F07 dispersion/CC-convergence flags. + "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), + "cc_converged": getattr(result, "cc_converged", None), + "dispersion_applied": getattr(result, "dispersion_applied", None), + "solvent": getattr(result, "solvent", None), + "gpu_used": bool(getattr(result, "gpu_used", False)), + "gpu_name": getattr(result, "gpu_name", None), + "density_fit": bool(getattr(result, "density_fit", False)), } @@ -115,6 +129,23 @@ def freq_result_payload(result, molecule) -> Dict[str, Any]: displacements = np.asarray(result.displacements).tolist() except Exception: displacements = None + _thermo = getattr(result, "thermo", None) + _thermo_payload = ( + { + "zpve_hartree": _thermo.zpve_hartree, + "H_hartree": _thermo.H_hartree, + "S_jmol": _thermo.S_jmol, + "G_hartree": _thermo.G_hartree, + "temperature_k": _thermo.temperature_k, + # AUDIT F18 — pressure and the thermo model itself were never + # recorded anywhere; both are fixed by the harmonic-oscillator/ + # rigid-rotor/ideal-gas model at 1 atm used in freq_calc.py. + "pressure_atm": 1.0, + "approximation": "ideal_gas_rigid_rotor_harmonic_oscillator", + } + if _thermo is not None + else None + ) return { "calc_type": "frequency", "energy_hartree": result.energy_hartree, @@ -126,6 +157,8 @@ def freq_result_payload(result, molecule) -> Dict[str, Any]: "formula": result.formula, # M-UX2 UXP2.10 — see session_result_payload's matching field. "scf_variant": getattr(result, "scf_variant", "") or None, + # AUDIT F12 — was never serialized, though FreqResult carries it. + "density_fit": bool(getattr(result, "density_fit", False)), "spectra": { "ir": { "frequencies_cm1": list(result.frequencies_cm1), @@ -133,6 +166,10 @@ def freq_result_payload(result, molecule) -> Dict[str, Any]: "raman_activities": list(getattr(result, "raman_activities", []) or []), "zpve_hartree": result.zpve_hartree, "displacements": displacements, + # AUDIT F18 — thermo (H, S, G) was computed by freq_calc.py + # but discarded here; the saved JSON had only frequencies, + # intensities, activities, displacements, and ZPVE. + "thermo": _thermo_payload, }, "molecule": { "atoms": list(molecule.atoms), @@ -157,6 +194,8 @@ def tddft_result_payload(result) -> Dict[str, Any]: "formula": result.formula, # M-UX2 UXP2.10 — see session_result_payload's matching field. "scf_variant": getattr(result, "scf_variant", "") or None, + # AUDIT F12 — was never serialized, though TDDFTResult carries it. + "density_fit": bool(getattr(result, "density_fit", False)), "spectra": { "uv_vis": { "excitation_energies_ev": list(result.excitation_energies_ev), @@ -179,6 +218,8 @@ def nmr_result_payload(result) -> Dict[str, Any]: "formula": result.formula, # M-UX2 UXP2.10 — see session_result_payload's matching field. "scf_variant": getattr(result, "scf_variant", "") or None, + # AUDIT F12 — was never serialized, though NMRResult carries it. + "density_fit": bool(getattr(result, "density_fit", False)), "spectra": { "nmr": { "atom_symbols": list(result.atom_symbols), diff --git a/quantui/calculator.py b/quantui/calculator.py index af9bf8b..14b7869 100644 --- a/quantui/calculator.py +++ b/quantui/calculator.py @@ -92,10 +92,39 @@ def generate_calculation_script(self, output_path: Path) -> str: formula = self.molecule.get_formula() job_name = f"{formula}_{self.method}_{self.basis}" + # AUDIT F06 — embed the resolved ECP mapping so a heavy-element + # system (e.g. NaH/LANL2DZ) reproduces the in-app electron count + # instead of silently running all-electron with only the basis set. + # + # ecp_for_basis() imports pyscf (to call gto.basis.load_ecp()), but + # generating this script is a platform-independent, PySCF-free + # operation by design (the module docstring: "students can + # download and run independently" — e.g. Windows without PySCF + # generating a script meant to run on a Linux/WSL machine or + # cluster). Windows CI caught this: a bare pyscf import here broke + # script generation itself wherever PySCF isn't installed, not + # just execution. Fall back to {} (this platform's own pre-F06 + # behavior — never worse than before the fix) when PySCF is + # unavailable locally; any machine with PySCF still gets the + # correct ECP mapping. + try: + from .inorganic_guards import ecp_for_basis + + ecp = ecp_for_basis(self.basis, self.molecule.atoms) + except ImportError: + logger.warning( + "PySCF not installed on this machine; exported script's " + "ECP mapping defaults to {} (all-electron). Heavy-element " + "systems (e.g. LANL2DZ/def2 on Na and heavier) may need " + "mol.ecp set manually before running the script." + ) + ecp = {} + script_content = config.PYSCF_SCRIPT_TEMPLATE.format( job_name=job_name, method=self.method, basis=self.basis, + ecp=repr(ecp), geometry=geometry, charge=self.molecule.charge, spin=spin, diff --git a/quantui/checkpoint.py b/quantui/checkpoint.py index 691722b..a4565d7 100644 --- a/quantui/checkpoint.py +++ b/quantui/checkpoint.py @@ -136,10 +136,29 @@ class CalcIdentity: multiplicity: int = 1 atom_symbols: tuple = () coords: tuple = () + # AUDIT F10 — calc-type-specific discriminators that don't fit the + # generic fields above but still make two runs genuinely different + # calculations. The only current use is PES scan's (scan_type, + # *atom_indices): "calc_type" alone was just "pes_scan" for every + # scan configuration, so a bond scan and an angle scan of the same + # starting molecule/method/basis produced the SAME resume_key — + # resuming one could silently reuse the other's cached points (see + # points.jsonl's own per-point scan_type/atom_indices check in + # pes_scan.py for the defense-in-depth layer under this one). Included + # in resume_key (an exact-match requirement) but deliberately excluded + # from warm_start_key: an SCF density is still a good initial guess + # across different scan configurations of the same system. + extra: tuple = () @classmethod def from_molecule( - cls, molecule: Any, *, calc_type: str, method: str, basis: str + cls, + molecule: Any, + *, + calc_type: str, + method: str, + basis: str, + extra: tuple = (), ) -> CalcIdentity: """Build an identity from a :class:`~quantui.molecule.Molecule`.""" coords = getattr(molecule, "coordinates", None) @@ -157,6 +176,7 @@ def from_molecule( multiplicity=int(getattr(molecule, "multiplicity", 1) or 1), atom_symbols=tuple(str(a) for a in (getattr(molecule, "atoms", []) or [])), coords=coord_rows, + extra=tuple(str(e) for e in extra), ) @property @@ -190,6 +210,7 @@ def resume_key(self) -> str: self.warm_start_key, self.calc_type, _coords_digest(self.coords), + ",".join(self.extra), ] return hashlib.sha256("|".join(parts).encode("utf-8")).hexdigest()[:16] diff --git a/quantui/config.py b/quantui/config.py index f12bde6..83a0510 100644 --- a/quantui/config.py +++ b/quantui/config.py @@ -101,11 +101,13 @@ }, "wB97X-D": { "type": "dft", - "label": "wB97X-D — Range-Separated Hybrid + D3 Dispersion", + "label": "wB97X-D — Range-Separated Hybrid + Built-in Dispersion", "description": ( - "Range-separated hybrid functional with empirical D3 dispersion correction. " - "Excellent for non-covalent interactions, charge-transfer excitations, " - "and systems where long-range exchange matters." + "Range-separated hybrid functional (Chai & Head-Gordon, 2008) with its " + "own empirical dispersion correction baked into the fit — not an " + "externally applied Grimme D3 correction. Excellent for non-covalent " + "interactions, charge-transfer excitations, and systems where " + "long-range exchange matters." ), "use_for": "Non-covalent interactions, excited states, large organic molecules.", }, @@ -639,6 +641,11 @@ def main(): {geometry} ''' mol.basis = '{basis}' + # AUDIT F06 — the resolved ECP mapping (e.g. {{'Na': 'LANL2DZ'}} for + # LANL2DZ/def2 heavy elements, {{}} for an all-electron basis). Without + # this a heavy-element system runs all-electron here — a different + # Hamiltonian than the in-app calculation, not just numerical noise. + mol.ecp = {ecp} mol.charge = {charge} mol.spin = {spin} mol.verbose = 4 # Detailed output @@ -658,23 +665,33 @@ def main(): try: method = '{method}' - # Display name → PySCF xc string + external D3 dispersion. Matches - # quantui/session_calc.py resolve_xc + maybe_apply_d3. Important - # for methods that PySCF doesn't accept directly (notably - # wB97X-D — on dftd3's black-list; PBE-D3 — D3 must be applied - # externally via pyscf.dftd3). + # Display name → PySCF xc string + external D3 dispersion where + # needed. Matches quantui/session_calc.py resolve_xc + maybe_apply_d3. + # wB97X-D maps to its full LibXC name (the actual Chai/Head-Gordon + # 2008 functional, built-in dispersion) because PySCF's short-alias + # parser black-lists 'wb97x-d'/'wb97x_d' as ambiguous; PBE-D3 needs + # Grimme D3 applied externally via pyscf.dftd3. _XC_ALIAS = {{ 'M06-L': 'm06l', - 'wB97X-D': 'wb97x', + 'wB97X-D': 'hyb_gga_xc_wb97x_d', 'CAM-B3LYP': 'camb3lyp', 'PBE-D3': 'pbe', }} - _NEEDS_D3 = {{'PBE-D3', 'wB97X-D'}} + _NEEDS_D3 = {{'PBE-D3'}} if method == 'RHF': mf = scf.RHF(mol) elif method == 'UHF': mf = scf.UHF(mol) + elif method in ('MP2', 'CCSD', 'CCSD(T)'): + # AUDIT F13 — post-HF methods used to fall into the DFT branch + # below, setting mf.xc = 'MP2'/'CCSD'/'CCSD(T)' directly and + # failing with "LibXCFunctional: name '...' not found". The + # reference is scf.RHF(mol) regardless of spin — a factory + # that dispatches to true RHF for closed-shell (mol.spin == 0) + # and to ROHF for open-shell — matching + # quantui/session_calc.py's post-HF reference dispatch exactly. + mf = scf.RHF(mol) else: # DFT: auto-select RKS/UKS based on spin mf = dft.RKS(mol) if mol.spin == 0 else dft.UKS(mol) @@ -691,13 +708,66 @@ def main(): energy = _run_scf_with_rescue(mf) - if mf.converged: + # AUDIT F13 — MP2/CCSD/CCSD(T) correlation, mirroring + # quantui/session_calc.py including its AUDIT F07 convergence + # gating: post-HF work only runs on a converged reference, and + # CCSD(T) triples only run if the CCSD amplitudes themselves + # converged, rather than reporting a correlation "correction" on + # top of a wrong or unconverged Hamiltonian. + mp2_correlation = None + ccsd_correlation = None + ccsd_t_correction = None + cc_converged = None + if method == 'MP2': + if mf.converged: + from pyscf import mp as _mp + _mp2 = _mp.MP2(mf) + _e_corr, _ = _mp2.kernel() + mp2_correlation = _e_corr + energy += _e_corr + else: + print("Skipping MP2 -- reference SCF did not converge.") + elif method in ('CCSD', 'CCSD(T)'): + if mf.converged: + from pyscf import cc as _cc + _ccsd = _cc.CCSD(mf) + _e_corr_ccsd, _, _ = _ccsd.kernel() + cc_converged = bool(_ccsd.converged) + ccsd_correlation = _e_corr_ccsd + energy += _e_corr_ccsd + if method == 'CCSD(T)': + if cc_converged: + _e_t = _ccsd.ccsd_t() + ccsd_t_correction = _e_t + energy += _e_t + else: + print( + "Skipping CCSD(T) triples -- CCSD amplitudes " + "did not converge." + ) + else: + print("Skipping CCSD -- reference SCF did not converge.") + + # Overall success requires the reference SCF, and (when a coupled- + # cluster method was requested) CCSD's own amplitude convergence. + _overall_converged = mf.converged and (cc_converged is not False) + + if _overall_converged: print() print("=" * 60) print("Calculation Results") print("=" * 60) print(f"SCF converged: Yes") print(f"Total energy: {{energy:.8f}} Ha") + if mp2_correlation is not None: + print(f" HF reference : {{energy - mp2_correlation:.8f}} Ha") + print(f" MP2 correlation : {{mp2_correlation:.8f}} Ha") + if ccsd_correlation is not None: + _hf_e = energy - ccsd_correlation - (ccsd_t_correction or 0.0) + print(f" HF reference : {{_hf_e:.8f}} Ha") + print(f" CCSD correlation : {{ccsd_correlation:.8f}} Ha") + if ccsd_t_correction is not None: + print(f" (T) triples : {{ccsd_t_correction:.8f}} Ha") mo_e = mf.mo_energy if not isinstance(mf.mo_energy, list) else mf.mo_energy[0] mo_o = mf.mo_occ if not isinstance(mf.mo_occ, list) else mf.mo_occ[0] n_occ = int((mo_o > 0).sum()) @@ -706,19 +776,37 @@ def main(): print(f"HOMO-LUMO gap: {{gap:.4f}} eV") print() - # Save results next to the script so the path is predictable + # Save results next to the script so the path is predictable. + # AUDIT additional-concerns — quantui.orbital_visualization. + # generate_cube_file() requires 'mol_atom'/'mol_basis' (raises + # ValueError without them: "Re-run the calculation with the + # updated script template") and reads an optional 'mo_occ' to + # infer charge/spin for charged/open-shell molecules. This + # template used to save neither, so a cube could never be + # generated from a standalone-exported result without manually + # re-running with a different template that doesn't exist. results_path = str(Path(__file__).parent / 'results.npz') np.savez(results_path, energy=energy, mo_energy=np.array(mf.mo_energy), mo_coeff=np.array(mf.mo_coeff), - converged=mf.converged) + mo_occ=np.array(mf.mo_occ), + mol_atom=mol.atom, + mol_basis=str(mol.basis), + converged=mf.converged, + mp2_correlation_hartree=mp2_correlation, + ccsd_correlation_hartree=ccsd_correlation, + ccsd_t_correction_hartree=ccsd_t_correction, + cc_converged=cc_converged) print(f"Results saved to {{results_path}}") print("=" * 60) sys.exit(0) else: - print("ERROR: SCF did not converge!") + if not mf.converged: + print("ERROR: SCF did not converge!") + else: + print("ERROR: CCSD did not converge!") sys.exit(1) except Exception as e: diff --git a/quantui/freq_calc.py b/quantui/freq_calc.py index 6b9c282..d665081 100644 --- a/quantui/freq_calc.py +++ b/quantui/freq_calc.py @@ -72,7 +72,13 @@ class FreqResult: Attributes: energy_hartree: SCF energy at the input geometry in Hartrees. homo_lumo_gap_ev: HOMO-LUMO gap in eV, or ``None``. - converged: ``True`` if the SCF converged. + converged: ``True`` only when BOTH the SCF converged AND the + Hessian/harmonic-analysis step actually completed (AUDIT F15) + — e.g. a ROHF reference whose analytic Hessian PySCF doesn't + support on this path used to report ``converged=True`` with + ``frequencies_cm1=[]``, since this flag came solely from the + SCF. A frequency calculation with no computed Hessian is not a + successful frequency analysis, whatever the reference SCF did. n_iterations: Number of SCF macro-iterations. method: Calculation method (e.g. ``'RHF'``, ``'B3LYP'``). basis: Basis set (e.g. ``'STO-3G'``). @@ -107,7 +113,16 @@ class FreqResult: provide ``norm_mode``. """ mo_energy_hartree: Optional[List] = None + """Orbital energies for the Energies panel's diagram, in Hartrees. + + AUDIT additional-concerns — for an open-shell (UHF/UKS) reference, + this is the ALPHA-channel orbital energies only; the beta channel is + extracted and then discarded (``_moe[0]`` on a 2-D ``mf.mo_energy``). + Not a complete open-shell orbital spectrum. + """ mo_occ: Optional[List] = None + """Orbital occupations matching ``mo_energy_hartree`` — same + alpha-only caveat for an open-shell reference.""" pyscf_mol_atom: Optional[List] = None pyscf_mol_basis: Optional[str] = None density_fit: bool = False @@ -437,7 +452,7 @@ def _status(msg: str) -> None: # M-UX2 UXP2.10 — capture before maybe_apply_d3 can wrap/rename it. scf_variant = type(mf).__name__ mf.xc = resolve_xc(method) - mf = maybe_apply_d3(mf, method, progress_stream=stream) + mf, _ = maybe_apply_d3(mf, method, progress_stream=stream) # Density fitting (RI), opt-in (M-DF). Off by default. Applied to the main # SCF; the per-displacement inner SCFs below get the same treatment so the @@ -500,6 +515,10 @@ def _status(msg: str) -> None: _moe = mf.mo_energy _moo = mf.mo_occ if isinstance(_moe, (list, _np_mo.ndarray)) and hasattr(_moe[0], "__len__"): + # AUDIT additional-concerns — open-shell (UHF/UKS): mo_energy + # is (2, n_mo), alpha then beta. Only the alpha channel is + # kept for the orbital-diagram fields below; see + # mo_energy_hartree/mo_occ's field docstrings above. _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() @@ -530,6 +549,12 @@ def _status(msg: str) -> None: zpve_hartree: float = 0.0 displacements: Optional[List] = None thermo_data: Optional[ThermoData] = None + # AUDIT F15 — SCF convergence and Hessian/harmonic-analysis completion + # are separate facts; a caught exception in the try block below (e.g. + # ROHF's Hessian being unavailable on this path) must not leave the + # overall FreqResult reading "converged" with an empty + # frequencies_cm1. + _hessian_completed = False try: hess_obj = mf.Hessian() @@ -556,6 +581,11 @@ def _status(msg: str) -> None: else: frequencies_cm1.append(float(f.real if hasattr(f, "real") else f)) + # AUDIT F15 — the Hessian was built and harmonic_analysis() ran; a + # real frequency result exists regardless of whether the optional + # IR/Raman/thermo enrichment below succeeds. + _hessian_completed = True + # ZPVE = ½ · Σ ν_i (positive modes only), converted cm⁻¹ → Hartree zpve_hartree = sum(0.5 * f * _CM1_TO_HARTREE for f in frequencies_cm1 if f > 0) @@ -711,7 +741,12 @@ def _displaced_scf_dipole() -> _np_ir.ndarray: # can pin the contract. from quantui import freq_ir_workers as _ir_par - _cpu_count = os.cpu_count() or 1 + # AUDIT additional-concerns — available_cpu_count() honors + # SLURM_CPUS_PER_TASK / cgroup affinity instead of the raw + # os.cpu_count() (whole-machine core count), so this run + # doesn't oversubscribe a SLURM allocation smaller than the + # host it landed on. + _cpu_count = _ir_par.available_cpu_count() _use_parallel = _ir_par.parallel_enabled_for_run( cpu_count=_cpu_count, displacement_count=_ir_total_solves, @@ -794,6 +829,9 @@ def _displaced_scf_dipole() -> _np_ir.ndarray: _dm0_handle.name, _threads_each, _ckpt_items_dir, + mol.ecp, # AUDIT F05 + _density_fit_used, # AUDIT F19 + scf_rescue, # AUDIT F19 ), ) as _pool: # Submit all and store futures keyed by task @@ -990,13 +1028,21 @@ def _tv(v): f"Missing H or S in thermo dict (keys: {sorted(_tout.keys())})" ) _H = _tv(_H_raw) - _S = _tv(_S_raw) # J/(mol·K) + # PySCF's thermo() returns S_tot in Eh/K, not J/(mol·K) — despite + # the misleading local variable name this used to carry. Convert + # to J/(mol·K) for storage/display, and use the Eh/K value + # (matching H_hartree's units) to compute G = H - T*S. The old + # code stored the raw Eh/K number as S_jmol, then divided by + # _HARTREE_TO_JMOL again when forming G — nearly canceling the + # entropy term's contribution to G (see AUDIT F01). + _S_hartree_per_k = _tv(_S_raw) + _S_jmol = _S_hartree_per_k * _HARTREE_TO_JMOL _zpve = _tv(_Z_raw) if _Z_raw is not None else zpve_hartree - _G = _H - 298.15 * _S / _HARTREE_TO_JMOL + _G = _H - 298.15 * _S_hartree_per_k thermo_data = ThermoData( zpve_hartree=_zpve, H_hartree=_H, - S_jmol=_S, + S_jmol=_S_jmol, G_hartree=_G, ) _status("Frequency backend complete.") @@ -1028,7 +1074,9 @@ def _tv(v): return FreqResult( energy_hartree=energy_hartree, homo_lumo_gap_ev=homo_lumo_gap_ev, - converged=converged, + # AUDIT F15 — overall success requires the Hessian/harmonic-analysis + # step to have actually completed, not just the reference SCF. + converged=converged and _hessian_completed, n_iterations=n_iterations, method=method, basis=basis, diff --git a/quantui/freq_ir_workers.py b/quantui/freq_ir_workers.py index d42d8d8..6952ef0 100644 --- a/quantui/freq_ir_workers.py +++ b/quantui/freq_ir_workers.py @@ -56,6 +56,9 @@ def init_worker( dm0_pickle_path: str, omp_threads: int, checkpoint_items_dir: str | None = None, + ecp: dict | None = None, + density_fit: bool = False, + scf_rescue: bool = True, ) -> None: """ProcessPoolExecutor worker initializer. @@ -94,6 +97,32 @@ def init_worker( :func:`quantui.checkpoint.mark_item_done_at`. ``None`` when the run has no checkpoint (checkpointing is always optional — see :mod:`quantui.checkpoint`'s "never break a calculation" rule). + ecp: + The reference molecule's ``mol.ecp`` mapping (AUDIT F05) — + ``{element: basis}`` for elements whose basis carries an effective + core potential (e.g. LANL2DZ/def2 on heavy atoms), or ``{}``/``None`` + for an all-electron basis. Without this, a worker rebuilding the + ``Mole`` from ``atom_str``/``basis``/``charge``/``spin`` alone would + silently run the ECP atoms all-electron instead — a different + Hamiltonian (more electrons, no core potential), not just numerical + noise. See :func:`quantui.inorganic_guards.ecp_for_basis`. + density_fit: + AUDIT F19 — whether the *reference* SCF used density fitting + (M-DF). The serial IR loop matches this via + ``try_density_fit(mf, enabled=density_fit_used)`` before every + displaced SCF; this worker previously had no such parameter at + all, so every parallel displacement ran without density fitting + even when the reference (and the serial fallback) used it — + silently changing the numerical approximation, not just its + speed, whenever the user opted into parallel IR on a fitted + calculation. + scf_rescue: + AUDIT F19 — whether :func:`quantui.scf_robust.run_scf_with_rescue` + may apply its convergence-rescue ladder for this run. The serial + loop threads the caller's ``scf_rescue`` flag through + (``run_scf_with_rescue(_mf_d, dm0=_dm0, rescue=scf_rescue)``); this + worker previously hardcoded the rescue default (``True``) + regardless of what the caller requested. """ # Order matters: set env vars before any NumPy / PySCF import. threads = str(int(omp_threads)) @@ -115,6 +144,9 @@ def init_worker( xc=xc, dm0=dm0, checkpoint_items_dir=checkpoint_items_dir, + ecp=ecp or {}, + density_fit=bool(density_fit), + scf_rescue=bool(scf_rescue), ) @@ -157,6 +189,10 @@ def run_displaced_scf(item_id: str, coords_bohr_flat) -> Any: mol = gto.Mole() mol.atom = state["atom_str"] mol.basis = state["basis"] + # AUDIT F05 — without this, a heavy-element ECP system (e.g. + # NaH/LANL2DZ) silently runs all-electron here: a different + # Hamiltonian than the reference calculation, not just numerical noise. + mol.ecp = state.get("ecp") or {} mol.charge = state["charge"] mol.spin = state["spin"] mol.verbose = 0 @@ -183,9 +219,20 @@ def run_displaced_scf(item_id: str, coords_bohr_flat) -> Any: else: mf = scf.UHF(mol) if dm0_is_unrestricted else scf.RHF(mol) mf.verbose = 0 + # AUDIT F19 — match the reference SCF's density-fitting choice, exactly + # like the serial loop's `_try_density_fit(_mf_d, enabled=_density_fit_used)` + # (freq_calc.py). Without this, enabling parallel IR silently ran every + # displaced SCF without density fitting whenever the reference was fitted + # — a different numerical approximation, not merely a speed difference. + from .density_fitting import try_density_fit + + mf, _ = try_density_fit(mf, enabled=bool(state.get("density_fit", False))) from .scf_robust import run_scf_with_rescue - run_scf_with_rescue(mf, dm0=dm0) + # AUDIT F19 — honor the caller's scf_rescue choice instead of always + # taking run_scf_with_rescue's default (True); the serial loop already + # threads scf_rescue through as `rescue=scf_rescue`. + run_scf_with_rescue(mf, dm0=dm0, rescue=bool(state.get("scf_rescue", True))) dipole = np.array(mf.dip_moment(verbose=0)) items_dir = state.get("checkpoint_items_dir") @@ -200,6 +247,45 @@ def run_displaced_scf(item_id: str, coords_bohr_flat) -> Any: return dipole +def available_cpu_count() -> int: + """CPU budget for sizing the parallel IR/Raman worker pool. + + AUDIT additional-concerns — the driver used to size the worker pool + from a bare ``os.cpu_count() or 1``, which reports the WHOLE + machine's core count regardless of any SLURM allocation or cgroup/ + container CPU limit. On a 128-core host with an 8-core SLURM + allocation and 18 displacements, that picked 18 workers (min(64, 18)) + x 7 threads each — wildly oversubscribing the actual allocation. + + Precedence, most-authoritative first: + + 1. ``SLURM_CPUS_PER_TASK`` — the CPU count SLURM itself assigned to + this task. Authoritative even when the cluster does not enforce + cgroup CPU limits (many don't), which is exactly the case + ``os.sched_getaffinity``/``os.cpu_count()`` cannot see. + 2. ``os.sched_getaffinity(0)`` (Linux only) — respects a cgroup or + container CPU limit when one IS enforced and SLURM's own env var + is absent (e.g. a non-SLURM containerized deployment). + 3. ``os.cpu_count()`` — final fallback (Windows/macOS, or a + restricted environment without ``sched_getaffinity``). + + Every path floors at 1; never raises. + """ + slurm_cpus = os.environ.get("SLURM_CPUS_PER_TASK") + if slurm_cpus is not None: + try: + n = int(slurm_cpus) + if n > 0: + return n + except ValueError: + pass + try: + return len(os.sched_getaffinity(0)) # type: ignore[attr-defined] + except (AttributeError, OSError, NotImplementedError): + pass + return os.cpu_count() or 1 + + def freq_parallel_opt_in() -> bool: """Return whether parallel IR displacements are enabled for the next run.""" return _freq_parallel_opt_in() diff --git a/quantui/freq_raman_workers.py b/quantui/freq_raman_workers.py index f088561..43ba369 100644 --- a/quantui/freq_raman_workers.py +++ b/quantui/freq_raman_workers.py @@ -24,6 +24,8 @@ def init_raman_worker( dm0_is_unrestricted: bool, density_fit_used: bool, checkpoint_items_dir: str | None = None, + ecp: dict | None = None, + scf_rescue: bool = True, ) -> None: """Worker initializer — same threading discipline as IR workers. @@ -33,6 +35,14 @@ def init_raman_worker( polarizability is durably recorded via :func:`quantui.checkpoint.mark_item_done_at`. ``None`` when the run has no checkpoint. + + ``ecp`` (AUDIT F05): the reference molecule's ``mol.ecp`` mapping — see + :func:`quantui.freq_ir_workers.init_worker`'s docstring. Without it, a + heavy-element ECP system runs all-electron in this worker instead. + + ``scf_rescue`` (AUDIT F19): whether ``run_scf_with_rescue`` may apply + its convergence-rescue ladder. Previously hardcoded to the default + (``True``) regardless of the caller's ``scf_rescue`` choice. """ import os import pickle @@ -56,6 +66,8 @@ def init_raman_worker( dm0_is_unrestricted=bool(dm0_is_unrestricted), density_fit_used=bool(density_fit_used), checkpoint_items_dir=checkpoint_items_dir, + ecp=ecp or {}, + scf_rescue=bool(scf_rescue), ) @@ -92,6 +104,9 @@ def run_displaced_polarizability(item_id: str, coords_bohr_flat) -> list[list[fl mol = gto.Mole() mol.atom = state["atom_str"] mol.basis = state["basis"] + # AUDIT F05 — without this, a heavy-element ECP system (e.g. + # NaH/LANL2DZ) silently runs all-electron here. + mol.ecp = state.get("ecp") or {} mol.charge = state["charge"] mol.spin = state["spin"] mol.verbose = 0 @@ -113,7 +128,9 @@ def run_displaced_polarizability(item_id: str, coords_bohr_flat) -> list[list[fl mf, _ = _try_density_fit(mf, enabled=bool(state.get("density_fit_used"))) from .scf_robust import run_scf_with_rescue - run_scf_with_rescue(mf, dm0=dm0) + # AUDIT F19 — honor the caller's scf_rescue choice instead of always + # taking run_scf_with_rescue's default (True). + run_scf_with_rescue(mf, dm0=dm0, rescue=bool(state.get("scf_rescue", True))) pol_mod = _polarizability_module(mol, dm0_is_unrestricted) alpha = np.asarray(pol_mod.polarizability(pol_mod.Polarizability(mf)), dtype=float) diff --git a/quantui/ir_plot.py b/quantui/ir_plot.py index 6060574..8161d4c 100644 --- a/quantui/ir_plot.py +++ b/quantui/ir_plot.py @@ -9,18 +9,66 @@ from quantui.ir_plot import plot_ir_spectrum fig = plot_ir_spectrum(result.frequencies_cm1, result.ir_intensities) fig = plot_ir_spectrum(freqs, intensities, mode="broadened", fwhm=30.0) + +AUDIT additional-concerns — broadened-mode normalization +---------------------------------------------------------- +The "broadened" Lorentzian kernel below is HEIGHT-normalized: each peak's +value at its own center equals the supplied intensity, exactly like the +stick plot, so switching FWHM only changes peak width, never peak height +— the same convention used by Gaussian/ORCA-style broadened spectra and +by :mod:`quantui.raman_plot` and the UV-Vis broadening in +app_visualization.py. This is deliberate, NOT a bug: it is not an +area-normalized spectral density, so the AREA under a broadened peak +scales with FWHM (∫ γ²/((x-x0)²+γ²) dx = πγ) even though the peak height +does not. A true physical spectral-density plot (area proportional to +the supplied intensity, independent of the chosen FWHM) would need an +area-normalized Lorentzian (dividing by πγ) and different y-axis units — +intentionally out of scope here, since the height-preserving convention +is what students and most external QM software display by default. """ from __future__ import annotations -from typing import List, Optional +from typing import List, Optional, cast import numpy as np import plotly.graph_objects as go -# x-axis range is low → high wavenumber (user-facing convention in QuantUI) -_XRANGE = [400, 4000] -_XGRID = np.arange(400, 4001, 1.0) # 1 cm⁻¹ resolution for broadened mode +# x-axis range is low → high wavenumber (user-facing convention in QuantUI). +# AUDIT additional-concerns — this used to be the FIXED plot range/grid +# regardless of the actual data: a real water/STO-3G calculation has O-H +# stretches at 4486.7/4788.3 cm⁻¹, which used to be clipped off the right +# edge in stick mode (Plotly's xaxis.range) and never even entered the +# broadened kernel (evaluated only on this fixed grid) — computed modes +# silently invisible. ``_default_xrange``/``_grid_for_range`` below widen +# this default window to always cover every real (positive) frequency +# actually present, while leaving the familiar 400–4000 cm⁻¹ look +# untouched for the common case where every mode already falls inside it. +_DEFAULT_XRANGE = [400, 4000] + + +def _default_xrange(freqs_real: tuple) -> List[float]: + """[xmin, xmax] covering 400–4000 cm⁻¹ AND every real frequency present. + + A fixed margin keeps a peak sitting exactly at the edge from being + clipped by the axis border or cut off mid-lineshape in broadened mode. + """ + if not freqs_real: + return [float(_DEFAULT_XRANGE[0]), float(_DEFAULT_XRANGE[1])] + margin = 100.0 + lo = min(_DEFAULT_XRANGE[0], min(freqs_real) - margin) + hi = max(_DEFAULT_XRANGE[1], max(freqs_real) + margin) + return [float(lo), float(hi)] + + +def _grid_for_range(xrange: List[float]) -> np.ndarray: + # numpy's stubs (as pinned: mypy~=1.10.0) resolve this call to `Any` + # rather than `ndarray[Any, dtype[Any]]` for non-literal float bounds + # — pyscf has no type stubs (ignore_missing_imports) for the same + # underlying reason elsewhere in this codebase, and the fix there is + # the same: an explicit cast documents the real, known return type + # instead of silencing the check. + return cast(np.ndarray, np.arange(xrange[0], xrange[1] + 1.0, 1.0)) def plot_ir_spectrum( @@ -48,11 +96,13 @@ def plot_ir_spectrum( in a :class:`~plotly.graph_objects.FigureWidget`. """ real_pairs = [(f, i) for f, i in zip(frequencies, intensities) if f > 0] + freqs_real_for_range = tuple(f for f, _ in real_pairs) + xrange = _default_xrange(freqs_real_for_range) _base_layout = dict( xaxis=dict( title="Wavenumber (cm⁻¹)", - range=_XRANGE, + range=xrange, showgrid=True, gridcolor="#e5e7eb", ), @@ -78,14 +128,15 @@ def plot_ir_spectrum( freqs_real, ints_real = zip(*real_pairs) if mode == "broadened": + _xgrid = _grid_for_range(xrange) half_gamma = fwhm / 2.0 - y_broad = np.zeros_like(_XGRID) + y_broad = np.zeros_like(_xgrid) for nu0, inten in zip(freqs_real, ints_real): - y_broad += inten * half_gamma**2 / ((_XGRID - nu0) ** 2 + half_gamma**2) + y_broad += inten * half_gamma**2 / ((_xgrid - nu0) ** 2 + half_gamma**2) fig.add_trace( go.Scatter( - x=_XGRID, + x=_xgrid, y=y_broad, mode="lines", line=dict(color="#2563eb", width=1.5), diff --git a/quantui/log_utils.py b/quantui/log_utils.py index 2399db9..4c8c692 100644 --- a/quantui/log_utils.py +++ b/quantui/log_utils.py @@ -463,12 +463,17 @@ def format_log_footer( zpve = getattr(result, "zpve_hartree", None) n_steps = getattr(result, "n_steps", None) # OptResult - # Convergence line + # Convergence line. AUDIT F07/F08/F15 — .converged on several result + # types now folds in more than the reference SCF (CCSD's own + # amplitude convergence, TD-DFT's per-root convergence, the + # Hessian/harmonic-analysis step actually completing), so this must + # not claim "SCF" specifically — that would misleadingly blame the + # reference SCF for e.g. a Hessian failure with a converged SCF. if converged is not None: tick = "✓" if converged else "✗" conv_word = "converged" if converged else "did NOT converge" iter_str = f" | Iterations: {n_iter}" if n_iter is not None else "" - lines.append(f" {tick} SCF {conv_word}{iter_str}") + lines.append(f" {tick} Result {conv_word}{iter_str}") if n_steps is not None: lines.append(f" Geometry optimization: {n_steps} steps") diff --git a/quantui/nmr_calc.py b/quantui/nmr_calc.py index 9db3fac..1523700 100644 --- a/quantui/nmr_calc.py +++ b/quantui/nmr_calc.py @@ -375,7 +375,7 @@ def _run_nmr_calc_body( # M-UX2 UXP2.10 — capture before maybe_apply_d3 can wrap/rename it. scf_variant = type(mf).__name__ mf.xc = resolve_xc(method) - mf = maybe_apply_d3(mf, method, progress_stream=stream) + mf, _ = maybe_apply_d3(mf, method, progress_stream=stream) # Density fitting (RI), opt-in (M-DF). Off by default. See DF.5: DF shifts # absolute shieldings, but chemical shifts are differences so the error diff --git a/quantui/optimizer.py b/quantui/optimizer.py index 57adf4e..0e85a8c 100644 --- a/quantui/optimizer.py +++ b/quantui/optimizer.py @@ -120,6 +120,10 @@ def __init__( self.status_label = status_label self.expected_steps = expected_steps # history-based ~N prior self._eval_count = 0 + # AUDIT F04 — None (method doesn't use D3), True (applied at + # every step so far), or False (pyscf.dftd3 unavailable at some + # step — sticky once seen, since it can't un-happen mid-run). + self.dispersion_applied: Optional[bool] = None def calculate( self, @@ -197,7 +201,16 @@ def calculate( mf = dft.RKS(mol) if mol.spin == 0 else dft.UKS(mol) mf.xc = resolve_xc(self.method) - mf = maybe_apply_d3(mf, self.method) + # AUDIT F04 — this call used to omit progress_stream + # entirely, so a missing pyscf.dftd3 gave NO warning + # anywhere on the optimizer path (unlike every other DFT + # entry point). maybe_apply_d3 now also always logs, but + # pass the stream too so the user sees it in-app. + mf, _dispersion_applied = maybe_apply_d3( + mf, self.method, progress_stream=self.progress_stream + ) + if self.dispersion_applied is not False: + self.dispersion_applied = _dispersion_applied # Density fitting (RI), opt-in (M-DF). Off by default; applies to # every SCF in the optimization when the user enables it. @@ -227,6 +240,27 @@ def _scf_progress(envs, _k=_k) -> None: run_scf_with_rescue(mf, rescue=self.scf_rescue, stream=self.progress_stream) + # AUDIT F09 — BFGS previously accepted whatever gradient came + # back regardless of mf.converged, so it could satisfy its force + # criterion using an invalid electronic solution (verified: an + # H2 optimization near its minimum, with SCF limited to one + # cycle and rescue disabled, reported converged=True after + # three steps despite all four SCF evaluations being + # unconverged). run_scf_with_rescue has already exhausted every + # rescue stage by this point, so an unconverged mf here means + # this step's energy/forces are not physically meaningful — + # raise rather than hand them to ASE, which would silently bake + # them into the optimization trajectory (and, via BFGS's + # Hessian update, corrupt every subsequent step too). + if not bool(getattr(mf, "converged", False)): + raise RuntimeError( + f"SCF did not converge at optimization step " + f"{self._eval_count} — the resulting energy/forces are " + "not physically meaningful. Try scf_rescue=True " + "(default), a different starting geometry, or a " + "different basis/method." + ) + # Save final SCF state for orbital visualization self._last_mf = mf self._last_atom_list = _atom_list_for_cube @@ -289,6 +323,10 @@ class OptimizationResult: pyscf_mol_atom: Optional[Any] = None # atom list at final geometry (Angstrom) pyscf_mol_basis: Optional[str] = None density_fit: bool = False + # AUDIT F04 — mirrors SessionResult.dispersion_applied: None (method + # doesn't use D3), True/False (does, and pyscf.dftd3 was/wasn't + # importable during the optimization). + dispersion_applied: Optional[bool] = None # Final-geometry Mulliken / dipole — same fields SessionResult carries so # the Populations Analysis panel activates after a Geometry Opt too. atom_symbols: Optional[List[str]] = None @@ -868,6 +906,7 @@ def _report_opt_fraction() -> None: mulliken_charges=_opt_mulliken, dipole_moment_debye=_opt_dipole, dipole_vector_debye=_opt_dipole_vec, + dispersion_applied=getattr(atoms.calc, "dispersion_applied", None), ) diff --git a/quantui/orbital_visualization.py b/quantui/orbital_visualization.py index 6a0aa61..0fa391d 100644 --- a/quantui/orbital_visualization.py +++ b/quantui/orbital_visualization.py @@ -543,7 +543,9 @@ def orbital_summary_html(info: OrbitalInfo) -> str: def infer_charge_and_spin( - mol_atom: Optional[list], mo_occ: Optional[np.ndarray | list] + mol_atom: Optional[list], + mo_occ: Optional[np.ndarray | list], + basis: Optional[str] = None, ) -> Tuple[int, int]: """Infer ``(charge, spin)`` for a ``gto.Mole`` from atoms + MO occupations. @@ -556,12 +558,22 @@ def infer_charge_and_spin( 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). + - ``spin`` (PySCF's ``2S = n_alpha - n_beta``) is ``n_alpha - n_beta`` + when ``mo_occ`` is 2-D (UHF/UKS, shape ``(2, n_mo)``). When ``mo_occ`` + is 1-D, it is the count of singly-occupied orbitals (AUDIT F14): a + 1-D array is NOT necessarily closed-shell — ROHF is 1-D too (values + 2/1/0), and by the standard ROHF convention every singly-occupied + orbital is alpha, so that count IS 2S directly. This is 0 for a + genuine closed-shell RHF/RKS/MP2/CCSD/CCSD(T) reference (no singly- + occupied orbitals) and correct for ROHF, without needing to know + which SCF variant actually produced ``mo_occ``. + - ``charge`` is the effective nuclear charge (sum of atomic numbers in + ``mol_atom``, minus each ECP's core-electron count when ``basis`` is + given — AUDIT F14: a bare atomic-number sum overcounts an ECP system + by however many core electrons the ECP replaced, since ``mo_occ`` + only counts the explicit/valence electrons the calculation used) + 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. @@ -577,10 +589,24 @@ def infer_charge_and_spin( spin = int(round(n_alpha - n_beta)) n_electrons = n_alpha + n_beta else: - spin = 0 + spin = int(np.sum(np.isclose(occ, 1.0))) n_electrons = float(occ.sum()) - nuclear_charge = sum(ATOMIC_NUMBERS.get(sym, 0) for sym, _ in mol_atom) + nuclear_charge = 0 + for sym, _pos in mol_atom: + z = ATOMIC_NUMBERS.get(sym, 0) + core_electrons = 0 + if basis: + try: + from pyscf import gto as _gto + + _ecp_data = _gto.basis.load_ecp(basis, sym) + if _ecp_data: + core_electrons = int(_ecp_data[0]) + except Exception: + core_electrons = 0 + nuclear_charge += z - core_electrons + charge = int(round(nuclear_charge - n_electrons)) return charge, spin @@ -720,7 +746,7 @@ def generate_cube_file( for tok in atom_str.replace(";", "\n").splitlines() if tok.strip() ] - charge, spin = infer_charge_and_spin(parsed_atoms, mo_occ) + charge, spin = infer_charge_and_spin(parsed_atoms, mo_occ, basis=basis_str) mol = gto.M( atom=atom_str, basis=basis_str, unit="Angstrom", charge=charge, spin=spin diff --git a/quantui/pes_scan.py b/quantui/pes_scan.py index 55eebcd..4181419 100644 --- a/quantui/pes_scan.py +++ b/quantui/pes_scan.py @@ -26,7 +26,7 @@ import math import sys from dataclasses import dataclass -from typing import IO, Any, Dict, List, Optional +from typing import IO, Any, Dict, List, Optional, Sequence from .ase_bridge import ASE_AVAILABLE, atoms_to_molecule, molecule_to_atoms from .molecule import Molecule @@ -175,6 +175,8 @@ def _reuse_scan_point( value: float, atoms: Any, molecule: Molecule, + scan_type: str, + atom_indices: Sequence[int], ) -> Optional[tuple]: """Return ``(energy_hartree, molecule)`` for a reusable point, else ``None``. @@ -186,12 +188,30 @@ def _reuse_scan_point( Returns ``None`` for anything questionable (no record, a coordinate value that no longer matches, malformed geometry). Recomputing a point is cheap next to trusting a mismatched one. + + AUDIT F10: a cached point's ``value`` alone does not identify WHICH + coordinate it was computed for — a banked O-H bond-scan point at 1.0 A + and a fresh H-H bond scan targeting 1.0 A compared equal on ``value`` + alone, so the O-H point was reused and returned as if it were the H-H + point (its real geometry, e.g. an O-H distance of 1.0 A, has whatever + H-H distance that geometry happens to have — not 1.0 A). Now the + record's own ``scan_type``/``atom_indices`` must match the current + scan's, AND the coordinate actually measured from the restored geometry + (not just the record's self-reported ``value``) must match the target + — catching a stale/mismatched record even if ``scan_type``/ + ``atom_indices`` were themselves corrupted or from an older checkpoint + schema that didn't record them at all. """ if not record: return None try: if abs(float(record["value"]) - float(value)) > 1e-9: return None + if str(record.get("scan_type", scan_type)) != str(scan_type): + return None + _rec_indices = record.get("atom_indices") + if _rec_indices is not None and list(_rec_indices) != list(atom_indices): + return None energy_ha = float(record["energy_hartree"]) symbols = [str(a) for a in record["atoms"]] coords = [[float(c) for c in row] for row in record["coordinates"]] @@ -210,6 +230,32 @@ def _reuse_scan_point( atoms.set_positions(coords) except Exception: # noqa: BLE001 — a stale live geometry is not fatal logger.debug("could not restore ASE positions for a reused scan point") + return None + + # Cross-check: measure the actual coordinate from the restored geometry + # itself, independent of anything the record claims about itself. + try: + i1, i2 = atom_indices[0], atom_indices[1] + if scan_type == "bond": + _actual = float(atoms.get_distance(i1, i2)) + _diff = abs(_actual - float(value)) + _tol = 1e-3 + else: + if scan_type == "angle": + _actual = float(atoms.get_angle(i1, i2, atom_indices[2])) + else: # dihedral + _actual = float( + atoms.get_dihedral(i1, i2, atom_indices[2], atom_indices[3]) + ) + # Angles/dihedrals wrap at 360 degrees (e.g. -170 and 190 are the + # same angle) — compare via the shortest angular distance. + _diff = abs((_actual - float(value) + 180.0) % 360.0 - 180.0) + _tol = 1e-2 + if _diff > _tol: + return None + except Exception: # noqa: BLE001 — malformed geometry: don't trust it + return None + return energy_ha, point_molecule @@ -420,7 +466,14 @@ def run_pes_scan( # (points already done / total) for the self-correcting time estimate. from .log_utils import emit_progress, emit_status - _reused = _reuse_scan_point(_cached_points.get(step_num), val, atoms, molecule) + _reused = _reuse_scan_point( + _cached_points.get(step_num), + val, + atoms, + molecule, + scan_type, + atom_indices, + ) if _reused is not None: _energy_ha, _mol_at_point = _reused energies_hartree.append(_energy_ha) @@ -517,6 +570,12 @@ def run_pes_scan( { "index": step_num, "value": float(val), + # AUDIT F10 — identifies WHICH coordinate this point + # belongs to, so a later scan of a different bond/ + # angle/dihedral (or different atoms) can't have a + # coincidentally-matching "value" reuse this point. + "scan_type": scan_type, + "atom_indices": list(atom_indices), "energy_hartree": float(e_ha), "ok": bool(ok), "atoms": list(mol_at_point.atoms), diff --git a/quantui/raman_calc.py b/quantui/raman_calc.py index 2f83b3f..52d721c 100644 --- a/quantui/raman_calc.py +++ b/quantui/raman_calc.py @@ -244,7 +244,10 @@ def _displaced_alpha(atom_idx: int, ax: int, sign: int) -> np.ndarray: from quantui import freq_ir_workers as _ir_par from quantui import freq_raman_workers as _ram_par - _cpu_count = os.cpu_count() or 1 + # AUDIT additional-concerns — see freq_calc.py's matching comment: + # available_cpu_count() honors SLURM_CPUS_PER_TASK / cgroup affinity + # instead of the whole machine's os.cpu_count(). + _cpu_count = _ir_par.available_cpu_count() _use_parallel = _ir_par.parallel_enabled_for_run( cpu_count=_cpu_count, displacement_count=_total, @@ -313,6 +316,8 @@ def _displaced_alpha(atom_idx: int, ax: int, sign: int) -> np.ndarray: dm0_is_unrestricted, density_fit_used, _ckpt_items_dir, + mol.ecp, # AUDIT F05 + scf_rescue, # AUDIT F19 ), ) as _pool: _futs = { @@ -377,7 +382,16 @@ def _displaced_alpha(atom_idx: int, ax: int, sign: int) -> np.ndarray: mol.set_geom_(_coords0, unit="Bohr") mol.verbose = _mol_v - dalpha_ang = dalpha / _BOHR_TO_ANG + # dalpha is d(alpha[a0^3]) / d(x[Bohr]) — polarizability in atomic units + # (a0^3), displacement in Bohr. Converting to d(alpha[A^3]) / d(x[A]) + # needs the numerator rescaled by BOHR_TO_ANGSTROM**3 (a0^3 -> A^3) *and* + # the denominator by BOHR_TO_ANGSTROM (Bohr -> A): a net factor of + # BOHR_TO_ANGSTROM**2. The old code divided by a single + # BOHR_TO_ANGSTROM, rescaling only the denominator and leaving the + # numerator in a0^3 instead of A^3 — a missing BOHR_TO_ANGSTROM**3 + # factor in the derivative, which becomes BOHR_TO_ANGSTROM**6 once + # squared into the Raman activity: ~45.54x too large (AUDIT F02). + dalpha_ang = dalpha * (_BOHR_TO_ANG**2) nm = np.asarray(displacements, dtype=float) if nm.ndim == 2: nm = nm.reshape(nm.shape[0], _n_atoms, 3) diff --git a/quantui/raman_plot.py b/quantui/raman_plot.py index 617019e..eb5d728 100644 --- a/quantui/raman_plot.py +++ b/quantui/raman_plot.py @@ -8,6 +8,11 @@ from quantui.raman_plot import plot_raman_spectrum fig = plot_raman_spectrum(result.frequencies_cm1, result.raman_activities) + +AUDIT additional-concerns — see ir_plot.py's module docstring: the +broadened-mode Lorentzian here is the same height-normalized convention +(peak height = supplied activity, area scales with FWHM), deliberately, +not a bug. """ from __future__ import annotations @@ -17,7 +22,7 @@ import numpy as np import plotly.graph_objects as go -from quantui.ir_plot import _XGRID, _XRANGE +from quantui.ir_plot import _default_xrange, _grid_for_range def plot_raman_spectrum( @@ -30,11 +35,16 @@ def plot_raman_spectrum( ) -> go.Figure: """Return a Plotly figure for the Raman scattering spectrum.""" real_pairs = [(f, a) for f, a in zip(frequencies, activities) if f > 0] + # AUDIT additional-concerns — see ir_plot.py: the x-range must cover + # every real (positive) frequency present, not just the fixed + # 400-4000 cm⁻¹ default, or a real high-frequency mode (O-H stretches + # routinely sit above 4000 cm⁻¹) is silently clipped off the plot. + xrange = _default_xrange(tuple(f for f, _ in real_pairs)) _base_layout = dict( xaxis=dict( title="Wavenumber (cm⁻¹)", - range=_XRANGE, + range=xrange, showgrid=True, gridcolor="#e5e7eb", ), @@ -60,14 +70,15 @@ def plot_raman_spectrum( freqs_real, acts_real = zip(*real_pairs) if mode == "broadened": + _xgrid = _grid_for_range(xrange) half_gamma = fwhm / 2.0 - y_broad = np.zeros_like(_XGRID) + y_broad = np.zeros_like(_xgrid) for nu0, act in zip(freqs_real, acts_real): - y_broad += act * half_gamma**2 / ((_XGRID - nu0) ** 2 + half_gamma**2) + y_broad += act * half_gamma**2 / ((_xgrid - nu0) ** 2 + half_gamma**2) fig.add_trace( go.Scatter( - x=_XGRID, + x=_xgrid, y=y_broad, mode="lines", line=dict(color="#059669", width=1.5), diff --git a/quantui/reorganization_energy.py b/quantui/reorganization_energy.py index 3287b26..97c6b47 100644 --- a/quantui/reorganization_energy.py +++ b/quantui/reorganization_energy.py @@ -31,6 +31,26 @@ the rest of QuantUI: :func:`quantui.optimizer.optimize_geometry` for the relaxations and :func:`quantui.session_calc.run_in_session` for the single-point cross evaluations. + +Scope and approximations (AUDIT additional-concerns) +----------------------------------------------------- +The core four-point energy differences above are exact given the four +single-point energies; the two approximations worth being explicit about +are in how those energies are obtained: + +* **Ion spin state.** :func:`_ion_multiplicity` picks the minimal valid + multiplicity from electron-count parity (1 for even, 2 for odd) — a + convenient default, NOT a determination of the true ground-state spin. + This is fine for most organic radicals/closed-shell ions but can be + wrong for a transition-metal ion, where the actual ground state may be + higher-spin. Users wanting a specific (e.g. high-spin) ion state must + build that calculation manually. +* **Solvent scope.** An optional PCM ``solvent`` (see + :func:`run_reorganization_energy`) applies only to the four single-point + energies, not to the geometry relaxations that produce ``R_neutral``/ + ``R_ion`` — those are always gas-phase optimizations. λ therefore mixes + a solvent-phase energy with a gas-phase-optimized geometry, not a fully + solvent-consistent (geometry-and-energy) treatment. """ from __future__ import annotations diff --git a/quantui/session_calc.py b/quantui/session_calc.py index d91360f..2a55cb5 100644 --- a/quantui/session_calc.py +++ b/quantui/session_calc.py @@ -47,7 +47,12 @@ class SessionResult: energy_hartree: Total SCF energy in Hartrees. homo_lumo_gap_ev: HOMO-LUMO gap in electronvolts, or ``None`` if the gap cannot be determined (e.g. open-shell UHF with complex orbital - occupations, or too few occupied orbitals). + occupations, or too few occupied orbitals). AUDIT additional- + concerns: for a 2-D ``mo_energy`` (UHF/UKS), this is the ALPHA- + channel gap only — the beta channel is not computed or reported + here. The result card labels this "HOMO-LUMO gap (α)" for an + open-shell reference; treat it as a single-channel descriptor, + not a complete open-shell orbital spectrum. converged: ``True`` if the SCF iterations reached the convergence threshold; ``False`` if the maximum iteration count was hit. n_iterations: Number of SCF macro-iterations completed. May be @@ -91,6 +96,12 @@ class SessionResult: # method is ``"CCSD(T)"``. ``None`` for plain CCSD. Again, included in # ``energy_hartree`` when set. ccsd_t_correction_hartree: Optional[float] = None + # Whether the CCSD amplitude iterations themselves converged (AUDIT + # F07). ``None`` unless method is ``"CCSD"``/``"CCSD(T)"``. ``converged`` + # above already folds this in (False whenever this is False), but a + # caller inspecting *why* wants this separated from the HF reference's + # own convergence. + cc_converged: Optional[bool] = None # GPU offload status. ``gpu_used`` is True only when the # SCF object was successfully migrated to gpu4pyscf for this run. # ``gpu_name`` carries the CUDA device name when ``gpu_used`` is True so @@ -101,6 +112,13 @@ class SessionResult: # ``False`` for exact four-centre integrals (the default) and for the # post-HF paths, which are never fitted here. density_fit: bool = False + # Whether Grimme D3 dispersion was actually applied (AUDIT F04). ``None`` + # when the method doesn't use D3 (e.g. RHF, or a functional whose + # dispersion is built in, like wB97X-D); ``True``/``False`` when it does + # and ``pyscf.dftd3`` was/wasn't importable. A ``False`` result is + # missing its dispersion correction even though ``method`` still reads + # e.g. "PBE-D3" — see :func:`maybe_apply_d3` and :meth:`summary`. + dispersion_applied: Optional[bool] = None solvent: Optional[str] = None mo_energy_hartree: Optional[Any] = None # np.ndarray (n_mo,) or (2, n_mo) UHF mo_occ: Optional[Any] = None # np.ndarray (n_mo,) or (2, n_mo) UHF @@ -136,6 +154,16 @@ def summary(self) -> str: ] if self.homo_lumo_gap_ev is not None: lines.append(f" HOMO-LUMO gap : {self.homo_lumo_gap_ev:.4f} eV") + if self.dispersion_applied is False: + lines.append( + f" ⚠️ {self.method} requires D3 dispersion, but pyscf.dftd3 " + "was unavailable — this result has NO dispersion correction." + ) + if self.cc_converged is False: + lines.append( + " ⚠️ CCSD amplitude iterations did NOT converge — the " + "correlation energy above is unreliable." + ) lines += [ "=" * 60, ( @@ -155,22 +183,33 @@ def summary(self) -> str: # Maps QuantUI display names → PySCF xc strings where they differ. # -# ``wB97X-D`` is a special case: PySCF + dftd3 cannot compose -# ``mf.xc = "wb97x-d"`` cleanly (it's on dftd3's black-list — see -# pyscf/pyscf#2069). The workaround that matches what our UI label -# already claims ("wB97X-D — Range-Separated Hybrid + D3 Dispersion") -# is to use the bare ``wb97x`` functional and apply D3 via dftd3 -# externally — same pattern as PBE-D3 below. This is D3, not the -# original Chai 2008 D2; the empirical dispersion energies differ by -# a few percent for most systems but the functional family is the same. +# ``wB97X-D`` is a special case, but NOT the one this table used to assume +# (AUDIT F03). PySCF rejects ``mf.xc = "wb97x-d"`` — but not because it needs +# an external dispersion correction composed on: PySCF's own xc_code parser +# (``pyscf.scf.dispersion.parse_dft``) black-lists the short "wb97x-d" / +# "wb97x_d" spellings specifically because they're ambiguous between the +# original Chai & Head-Gordon (2008) wB97X-D functional (its own built-in +# empirical dispersion, baked into the fit, no Grimme correction needed) and +# a Grimme-D3-corrected bare wB97X. Aliasing to bare ``wb97x`` and applying +# external Grimme D3 (as this table previously did) silently calculates a +# *different* functional: wb97x has omega=0.3 range separation, wb97x-d has +# omega=0.2 and different short-range exact exchange (confirmed via +# ``pyscf.dft.libxc.rsh_coeff``) — not just a different dispersion model. +# +# The actual wB97X-D functional is available directly under its full LibXC +# name, which PySCF's short-alias black-list does not intercept, and needs +# no external D3 wrapper (see ``_NEEDS_D3`` below): _XC_ALIAS: Dict[str, str] = { "M06-L": "m06l", - "wB97X-D": "wb97x", # bare functional; D3 applied via _NEEDS_D3 + "wB97X-D": "hyb_gga_xc_wb97x_d", # true Chai/Head-Gordon 2008 functional "CAM-B3LYP": "camb3lyp", "PBE-D3": "pbe", # base functional; D3 applied separately } # Methods that require Grimme D3 dispersion correction via pyscf.dftd3. -_NEEDS_D3: frozenset = frozenset({"PBE-D3", "wB97X-D"}) +# wB97X-D is NOT here: its dispersion is already part of the XC functional +# itself (see _XC_ALIAS comment above) — wrapping it in pyscf.dftd3 would +# double-count dispersion under a method that already includes its own. +_NEEDS_D3: frozenset = frozenset({"PBE-D3"}) def resolve_xc(method: str) -> str: @@ -207,18 +246,31 @@ def needs_d3(method: str) -> bool: def maybe_apply_d3(mf, method: str, progress_stream=None): """Wrap ``mf`` in ``pyscf.dftd3.dftd3(mf)`` if ``method`` requires D3. - Returns the (possibly wrapped) mf object. On ``pyscf.dftd3`` - ImportError, returns the original ``mf`` unmodified and surfaces - a warning via ``progress_stream`` (if provided) so the user sees - that the result is missing the dispersion correction. + Returns ``(mf, dispersion_applied)``: the (possibly wrapped) mf object, + and whether the D3 wrapper was actually applied. ``dispersion_applied`` + is ``True`` when D3 was applied, ``False`` when the method needs D3 but + ``pyscf.dftd3`` is unavailable (AUDIT F04 — the result is silently + missing its dispersion correction; callers should record this rather + than keep reporting the original method label as if uncorrected = + corrected), and ``None`` when the method doesn't use D3 at all. + + On ``pyscf.dftd3`` ImportError, always logs a warning (so every call + site is visible in logs even without a progress stream — the optimizer + path used to call this with no stream and so surfaced nothing at all), + and additionally surfaces the warning via ``progress_stream`` when one + is provided. """ if not needs_d3(method): - return mf + return mf, None try: from pyscf import dftd3 as _dftd3 - return _dftd3.dftd3(mf) + return _dftd3.dftd3(mf), True except ImportError: + logger.warning( + "pyscf.dftd3 not available — running %s without D3 correction.", + method, + ) if progress_stream is not None: try: progress_stream.write( @@ -227,7 +279,7 @@ def maybe_apply_d3(mf, method: str, progress_stream=None): ) except Exception: # noqa: BLE001 — cleanup (stream may be closed) pass - return mf + return mf, False def run_in_session( @@ -456,6 +508,7 @@ def _run_session_calc_body( # --- Select SCF method --- method_upper = method.upper() + dispersion_applied: Optional[bool] = None if method_upper == "RHF": mf = scf.RHF(mol) scf_variant = type(mf).__name__ @@ -500,7 +553,9 @@ def _run_session_calc_body( # GOTCHAS.md. scf_variant = type(mf).__name__ mf.xc = resolve_xc(method) - mf = maybe_apply_d3(mf, method, progress_stream=progress_stream) + mf, dispersion_applied = maybe_apply_d3( + mf, method, progress_stream=progress_stream + ) # --- Density fitting (RI), opt-in (M-DF) --- # Applied to the freshly built SCF object, BEFORE the PCM wrap and the GPU @@ -623,25 +678,38 @@ def _run_session_calc_body( f"({method}/{basis}): {exc}" ) from exc + # AUDIT F07 — post-HF work (MP2/CCSD/CCSD(T)) needs a converged + # reference; running it on an unconverged SCF's orbitals produces a + # correlation "correction" on top of a wrong Hamiltonian, not a small + # numerical difference. Checked once, right after the reference SCF, + # before any post-HF method is even attempted. + scf_converged = bool(getattr(mf, "converged", False)) + # --- MP2 correlation energy (post-HF) --- mp2_correlation_hartree: Optional[float] = None if method_upper == "MP2": - try: - from pyscf import mp as _mp - - emit_status(stream, "Running MP2 correlation…") - _mp2 = _mp.MP2(mf) - # verbose=5 surfaces integral-transform / kernel milestones for - # the live status label during the correlation step. - _mp2.verbose = 5 - _mp2.stdout = stream - _e_corr, _ = _mp2.kernel() - mp2_correlation_hartree = float(_e_corr) - energy_hartree += float(_e_corr) - except Exception as exc: - raise RuntimeError( - f"MP2 correction failed for {molecule.get_formula()}: {exc}" - ) from exc + if not scf_converged: + emit_status( + stream, + "Skipping MP2 — reference SCF did not converge.", + ) + else: + try: + from pyscf import mp as _mp + + emit_status(stream, "Running MP2 correlation…") + _mp2 = _mp.MP2(mf) + # verbose=5 surfaces integral-transform / kernel milestones for + # the live status label during the correlation step. + _mp2.verbose = 5 + _mp2.stdout = stream + _e_corr, _ = _mp2.kernel() + mp2_correlation_hartree = float(_e_corr) + energy_hartree += float(_e_corr) + except Exception as exc: + raise RuntimeError( + f"MP2 correction failed for {molecule.get_formula()}: {exc}" + ) from exc # --- Coupled cluster correlation --- # CCSD adds singles + doubles excitations on top of the RHF reference; @@ -650,37 +718,65 @@ def _run_session_calc_body( # show the HF reference + correlation breakdown (mirrors the MP2 path). ccsd_correlation_hartree: Optional[float] = None ccsd_t_correction_hartree: Optional[float] = None + # AUDIT F07 — CCSD's own amplitude convergence, tracked separately from + # the HF reference's. The old code accepted _ccsd_obj.kernel()'s + # correlation energy unconditionally and never checked + # _ccsd_obj.converged, so a real one-iteration-limited non-convergence + # (verified: energy -75.007987575 Eh, converged=False) was reported as + # converged=True (from the HF reference alone). CCSD(T) also used to + # proceed to the triples correction regardless of CCSD's convergence. + cc_converged: Optional[bool] = None if method_upper in ("CCSD", "CCSD(T)"): - try: - from pyscf import cc as _cc - - emit_status(stream, "Running CCSD correlation…") - _ccsd_obj = _cc.CCSD(mf) - _ccsd_obj.verbose = 4 - _ccsd_obj.stdout = stream - _e_corr_ccsd, _t1, _t2 = _ccsd_obj.kernel() - ccsd_correlation_hartree = float(_e_corr_ccsd) - energy_hartree += float(_e_corr_ccsd) - except Exception as exc: - raise RuntimeError( - f"CCSD correction failed for {molecule.get_formula()}: {exc}" - ) from exc - if method_upper == "CCSD(T)": + if not scf_converged: + emit_status( + stream, + "Skipping CCSD — reference SCF did not converge.", + ) + else: try: - emit_status(stream, "Computing CCSD(T) triples…") + from pyscf import cc as _cc + + emit_status(stream, "Running CCSD correlation…") + _ccsd_obj = _cc.CCSD(mf) _ccsd_obj.verbose = 4 _ccsd_obj.stdout = stream - _e_t = _ccsd_obj.ccsd_t() - ccsd_t_correction_hartree = float(_e_t) - energy_hartree += float(_e_t) + _e_corr_ccsd, _t1, _t2 = _ccsd_obj.kernel() + cc_converged = bool(getattr(_ccsd_obj, "converged", False)) + ccsd_correlation_hartree = float(_e_corr_ccsd) + energy_hartree += float(_e_corr_ccsd) except Exception as exc: raise RuntimeError( - f"CCSD(T) triples correction failed " - f"for {molecule.get_formula()}: {exc}" + f"CCSD correction failed for {molecule.get_formula()}: {exc}" ) from exc + if method_upper == "CCSD(T)": + if not cc_converged: + emit_status( + stream, + "Skipping CCSD(T) triples — CCSD amplitudes did " + "not converge.", + ) + else: + try: + emit_status(stream, "Computing CCSD(T) triples…") + _ccsd_obj.verbose = 4 + _ccsd_obj.stdout = stream + _e_t = _ccsd_obj.ccsd_t() + ccsd_t_correction_hartree = float(_e_t) + energy_hartree += float(_e_t) + except Exception as exc: + raise RuntimeError( + f"CCSD(T) triples correction failed " + f"for {molecule.get_formula()}: {exc}" + ) from exc # --- Extract results from the mean-field object --- - converged = bool(getattr(mf, "converged", False)) + # AUDIT F07 — overall convergence must reflect every stage that ran: + # the HF reference, and (when requested) CCSD's own amplitude solve. + # A method that needed CC and didn't get a converged one is not a + # converged result, regardless of what the HF reference alone did. + converged = scf_converged + if method_upper in ("CCSD", "CCSD(T)"): + converged = scf_converged and bool(cc_converged) n_iterations = int(getattr(mf, "cycles", -1)) import numpy as _np @@ -729,6 +825,15 @@ def _to_numpy_array(arr: Any) -> Any: mulliken_charges: Optional[List[float]] = None dipole_moment_debye: Optional[float] = None + # AUDIT additional-concerns — for MP2/CCSD/CCSD(T), ``mf`` here is + # still the HF reference object (the post-HF correlation energy is + # computed separately and added to ``energy_hartree``; no correlated + # density is built for these methods). Both properties below are + # therefore HF-reference values even when method='CCSD(T)', NOT a + # correlated dipole/population — the result card labels them + # accordingly (_result_extra_rows' "HF reference" note) rather than + # presenting them as an unqualified property of the requested method. + # # 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"`` @@ -834,9 +939,11 @@ def _to_numpy_array(arr: Any) -> Any: mp2_correlation_hartree=mp2_correlation_hartree, ccsd_correlation_hartree=ccsd_correlation_hartree, ccsd_t_correction_hartree=ccsd_t_correction_hartree, + cc_converged=cc_converged, gpu_used=gpu_used, gpu_name=gpu_name, density_fit=density_fit_used, + dispersion_applied=dispersion_applied, solvent=solvent, mo_energy_hartree=_mo_energy_ha_arr, mo_occ=_mo_occ_arr, diff --git a/quantui/tddft_calc.py b/quantui/tddft_calc.py index 9104781..c0632a1 100644 --- a/quantui/tddft_calc.py +++ b/quantui/tddft_calc.py @@ -3,8 +3,15 @@ Computes vertical excitation energies and oscillator strengths using time-dependent density functional theory (TD-DFT). For Hartree-Fock -methods (RHF/UHF), falls back to TDHF (equivalent to CIS) and notes -this in the output. +methods (RHF/UHF), falls back to full TDHF (the random-phase +approximation, RPA — ``mf.TDHF()``) and notes this in the output. + +AUDIT additional-concerns — TDHF is NOT the same method as CIS. Full +TDHF/RPA includes the excitation/de-excitation (A/B block) coupling that +the Tamm-Dancoff approximation (TDA) drops; CIS is HF's TDA. This module +calls ``mf.TDHF()`` (full RPA), so its labels say TDHF/RPA rather than +CIS. See PySCF's own discussion of the distinction: +https://pyscf.org/user/tddft.html Platform notes -------------- @@ -56,7 +63,12 @@ class TDDFTResult: energy_hartree: Ground-state SCF energy in Hartrees. homo_lumo_gap_ev: HOMO-LUMO gap in eV from the ground-state SCF, or ``None``. - converged: ``True`` if the ground-state SCF converged. + converged: ``True`` only when BOTH the ground-state SCF converged + AND (if excited states were requested and the solve ran) every + requested TD root converged (AUDIT F08) — an SCF-only flag is + not overall success for a calculation whose deliverable is the + excited states. See ``td_converged``/``n_converged_states`` for + the per-root detail this folds together. n_iterations: Number of ground-state SCF macro-iterations. method: DFT functional or HF method used. basis: Basis set. @@ -78,6 +90,16 @@ class TDDFTResult: oscillator_strengths: List[float] = field(default_factory=list) nstates: int = 10 density_fit: bool = False + # AUDIT F08 — per-root Davidson convergence from the TD solver itself, + # distinct from the ground-state SCF's own converged flag above. None + # if the TD solve never ran (e.g. it raised before td.kernel() + # completed) or the installed PySCF doesn't expose td.converged. + td_converged: Optional[List[bool]] = None + # Count of roots whose td_converged flag is True — distinguishes + # "requested" (nstates), "returned" (len(excitation_energies_ev)), and + # "converged" root counts, since a Davidson solve can return energies + # for roots it never actually converged. + n_converged_states: Optional[int] = None # M-UX2 UXP2.10 — the actual PySCF class dispatched for the # ground-state SCF (e.g. "RHF", "UHF", "RKS", "UKS"); "" for an older # saved result. @@ -115,8 +137,10 @@ def run_tddft_calc( equations to compute the requested number of vertical excitation energies and their oscillator strengths. - When *method* is ``'RHF'`` or ``'UHF'``, the function uses TDHF (CIS) - rather than TD-DFT and writes a note to *progress_stream*. For a proper + When *method* is ``'RHF'`` or ``'UHF'``, the function uses full TDHF + (RPA — NOT the CIS/Tamm-Dancoff approximation; see the module + docstring) rather than TD-DFT, and writes a note to *progress_stream*. + For a proper UV-Vis simulation, a DFT functional such as ``'B3LYP'`` or ``'PBE0'`` is strongly recommended. @@ -234,7 +258,7 @@ def _run_tddft_calc_body( # M-UX2 UXP2.10 — capture before maybe_apply_d3 can wrap/rename it. scf_variant = type(mf).__name__ mf.xc = resolve_xc(method) - mf = maybe_apply_d3(mf, method, progress_stream=progress_stream) + mf, _ = maybe_apply_d3(mf, method, progress_stream=progress_stream) # Density fitting (RI), opt-in (M-DF). Off by default. TD-DFT is where the # measured win is largest (~1.6x on aspirin), so this is the primary target. @@ -245,9 +269,11 @@ def _run_tddft_calc_body( if using_hf and progress_stream is not None: try: progress_stream.write( - "\nNote: Using TDHF (CIS) for excited states — RHF/UHF was selected.\n" - "For a proper TD-DFT UV-Vis spectrum, use a DFT functional\n" - "such as B3LYP or PBE0 in the Method dropdown.\n\n" + "\nNote: Using TDHF/RPA for excited states — RHF/UHF was selected.\n" + "This is full TDHF (the random-phase approximation, with\n" + "excitation/de-excitation coupling), not the CIS/Tamm-Dancoff\n" + "approximation. For a proper TD-DFT UV-Vis spectrum, use a DFT\n" + "functional such as B3LYP or PBE0 in the Method dropdown.\n\n" ) except Exception: # noqa: BLE001 — cleanup (stream may be closed) pass @@ -297,11 +323,21 @@ def _run_tddft_calc_body( # ── TD-DFT / TDHF ──────────────────────────────────────────────────────── excitation_energies_ev: List[float] = [] oscillator_strengths: List[float] = [] + # AUDIT F08 — td.kernel() copies energies/oscillator strengths without + # checking td.converged (a per-root array from the Davidson solve). A + # real one-iteration-limited TDHF/6-31G water solve had + # converged=[False, False, False] yet reported success with three + # excitations. scf_converged is this function's SCF-only flag (the old + # sole source of `converged` below); td_converged/n_converged_states + # carry the TD solve's own per-root status. + scf_converged = converged + td_converged: Optional[List[bool]] = None + n_converged_states: Optional[int] = None try: emit_status( stream, - f"Solving {'TDHF (CIS)' if using_hf else 'TD-DFT'} " + f"Solving {'TDHF/RPA' if using_hf else 'TD-DFT'} " f"excited states ({nstates})…", ) td = mf.TDHF() if using_hf else mf.TDDFT() @@ -318,6 +354,11 @@ def _run_tddft_calc_body( osc = td.oscillator_strength() oscillator_strengths = [float(f) for f in osc] + _raw_td_converged = getattr(td, "converged", None) + if _raw_td_converged is not None: + td_converged = [bool(c) for c in _raw_td_converged] + n_converged_states = sum(td_converged) + except Exception as exc: logger.warning("TD-DFT/TDHF calculation failed: %s", exc) if progress_stream is not None: @@ -326,6 +367,18 @@ def _run_tddft_calc_body( except Exception: # noqa: BLE001 — cleanup (stream may be closed) pass + # AUDIT F08 — overall success requires every requested root to have + # actually converged, not just the ground-state SCF. A TD-DFT run whose + # entire purpose is the excited states is not "converged" if the + # Davidson solve raised before producing any roots, or if it returned + # roots that never converged. + converged = ( + scf_converged + and excitation_energies_ev != [] + and td_converged is not None + and all(td_converged) + ) + return TDDFTResult( energy_hartree=energy_hartree, homo_lumo_gap_ev=homo_lumo_gap_ev, @@ -339,4 +392,6 @@ def _run_tddft_calc_body( nstates=nstates, density_fit=density_fit_used, scf_variant=scf_variant, + td_converged=td_converged, + n_converged_states=n_converged_states, ) diff --git a/tests/test_app.py b/tests/test_app.py index 29b6d8e..bf8ec28 100644 --- a/tests/test_app.py +++ b/tests/test_app.py @@ -1265,6 +1265,42 @@ def test_solvent_shown_in_format_result(self): assert "Ethanol" in html assert "PCM" in html + @pytest.mark.parametrize( + "calc_type", ["Frequency", "UV-Vis (TD-DFT)", "NMR Shielding", "PES Scan"] + ) + def test_solvent_disabled_for_unsupported_calc_type(self, calc_type): + """AUDIT F11 — run_freq_calc/run_tddft_calc/run_nmr_calc/run_pes_scan + don't accept a solvent argument at all, so the checkbox must not be + left checkable (and checked) for these — that used to be a + silent no-op.""" + app = QuantUIApp() + app.solvent_cb.value = True + app.calc_type_dd.value = calc_type + assert app.solvent_cb.value is False + assert app.solvent_cb.disabled is True + # solvent_dd follows solvent_cb via the existing observer. + assert app.solvent_dd.layout.display == "none" + + @pytest.mark.parametrize( + "calc_type", ["Single Point", "Geometry Opt", "Reorganization Energy"] + ) + def test_solvent_enabled_for_supported_calc_type(self, calc_type): + app = QuantUIApp() + app.calc_type_dd.value = "Frequency" # disables it + app.calc_type_dd.value = calc_type # switching back must re-enable + assert app.solvent_cb.disabled is False + + def test_solvent_checkbox_re_enables_after_switching_back(self): + app = QuantUIApp() + app.solvent_cb.value = True + app.calc_type_dd.value = "PES Scan" + assert app.solvent_cb.disabled is True + app.calc_type_dd.value = "Single Point" + assert app.solvent_cb.disabled is False + # Re-enabling must not silently re-check it — the user unchecked + # nothing; the app did, and switching back doesn't restore intent. + assert app.solvent_cb.value is False + # --------------------------------------------------------------------------- # M-CAL — Calibration UI widgets @@ -1462,6 +1498,35 @@ def test_no_sto3g_warning_for_631g(self): html = app._format_nmr_result(self._make_nmr(basis="6-31G*")) assert "qualitative" not in html + def test_fallback_reference_shows_warning(self): + """AUDIT F17 — a fallback substitution must be visible on the card, + not just silently carried on the result object.""" + from quantui.nmr_calc import NMRResult + + app = QuantUIApp() + result = NMRResult( + atom_symbols=["O", "H", "H"], + shielding_iso_ppm=[320.1, 28.5, 28.5], + chemical_shifts_ppm={1: 3.22, 2: 3.22}, + method="CAM-B3LYP", + basis="6-31G*", + formula="H2O", + converged=True, + reference_key="B3LYP/6-31G*", + is_fallback_reference=True, + ) + html = app._format_nmr_result(result) + assert "B3LYP/6-31G*" in html + assert "⚠" in html + + def test_exact_match_reference_shows_no_warning(self): + result = self._make_nmr() + result.reference_key = "B3LYP/6-31G*" + result.is_fallback_reference = False + app = QuantUIApp() + html = app._format_nmr_result(result) + assert "⚠ No reference" not in html + def test_not_converged_shows_warning(self): app = QuantUIApp() html = app._format_nmr_result(self._make_nmr(converged=False)) @@ -1484,6 +1549,41 @@ def test_no_hc_atoms_shows_empty_message(self): assert "No ¹H or ¹³C" in html +class TestNmrSavePersistsFallbackReference: + """AUDIT F17 — the local NMR save path used to omit reference_key and + is_fallback_reference from the saved spectra, even though the backend + (nmr_calc.run_nmr_calc) already computes both and the batch NMR + serializer already includes them. + """ + + def test_local_nmr_save_includes_reference_metadata(self): + from quantui.nmr_calc import NMRResult + + app = QuantUIApp() + app._set_molecule(_water()) + app.calc_type_dd.value = "NMR Shielding" + mock_result = NMRResult( + atom_symbols=["O", "H", "H"], + shielding_iso_ppm=[320.1, 28.5, 28.5], + chemical_shifts_ppm={1: 3.22, 2: 3.22}, + method="CAM-B3LYP", + basis="6-31G*", + formula="H2O", + converged=True, + reference_key="B3LYP/6-31G*", + is_fallback_reference=True, + ) + with patch("quantui.nmr_calc.run_nmr_calc", return_value=mock_result): + with patch("quantui.save_result") as mock_save: + app._do_run() + + mock_save.assert_called_once() + _, kwargs = mock_save.call_args + nmr_spectra = kwargs["spectra"]["nmr"] + assert nmr_spectra["reference_key"] == "B3LYP/6-31G*" + assert nmr_spectra["is_fallback_reference"] is True + + # --------------------------------------------------------------------------- # M-IR — IR Spectrum accordion widgets # --------------------------------------------------------------------------- @@ -3020,6 +3120,64 @@ def _fake_generate(_atom, _basis, _coeff, _idx, out_path, **_kwargs): mock_gen.assert_called_once() mock_plot.assert_called_once() + def test_render_orbital_isosurface_uses_snapshotted_method_not_live_dropdown( + self, tmp_path + ): + """AUDIT additional-concerns — cube provenance used to read the + LIVE Method dropdown at Generate time, not necessarily what + actually produced the stored mo_coeff. Loads orbitals from a + result computed with method='B3LYP' (setting + app._last_orb_method via show_orbital_diagram), then changes the + dropdown to 'RHF' before generating — the cube's method label + must still say 'B3LYP', from the result's own immutable + provenance, not the now-mismatched dropdown. + """ + from unittest.mock import MagicMock + + import numpy as np + + app = QuantUIApp() + app._last_result_dir = tmp_path + + result = MagicMock() + result.formula = "H2O" + result.method = "B3LYP" + result.mo_energy_hartree = np.array([-1.5, -0.8, 0.2, 0.9]) + result.mo_occ = np.array([2.0, 2.0, 0.0, 0.0]) + result.mo_coeff = [[1.0, 0.0], [0.0, 1.0]] + result.pyscf_mol_atom = [["H", [0.0, 0.0, 0.0]]] + result.pyscf_mol_basis = "sto-3g" + app._show_orbital_diagram(result) + assert app._last_orb_method == "B3LYP" + + # Simulate the dropdown changing after the orbitals were loaded — + # e.g. the user tries a different method, or a different History + # result populated the Results panel without touching orbitals. + app.method_dd.value = "RHF" + + app._resolve_backend = lambda task: "plotlymol" + captured: dict[str, object] = {} + + def _fake_generate(_atom, _basis, _coeff, _idx, out_path, **kwargs): + captured["method"] = kwargs.get("method") + out_path.write_text("cube", encoding="utf-8") + return out_path + + with ( + patch( + "quantui.orbital_visualization.generate_cube_from_arrays", + side_effect=_fake_generate, + ), + patch( + "quantui.orbital_visualization.plot_cube_isosurface", + return_value=MagicMock(), + ), + patch("plotly.io.to_html", return_value="
iso
"), + ): + app._render_orbital_isosurface("HOMO") + + assert captured["method"] == "B3LYP" + def test_render_orbital_isosurface_py3dmol_path(self, tmp_path): # When the backend resolves to py3Dmol, the renderer is the py3Dmol # cube path (not Plotly), and the cube is still saved to disk. diff --git a/tests/test_app_formatters.py b/tests/test_app_formatters.py index 1a92cba..76b53e4 100644 --- a/tests/test_app_formatters.py +++ b/tests/test_app_formatters.py @@ -346,6 +346,78 @@ def test_format_past_result_renders_mp2_breakdown(): assert "MP2 correlation" in html +def test_format_past_result_restores_frequency_thermo(): + """AUDIT F18 — persisted thermochemistry must be restored into the + History card, not just silently saved and never shown again.""" + data = { + "calc_type": "frequency", + "converged": True, + "homo_lumo_gap_ev": 27.0, + "energy_hartree": -74.933498241, + "energy_ev": -2039.29, + "n_iterations": 12, + "timestamp": "2026-06-10_15-48-02-285574", + "formula": "H2O", + "method": "RHF", + "basis": "STO-3G", + "spectra": { + "ir": { + "frequencies_cm1": [1600.0, 3700.0, 3800.0], + "ir_intensities": [10.0, 5.0, 5.0], + "raman_activities": [], + "zpve_hartree": 0.021, + "thermo": { + "zpve_hartree": 0.021, + "H_hartree": -74.933498241, + "S_jmol": 188.538424, + "G_hartree": -74.954908540, + "temperature_k": 298.15, + "pressure_atm": 1.0, + "approximation": "ideal_gas_rigid_rotor_harmonic_oscillator", + }, + }, + "molecule": { + "atoms": ["O", "H", "H"], + "coords": [[0, 0, 0], [0.96, 0, 0], [0, 0.96, 0]], + "charge": 0, + "multiplicity": 1, + }, + }, + } + html = format_past_result(data) + assert "Thermochemistry at 298 K" in html + assert "-74.933498 Ha" in html # H + assert "188.54 J" in html # S + assert "-74.954909 Ha" in html # G + + +def test_format_past_result_frequency_without_thermo_omits_section(): + """A saved result predating this fix (or a thermo-less run) has no + 'thermo' key at all — that must be a silent no-op, not an error.""" + data = { + "calc_type": "frequency", + "converged": True, + "homo_lumo_gap_ev": 27.0, + "energy_hartree": -74.9, + "energy_ev": -2039.0, + "n_iterations": 12, + "timestamp": "2026-06-10_15-48-02-285574", + "formula": "H2O", + "method": "RHF", + "basis": "STO-3G", + "spectra": { + "ir": { + "frequencies_cm1": [], + "ir_intensities": [], + "raman_activities": [], + "zpve_hartree": 0.0, + }, + }, + } + html = format_past_result(data) + assert "Thermochemistry" not in html + + def test_format_past_result_hf_dft_has_no_breakdown(): data = { "calc_type": "single_point", diff --git a/tests/test_backends_slurm.py b/tests/test_backends_slurm.py index 920c05c..26c2f36 100644 --- a/tests/test_backends_slurm.py +++ b/tests/test_backends_slurm.py @@ -3,6 +3,7 @@ """ import json +import shlex import subprocess import sys import time @@ -165,6 +166,119 @@ def test_dispatch_records_submit_timestamp(self, mock_run, slurm_backend): assert since < 5 +class TestSlurmBackendWorkerCommandQuoting: + """AUDIT F21 — the generated worker command is embedded as shell text + in the sbatch script (not passed as an argv list), so an unquoted path + containing a space splits into multiple shell arguments. The audit's + own reproduction: ``/tmp/audit folder/request.json`` became + ``--request /tmp/audit`` plus a stray ``folder/request.json`` token. + """ + + def test_request_path_with_space_stays_one_argument(self, slurm_backend): + request_path = Path("/tmp/audit folder/request.json") + staging_dir = Path("/tmp/audit folder/staging") + cmd = slurm_backend._worker_command(request_path, staging_dir) + tokens = shlex.split(cmd) + assert tokens[-2] == "--request" + assert tokens[-1] == str(request_path) + assert ( + len(tokens) == tokens.index("--request") + 2 + ), f"request path split into extra shell tokens: {tokens}" + + def test_apptainer_branch_quotes_staging_and_request_paths(self, tmp_path): + registry = JobRegistry( + jobs_root=tmp_path / "jobs", staging_root=tmp_path / "staging" + ) + backend = SlurmBackend( + registry=registry, + partition="test", + use_apptainer=True, + apptainer_image="/opt/images/quantui image.sif", + ) + request_path = Path("/tmp/audit folder/request.json") + staging_dir = Path("/tmp/audit folder/staging") + cmd = backend._worker_command(request_path, staging_dir) + tokens = shlex.split(cmd) + assert str(staging_dir) in tokens + assert "/opt/images/quantui image.sif" in tokens + assert tokens[-1] == str(request_path) + + def test_dispatch_with_space_in_staging_path_produces_parseable_script( + self, mock_slurm_env, tmp_path + ): + """End-to-end: a staging root containing a space must still produce + a submit.slurm whose worker-command line survives a real shell + tokenization pass with the request path intact as one argument.""" + spaced_root = tmp_path / "audit folder" + registry = JobRegistry( + jobs_root=spaced_root / "jobs", + staging_root=spaced_root / "staging", + ) + backend = SlurmBackend(registry=registry, partition="test", use_apptainer=False) + + with patch( + "quantui.backends.slurm.subprocess.run", + **{ + "return_value.stdout": "Submitted batch job 777888\n", + "return_value.stderr": "", + "return_value.returncode": 0, + }, + ): + rid = backend.dispatch(_request("spaced001")) + + record = backend.registry.load(rid) + slurm_text = (record.staging_path / "submit.slurm").read_text() + worker_line = next( + line + for line in slurm_text.splitlines() + if "quantui.backends.worker" in line + ) + tokens = shlex.split(worker_line) + request_arg = tokens[tokens.index("--request") + 1] + assert request_arg == str(record.staging_path / "request.json") + assert Path(request_arg).exists() + + def test_output_and_error_directives_are_quoted_for_spaced_paths( + self, mock_slurm_env, tmp_path + ): + """SBATCH directive lines are also parsed word-by-word by sbatch; + an unquoted --output/--error path with a space in it is the same + class of bug as the worker command, just one line up.""" + spaced_root = tmp_path / "audit folder" + registry = JobRegistry( + jobs_root=spaced_root / "jobs", + staging_root=spaced_root / "staging", + ) + backend = SlurmBackend(registry=registry, partition="test", use_apptainer=False) + + with patch( + "quantui.backends.slurm.subprocess.run", + **{ + "return_value.stdout": "Submitted batch job 777999\n", + "return_value.stderr": "", + "return_value.returncode": 0, + }, + ): + rid = backend.dispatch(_request("spaced002")) + + record = backend.registry.load(rid) + slurm_text = (record.staging_path / "submit.slurm").read_text() + output_line = next( + line + for line in slurm_text.splitlines() + if line.startswith("#SBATCH --output=") + ) + error_line = next( + line + for line in slurm_text.splitlines() + if line.startswith("#SBATCH --error=") + ) + assert ( + output_line == f'#SBATCH --output="{record.staging_path / "slurm-%j.out"}"' + ) + assert error_line == f'#SBATCH --error="{record.staging_path / "slurm-%j.err"}"' + + class TestSlurmBackendReconcile: def test_reconcile_stale_record_without_slurm_id(self, slurm_backend, monkeypatch): monkeypatch.setenv("QUANTUI_SLURM_STALE_NO_ID_S", "60") diff --git a/tests/test_backends_worker.py b/tests/test_backends_worker.py index 72ad658..1ac2576 100644 --- a/tests/test_backends_worker.py +++ b/tests/test_backends_worker.py @@ -50,6 +50,46 @@ def test_unsupported_calc_type_returns_error(self, staging): assert outcome.status == "error" assert outcome.error["code"] == "UNSUPPORTED_CAPABILITY" + @pytest.mark.parametrize( + "calc_type", ["geometry_opt", "frequency", "tddft", "nmr", "pes_scan"] + ) + def test_solvent_on_unsupported_calc_type_returns_error(self, staging, calc_type): + """AUDIT F11 — run_freq_calc/run_tddft_calc/run_nmr_calc/run_pes_scan/ + optimize_geometry don't accept a solvent argument at all; a + solvent set for one of these calc_types must fail the request + rather than silently run gas-phase. + """ + data = json.loads((staging / "request.json").read_text()) + data["calc_type"] = calc_type + data["solvent"] = "water" + (staging / "request.json").write_text(json.dumps(data)) + + outcome = run_worker_request(staging / "request.json") + assert outcome.status == "error" + assert outcome.error["code"] == "UNSUPPORTED_CAPABILITY" + assert "solvent" in outcome.error["user_message"].lower() + + @patch("quantui.session_calc.run_in_session") + def test_solvent_on_single_point_is_accepted(self, mock_run, staging): + """Sanity check: the calc_types that DO support solvent must not be + rejected by the new guard.""" + data = json.loads((staging / "request.json").read_text()) + data["solvent"] = "water" + (staging / "request.json").write_text(json.dumps(data)) + mock_run.return_value = SimpleNamespace( + energy_hartree=-1.12, + homo_lumo_gap_ev=10.0, + converged=True, + n_iterations=5, + method="RHF", + basis="STO-3G", + formula="H2", + ) + + outcome = run_worker_request(staging / "request.json") + assert outcome.status == "success" + assert mock_run.call_args.kwargs["solvent"] == "water" + @patch("quantui.session_calc.run_in_session") def test_single_point_success(self, mock_run, staging): mock_run.return_value = SimpleNamespace( @@ -228,7 +268,7 @@ class TestCheckpointWiring: most for a killed run at, say, point 20 of 25. """ - def _identity(self, *, calc_type: str): + def _identity(self, *, calc_type: str, extra: tuple = ()): from quantui.checkpoint import CalcIdentity from quantui.molecule import Molecule @@ -239,7 +279,7 @@ def _identity(self, *, calc_type: str): multiplicity=1, ) return CalcIdentity.from_molecule( - mol, calc_type=calc_type, method="RHF", basis="STO-3G" + mol, calc_type=calc_type, method="RHF", basis="STO-3G", extra=extra ) @patch("quantui.optimizer.optimize_geometry") @@ -368,7 +408,11 @@ def test_pes_scan_resumes_and_reports_points_already_computed( ): import json as _json - ckpt = self._identity(calc_type="pes_scan") + # AUDIT F10 — the request sets no scan_type/atom_indices options, so + # _run_pes_scan defaults to scan_type="bond", atom_indices=[0, 1]; + # the checkpoint identity here must match that exactly, since + # resume_key now includes them. + ckpt = self._identity(calc_type="pes_scan", extra=("bond", "0", "1")) from quantui.checkpoint import Checkpoint real_ckpt = Checkpoint(ckpt, root=staging / ".checkpoint") diff --git a/tests/test_calculator.py b/tests/test_calculator.py index 42c7d59..8f74638 100644 --- a/tests/test_calculator.py +++ b/tests/test_calculator.py @@ -185,6 +185,200 @@ def test_script_charged_molecule(self, tmp_path): assert "charge = 1" in script_content + def test_ecp_embedded_for_heavy_element_basis(self, tmp_path): + """The resolved ECP mapping must appear in the generated script for + a basis that carries one (LANL2DZ on Na), and be an explicit empty + dict for an all-electron basis — never simply absent. + + Requires PySCF (ecp_for_basis() looks up LANL2DZ's ECP table via + pyscf.gto.basis.load_ecp) — unlike script *generation* itself, + which must keep working without it (see + test_script_generation_works_without_pyscf_installed above); on a + machine without PySCF the mapping correctly degrades to {}, which + this test cannot verify one way or the other. + """ + pytest.importorskip("pyscf") + na_h = Molecule(["Na", "H"], [[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]) + calc = PySCFCalculation(na_h, method="RHF", basis="LANL2DZ") + script_content = calc.generate_calculation_script(tmp_path / "nah.py") + assert "mol.ecp = {'Na': 'LANL2DZ'}" in script_content + + water = Molecule( + ["O", "H", "H"], + [[0.0, 0.0, 0.0], [0.757, 0.587, 0.0], [-0.757, 0.587, 0.0]], + ) + calc_ae = PySCFCalculation(water, method="RHF", basis="6-31G") + script_ae = calc_ae.generate_calculation_script(tmp_path / "water.py") + assert "mol.ecp = {}" in script_ae + + def test_script_generation_works_without_pyscf_installed( + self, tmp_path, monkeypatch + ): + """CI regression — AUDIT F06's ecp_for_basis() imports pyscf, but + generate_calculation_script() must keep working on a machine that + doesn't have PySCF installed at all (e.g. Windows, no WSL): the + module's own docstring says the exported script is meant to be + downloaded and run independently, possibly on a different machine + than the one that generated it. This broke Windows CI outright + (ModuleNotFoundError at *generation* time, not just execution) + until this fallback was added — reproduced here by making the + import raise ImportError regardless of platform. + """ + import builtins + + real_import = builtins.__import__ + + def _no_pyscf(name, *args, **kwargs): + if name == "pyscf" or name.startswith("pyscf."): + raise ImportError("simulated: PySCF not installed") + return real_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", _no_pyscf) + + na_h = Molecule(["Na", "H"], [[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]) + calc = PySCFCalculation(na_h, method="RHF", basis="LANL2DZ") + script_content = calc.generate_calculation_script(tmp_path / "nah.py") + # Degrades to the pre-AUDIT-F06 (all-electron) mapping on this + # platform only — script generation itself must not raise. + assert "mol.ecp = {}" in script_content + assert (tmp_path / "nah.py").exists() + + @pytest.mark.slow + def test_exported_ecp_script_reproduces_in_app_electron_count(self, tmp_path): + """AUDIT F06 regression — executes the exported NaH/LANL2DZ script + (real PySCF subprocess, no QuantUI import) and confirms it reports + the correct 2-explicit-electron ECP calculation, not the wrong + 12-electron all-electron one the unfixed template produced. + """ + pytest.importorskip("pyscf") + import subprocess + import sys + + na_h = Molecule(["Na", "H"], [[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]) + calc = PySCFCalculation(na_h, method="RHF", basis="LANL2DZ") + script_path = tmp_path / "nah_rhf_lanl2dz.py" + calc.generate_calculation_script(script_path) + + proc = subprocess.run( + [sys.executable, str(script_path)], + cwd=tmp_path, + capture_output=True, + text=True, + timeout=120, + ) + assert proc.returncode == 0, proc.stdout + proc.stderr + assert "Number of electrons: 2" in proc.stdout + # The unfixed template reported 12 electrons and E ~= -20.13 Ha for + # this system; the correct ECP calculation converges near -0.7 Ha. + assert "Number of electrons: 12" not in proc.stdout + assert (tmp_path / "results.npz").exists() + + @pytest.mark.slow + def test_exported_npz_has_fields_the_cube_helper_needs(self, tmp_path): + """AUDIT additional-concerns — generate_cube_file() requires + 'mol_atom'/'mol_basis' (raises ValueError without them: "Re-run + the calculation with the updated script template" — a template + that never actually wrote them) and reads an optional 'mo_occ' to + infer charge/spin. The exported script used to save only energy/ + mo_energy/mo_coeff/converged, so a cube could never be generated + from a standalone-exported result. + """ + pytest.importorskip("pyscf") + import subprocess + import sys + + import numpy as np + + water = Molecule( + ["O", "H", "H"], + [[0.0, 0.0, 0.0], [0.757, 0.587, 0.0], [-0.757, 0.587, 0.0]], + ) + calc = PySCFCalculation(water, method="RHF", basis="STO-3G") + script_path = tmp_path / "water_rhf.py" + calc.generate_calculation_script(script_path) + + proc = subprocess.run( + [sys.executable, str(script_path)], + cwd=tmp_path, + capture_output=True, + text=True, + timeout=120, + ) + assert proc.returncode == 0, proc.stdout + proc.stderr + + npz = np.load(tmp_path / "results.npz", allow_pickle=True) + assert "mol_atom" in npz.files + assert "mol_basis" in npz.files + assert "mo_occ" in npz.files + assert str(npz["mol_basis"]) == "STO-3G" + assert "O" in str(npz["mol_atom"]) + + # The actual regression: generate_cube_file must not raise with + # this exported npz — it used to unconditionally raise ValueError + # ("does not contain 'mol_atom'/'mol_basis' keys") on any + # standalone export. + from quantui.orbital_visualization import generate_cube_file + + cube_path = generate_cube_file( + tmp_path / "results.npz", + 0, + tmp_path / "orbital0.cube", + nx=6, + ny=6, + nz=6, + ) + assert cube_path.exists() + + @pytest.mark.slow + @pytest.mark.parametrize("method", ["MP2", "CCSD", "CCSD(T)"]) + def test_exported_post_hf_script_runs_successfully(self, tmp_path, method): + """AUDIT F13 regression — the exported script used to fall into the + DFT branch for any non-RHF/UHF method, setting + mf.xc = 'MP2'/'CCSD'/'CCSD(T)' and failing with + "LibXCFunctional: name '...' not found". Executes the real + generated script (subprocess, no QuantUI import) for water/STO-3G. + """ + pytest.importorskip("pyscf") + import subprocess + import sys + + water = Molecule( + ["O", "H", "H"], + [[0.0, 0.0, 0.0], [0.757, 0.587, 0.0], [-0.757, 0.587, 0.0]], + ) + calc = PySCFCalculation(water, method=method, basis="STO-3G") + script_path = tmp_path / "water_post_hf.py" + calc.generate_calculation_script(script_path) + + proc = subprocess.run( + [sys.executable, str(script_path)], + cwd=tmp_path, + capture_output=True, + text=True, + timeout=120, + ) + # Note: PySCF's verbose=4 logging echoes the script's own source + # (including this docstring) into stdout, so don't substring-match + # error text there — exit code 0 plus the numeric checks below are + # the real assertion. + assert proc.returncode == 0, proc.stdout + proc.stderr + assert "Total energy" in proc.stdout + assert (tmp_path / "results.npz").exists() + + import numpy as np + + npz = np.load(tmp_path / "results.npz", allow_pickle=True) + # RHF/STO-3G water energy is ~-74.963; any real correlation energy + # must make the total more negative than that. + assert float(npz["energy"]) < -74.97 + if method == "MP2": + assert float(npz["mp2_correlation_hartree"]) < 0 + else: + assert float(npz["ccsd_correlation_hartree"]) < 0 + assert bool(npz["cc_converged"]) is True + if method == "CCSD(T)": + assert float(npz["ccsd_t_correction_hartree"]) < 0 + def test_script_creates_parent_directories(self, tmp_path): """Test that script creation makes parent directories.""" script_path = tmp_path / "nested" / "dir" / "calc.py" diff --git a/tests/test_checkpoint_wiring.py b/tests/test_checkpoint_wiring.py index 9551814..64e2575 100644 --- a/tests/test_checkpoint_wiring.py +++ b/tests/test_checkpoint_wiring.py @@ -104,17 +104,62 @@ def test_warm_start_is_on_by_default(self): ) def test_app_passes_the_checkpoint_to_every_long_calc(self): - """The three calc types that can be interrupted must all receive one.""" + """The calc types that can be interrupted must all receive one.""" import quantui.app as A src = Path(A.__file__).read_text(encoding="utf-8") - assert src.count("checkpoint=_ckpt") >= 3 + assert src.count("checkpoint=_ckpt") >= 4 def test_app_passes_resume_to_the_resumable_calc_types(self): import quantui.app as A src = Path(A.__file__).read_text(encoding="utf-8") - assert src.count("resume=_resume") >= 2 + assert src.count("resume=_resume") >= 3 + + def test_frequency_run_receives_the_actual_checkpoint_and_resume(self): + """AUDIT F16 — a source-text occurrence count can pass even if + Frequency's own call site never grew ``checkpoint=_ckpt`` — the + app opened ``_ckpt``/resolved ``_resume`` but never passed either + into ``run_freq_calc``. Drives the real dispatch (mocking only + ``run_freq_calc`` itself) and inspects the actual call arguments. + """ + from unittest.mock import patch + + from quantui.app import QuantUIApp + from quantui.freq_calc import FreqResult + from quantui.molecule import Molecule + + app = QuantUIApp() + app._set_molecule( + Molecule( + ["O", "H", "H"], + [[0.0, 0.0, 0.0], [0.757, 0.587, 0.0], [-0.757, 0.587, 0.0]], + ) + ) + app.calc_type_dd.value = "Frequency" + mock_result = FreqResult( + energy_hartree=-76.0, + homo_lumo_gap_ev=10.0, + converged=True, + n_iterations=8, + method="RHF", + basis="STO-3G", + formula="H2O", + frequencies_cm1=[1600.0, 3600.0, 3800.0], + ir_intensities=[1.0, 2.0, 3.0], + raman_activities=[], + zpve_hartree=0.02, + ) + with patch( + "quantui.freq_calc.run_freq_calc", return_value=mock_result + ) as mock_run: + with patch("quantui.save_result"): + app._do_run() + + mock_run.assert_called_once() + _, kwargs = mock_run.call_args + assert kwargs.get("checkpoint") is not None + assert "resume" in kwargs # ══ Warm-start selection ═════════════════════════════════════════════════════ @@ -222,16 +267,70 @@ def from_chk(self, path): class TestReuseScanPoint: class _FakeAtoms: + """Minimal ASE-Atoms-alike: implements exactly what + _reuse_scan_point needs, including the real geometry measurements + (AUDIT F10's cross-check) computed from whatever positions were + last set — not hardcoded, so a mismatched geometry is actually + caught rather than trivially passing. + """ + def __init__(self): self.positions = None def set_positions(self, coords): self.positions = coords + def _vec(self, i, j): + a, b = self.positions[i], self.positions[j] + return [b[k] - a[k] for k in range(3)] + + def get_distance(self, i1, i2): + import math + + v = self._vec(i1, i2) + return math.sqrt(sum(c * c for c in v)) + + def get_angle(self, i1, i2, i3): + import math + + v1 = self._vec(i2, i1) + v2 = self._vec(i2, i3) + dot = sum(a * b for a, b in zip(v1, v2)) + n1 = math.sqrt(sum(c * c for c in v1)) + n2 = math.sqrt(sum(c * c for c in v2)) + return math.degrees(math.acos(max(-1.0, min(1.0, dot / (n1 * n2))))) + + def get_dihedral(self, i1, i2, i3, i4): + import math + + p = self.positions + b1 = [p[i2][k] - p[i1][k] for k in range(3)] + b2 = [p[i3][k] - p[i2][k] for k in range(3)] + b3 = [p[i4][k] - p[i3][k] for k in range(3)] + + def cross(a, b): + return [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ] + + def norm(v): + return math.sqrt(sum(c * c for c in v)) + + n1 = cross(b1, b2) + n2 = cross(b2, b3) + m1 = cross(n1, [c / norm(b2) for c in b2]) + x = sum(a * b for a, b in zip(n1, n2)) + y = sum(a * b for a, b in zip(m1, n2)) + return math.degrees(math.atan2(y, x)) % 360.0 + def _record(self, **overrides): base = { "index": 1, "value": 1.25, + "scan_type": "bond", + "atom_indices": [0, 1], "energy_hartree": -1.5, "ok": True, "atoms": ["H", "H"], @@ -240,9 +339,11 @@ def _record(self, **overrides): base.update(overrides) return base - def _call(self, record, value=1.25): + def _call(self, record, value=1.25, scan_type="bond", atom_indices=(0, 1)): atoms = self._FakeAtoms() - result = pes_scan._reuse_scan_point(record, value, atoms, _FakeMolecule()) + result = pes_scan._reuse_scan_point( + record, value, atoms, _FakeMolecule(), scan_type, atom_indices + ) return result, atoms def test_reuses_a_matching_point(self): @@ -281,10 +382,64 @@ def test_rejects_an_empty_geometry(self): def test_carries_charge_and_multiplicity_from_the_live_molecule(self): atoms = self._FakeAtoms() molecule = _FakeMolecule(charge=-1, multiplicity=2) - result = pes_scan._reuse_scan_point(self._record(), 1.25, atoms, molecule) + result = pes_scan._reuse_scan_point( + self._record(), 1.25, atoms, molecule, "bond", (0, 1) + ) assert result[1].charge == -1 assert result[1].multiplicity == 2 + # ── AUDIT F10 — coordinate identity, not just the scalar value ────────── + + def test_rejects_a_different_scan_type_even_with_matching_value(self): + """The audit's exact scenario: a banked bond-scan point at 1.0 A + must not be reused for an angle scan whose target also happens to + be 1.0 (degrees, but the record's "value" field carries no unit).""" + record = self._record(value=1.0, scan_type="bond") + result, _ = self._call(record, value=1.0, scan_type="angle") + assert result is None + + def test_rejects_different_atom_indices_even_with_matching_scan_type(self): + """A bond scan of atoms (0,1) must not be reused for a bond scan of + atoms (0,2) — same scan_type, different coordinate.""" + record = self._record(atom_indices=[0, 1]) + result, _ = self._call(record, atom_indices=(0, 2)) + assert result is None + + def test_rejects_a_stored_geometry_whose_actual_coordinate_disagrees(self): + """Defense in depth: even if scan_type/atom_indices both match, the + stored geometry's OWN bond length must actually equal the target — + a record whose self-reported "value" doesn't match its real + coordinates (e.g. from a bug, or hand-edited checkpoint) must not + be trusted just because the label says so.""" + # atoms 0,1 are really 2.0 A apart, but the record claims 1.0 A. + record = self._record(value=1.0, coordinates=[[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]) + result, _ = self._call(record, value=1.0) + assert result is None + + def test_angle_scan_reuse_still_works(self): + """Sanity check that the new cross-check doesn't break a genuine + angle-scan reuse.""" + # H-O-H at 104.5 degrees, O at origin. + import math + + theta = math.radians(104.5) + coords = [ + [math.sin(theta / 2), math.cos(theta / 2), 0.0], + [0.0, 0.0, 0.0], + [-math.sin(theta / 2), math.cos(theta / 2), 0.0], + ] + record = self._record( + value=104.5, + scan_type="angle", + atom_indices=[0, 1, 2], + atoms=["H", "O", "H"], + coordinates=coords, + ) + result, _ = self._call( + record, value=104.5, scan_type="angle", atom_indices=(0, 1, 2) + ) + assert result is not None + # ══ Resume offer ═════════════════════════════════════════════════════════════ diff --git a/tests/test_freq_calc.py b/tests/test_freq_calc.py index e278a2c..fe958d1 100644 --- a/tests/test_freq_calc.py +++ b/tests/test_freq_calc.py @@ -176,6 +176,31 @@ def test_thermo_s_positive(self): if result.thermo is not None: assert result.thermo.S_jmol > 0 + @pyscf_only + @pytest.mark.slow + def test_thermo_matches_independent_pyscf_reference(self): + """AUDIT F01 regression — S_jmol/G_hartree against an independent + PySCF reference (not merely S > 0 / G < H), for RHF/STO-3G water at + the fixed geometry in ``_water()``. + + Before the F01 fix, PySCF's S_tot (returned in Eh/K) was stored + directly as S_jmol without converting to J/(mol*K), then divided by + _HARTREE_TO_JMOL a second time when forming G — deflating S_jmol by + ~2.6e6x and leaving G ~= H. Reference values below (S=188.538424 + J/(mol*K), G=-74.954908540 Eh) come from calling + pyscf.hessian.thermo.thermo() directly on the same RHF/STO-3G water + SCF object/frequencies, independent of quantui.freq_calc. + """ + from quantui.freq_calc import run_freq_calc + + result = run_freq_calc(_water(), method="RHF", basis="STO-3G") + assert result.thermo is not None + assert result.thermo.S_jmol == pytest.approx(188.538424, abs=0.01) + assert result.thermo.G_hartree == pytest.approx(-74.954908540, abs=1e-6) + # The old bug's error was ~8e-9 Eh (S_jmol deflated to ~7.18e-5); a + # correct calculation differs from H by orders of magnitude more. + assert result.thermo.H_hartree - result.thermo.G_hartree > 1e-3 + @pyscf_only @pytest.mark.slow def test_thermo_g_less_than_h(self): @@ -239,6 +264,37 @@ def test_post_hf_method_raises_value_error(self, method): run_freq_calc(_water(), method=method, basis="STO-3G") +# ============================================================================ +# AUDIT F15 — a failed Hessian must not read as a converged result +# ============================================================================ + + +class TestFreqResultReflectsHessianCompletion: + """A ROHF reference's analytic Hessian is unavailable on this path + (PySCF has no ROHF Hessian implementation here); the caught exception + used to leave FreqResult.converged reading whatever the reference SCF + alone reported, with frequencies_cm1=[] — a frequency calculation with + no computed Hessian is not a successful frequency analysis. + """ + + @pyscf_only + @pytest.mark.slow + def test_rohf_hessian_failure_reports_unconverged(self): + from quantui.freq_calc import run_freq_calc + from quantui.molecule import Molecule + + # Real OH doublet — RHF/STO-3G dispatches to ROHF for this + # open-shell molecule, matching the audit's exact reproduction. + oh = Molecule( + ["O", "H"], [[0.0, 0.0, 0.0], [0.0, 0.0, 0.97]], charge=0, multiplicity=2 + ) + result = run_freq_calc(oh, method="RHF", basis="STO-3G") + + assert result.scf_variant == "ROHF" + assert result.frequencies_cm1 == [] + assert result.converged is False + + # ============================================================================ # IR intensities — PySCF required # ============================================================================ @@ -355,6 +411,215 @@ def test_freq_ir_workers_dispatches_on_dm0_shape_not_spin(self): finally: os.unlink(tmp.name) + def test_ecp_omission_gives_a_different_hamiltonian(self): + """AUDIT F05 regression — the worker must run the reference's ECP + (e.g. LANL2DZ on Na), not silently fall back to all-electron. + + Without ``ecp``, NaH/LANL2DZ is an all-electron (12-electron) + calculation instead of the correct 2-explicit-electron ECP one — + a different Hamiltonian, not numerical noise. Confirms both the + electron-count claim directly (via the same ecp_for_basis mapping + the worker now receives) and that the worker's own SCF result + (the dipole it returns) differs materially between the two cases. + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto + + from quantui.freq_ir_workers import init_worker, run_displaced_scf + from quantui.inorganic_guards import ecp_for_basis + + atom_str = "Na 0 0 0; H 0 0 2.0" + basis = "LANL2DZ" + ecp = ecp_for_basis(basis, ["Na", "H"]) + assert ecp == {"Na": "LANL2DZ"} + + mol_with_ecp = gto.M( + atom=atom_str, basis=basis, ecp=ecp, charge=0, spin=0, verbose=0 + ) + mol_without_ecp = gto.M( + atom=atom_str, basis=basis, ecp={}, charge=0, spin=0, verbose=0 + ) + assert mol_with_ecp.nelectron == 2 + assert mol_without_ecp.nelectron == 12 + + coords = mol_with_ecp.atom_coords(unit="Bohr").flatten().tolist() + + def _run(ecp_arg): + tmp = tempfile.NamedTemporaryFile(delete=False, suffix=".pkl") + try: + pickle.dump(None, tmp) + tmp.close() + init_worker(atom_str, basis, 0, 0, None, tmp.name, 1, None, ecp_arg) + return run_displaced_scf("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + dip_with_ecp = _run(ecp) + dip_without_ecp = _run({}) + # A 2-electron vs 12-electron calculation on the same geometry + # produces a substantially different dipole, not a small + # numerical discrepancy. + assert abs(dip_with_ecp[2] - dip_without_ecp[2]) > 1.0 + + def test_density_fit_option_applied_to_displaced_scf(self, monkeypatch): + """AUDIT F19 regression — the worker had no density_fit parameter at + all; every displaced SCF ran without density fitting even when the + reference (and the serial fallback in freq_calc.py, which calls + ``_try_density_fit(_mf_d, enabled=_density_fit_used)``) used it — + silently changing the numerical approximation for parallel IR runs + on a fitted reference, not merely its speed. + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto, scf + + import quantui.density_fitting as density_fitting + 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.RHF(mol) + mf.kernel() + dm0 = mf.make_rdm1() + + seen_enabled = [] + _orig = density_fitting.try_density_fit + + def _spy(mf_arg, *, enabled=None, auxbasis=None): + seen_enabled.append(enabled) + return _orig(mf_arg, enabled=enabled, auxbasis=auxbasis) + + monkeypatch.setattr(density_fitting, "try_density_fit", _spy) + + 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, + None, # checkpoint_items_dir + None, # ecp + True, # density_fit + True, # scf_rescue + ) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + dip = run_displaced_scf("d000_x_+", coords) + assert len(dip) == 3 + finally: + os.unlink(tmp.name) + + assert seen_enabled == [True] + + def test_density_fit_defaults_off_for_an_older_caller(self, monkeypatch): + """Backward compatibility: a caller that predates this fix (only + positional args through ``ecp``) must still get density_fit=False, + matching the old always-off behavior exactly.""" + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto, scf + + import quantui.density_fitting as density_fitting + 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.RHF(mol) + mf.kernel() + dm0 = mf.make_rdm1() + + seen_enabled = [] + _orig = density_fitting.try_density_fit + + def _spy(mf_arg, *, enabled=None, auxbasis=None): + seen_enabled.append(enabled) + return _orig(mf_arg, enabled=enabled, auxbasis=auxbasis) + + monkeypatch.setattr(density_fitting, "try_density_fit", _spy) + + 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() + run_displaced_scf("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + assert seen_enabled == [False] + + def test_scf_rescue_option_honored_in_displaced_scf(self, monkeypatch): + """AUDIT F19 regression — the worker always called + ``run_scf_with_rescue(mf, dm0=dm0)``, taking its default + ``rescue=True`` regardless of what the caller requested. The + serial loop threads the caller's choice through as + ``rescue=scf_rescue``; this confirms the worker now does too. + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto, scf + + import quantui.scf_robust as scf_robust + 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.RHF(mol) + mf.kernel() + dm0 = mf.make_rdm1() + + seen_rescue = [] + _orig = scf_robust.run_scf_with_rescue + + def _spy(mf_arg, **kwargs): + seen_rescue.append(kwargs.get("rescue", True)) + return _orig(mf_arg, **kwargs) + + monkeypatch.setattr(scf_robust, "run_scf_with_rescue", _spy) + + 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, + None, # checkpoint_items_dir + None, # ecp + False, # density_fit + False, # scf_rescue + ) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + run_displaced_scf("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + assert seen_rescue == [False] + def test_worker_writes_its_own_checkpoint_record(self, tmp_path): """M-CHECKPOINT CHK.4.4's crash-safety claim, at the unit level: the worker itself durably records completion — not just the parent diff --git a/tests/test_freq_ir_workers.py b/tests/test_freq_ir_workers.py index 5aa55e0..3ef2aaa 100644 --- a/tests/test_freq_ir_workers.py +++ b/tests/test_freq_ir_workers.py @@ -17,6 +17,7 @@ from quantui.freq_ir_workers import ( _truthy, + available_cpu_count, freq_parallel_env_configured, freq_parallel_opt_in, parallel_enabled_for_run, @@ -87,6 +88,41 @@ def test_env_var_overrides_settings_off(self, monkeypatch, tmp_path): assert parallel_enabled_for_run(cpu_count=16, displacement_count=60) is False +class TestAvailableCpuCount: + """AUDIT additional-concerns — worker sizing must respect a SLURM + allocation/cgroup limit, not just report the whole host's core count. + """ + + def test_slurm_cpus_per_task_takes_precedence(self, monkeypatch): + monkeypatch.setenv("SLURM_CPUS_PER_TASK", "8") + assert available_cpu_count() == 8 + + def test_ignores_invalid_slurm_value(self, monkeypatch): + monkeypatch.setenv("SLURM_CPUS_PER_TASK", "not-a-number") + # Falls through to sched_getaffinity/cpu_count — just must not raise + # and must return a positive int. + assert available_cpu_count() >= 1 + + def test_ignores_non_positive_slurm_value(self, monkeypatch): + monkeypatch.setenv("SLURM_CPUS_PER_TASK", "0") + assert available_cpu_count() >= 1 + + def test_falls_back_to_affinity_or_cpu_count_without_slurm(self, monkeypatch): + import os + + monkeypatch.delenv("SLURM_CPUS_PER_TASK", raising=False) + result = available_cpu_count() + try: + expected = len(os.sched_getaffinity(0)) + except (AttributeError, OSError, NotImplementedError): + expected = os.cpu_count() or 1 + assert result == expected + + def test_never_returns_less_than_one(self, monkeypatch): + monkeypatch.delenv("SLURM_CPUS_PER_TASK", raising=False) + assert available_cpu_count() >= 1 + + class TestPickWorkerCount: """Worker count = ``min(cpu // 2, displacement_count)``, floored at 1.""" diff --git a/tests/test_freq_raman_workers.py b/tests/test_freq_raman_workers.py index 5cb0259..c09233b 100644 --- a/tests/test_freq_raman_workers.py +++ b/tests/test_freq_raman_workers.py @@ -2,6 +2,8 @@ from __future__ import annotations +import pytest + from quantui import freq_ir_workers as irw @@ -16,3 +18,160 @@ def test_raman_calc_imports_raman_workers(self): assert callable(rw.init_raman_worker) assert callable(rw.run_displaced_polarizability) + + +class TestRamanWorkerEcp: + def test_ecp_omission_gives_a_different_hamiltonian(self): + """AUDIT F05 regression — mirrors the IR worker's ECP fix + (test_freq_calc.py::TestIrIntensityUhfClosedShellDispatch), for the + Raman polarizability worker. Without ``ecp``, NaH/LANL2DZ runs as + an all-electron (12-electron) calculation instead of the correct + 2-explicit-electron ECP one, giving a materially different + polarizability tensor. + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto + + from quantui.freq_raman_workers import ( + init_raman_worker, + run_displaced_polarizability, + ) + from quantui.inorganic_guards import ecp_for_basis + + atom_str = "Na 0 0 0; H 0 0 2.0" + basis = "LANL2DZ" + ecp = ecp_for_basis(basis, ["Na", "H"]) + assert ecp == {"Na": "LANL2DZ"} + + mol = gto.M(atom=atom_str, basis=basis, ecp=ecp, charge=0, spin=0, verbose=0) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + + def _run(ecp_arg): + tmp = tempfile.NamedTemporaryFile(delete=False, suffix=".pkl") + try: + pickle.dump(None, tmp) + tmp.close() + init_raman_worker( + atom_str, + basis, + 0, + 0, + None, + tmp.name, + 1, + False, + False, + None, + ecp_arg, + ) + return run_displaced_polarizability("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + alpha_with_ecp = _run(ecp) + alpha_without_ecp = _run({}) + assert abs(alpha_with_ecp[2][2] - alpha_without_ecp[2][2]) > 1.0 + + +class TestRamanWorkerScfRescue: + def test_scf_rescue_option_honored(self, monkeypatch): + """AUDIT F19 regression — the Raman worker always called + ``run_scf_with_rescue(mf, dm0=dm0)``, taking its default + ``rescue=True`` regardless of what the caller (raman_calc.py's + ``scf_rescue`` parameter) requested. Confirms the worker now + threads that choice through. + """ + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto + + import quantui.scf_robust as scf_robust + from quantui.freq_raman_workers import ( + init_raman_worker, + run_displaced_polarizability, + ) + + 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) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + + seen_rescue = [] + _orig = scf_robust.run_scf_with_rescue + + def _spy(mf_arg, **kwargs): + seen_rescue.append(kwargs.get("rescue", True)) + return _orig(mf_arg, **kwargs) + + monkeypatch.setattr(scf_robust, "run_scf_with_rescue", _spy) + + tmp = tempfile.NamedTemporaryFile(delete=False, suffix=".pkl") + try: + pickle.dump(None, tmp) + tmp.close() + init_raman_worker( + atom_str, + "sto-3g", + 0, + 0, + None, + tmp.name, + 1, + False, # dm0_is_unrestricted + False, # density_fit_used + None, # checkpoint_items_dir + None, # ecp + False, # scf_rescue + ) + run_displaced_polarizability("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + assert seen_rescue == [False] + + def test_scf_rescue_defaults_true_for_an_older_caller(self, monkeypatch): + """Backward compatibility: a caller predating this fix (only + positional args through ``ecp``) must still get rescue=True, + matching the old always-on behavior exactly.""" + pytest.importorskip("pyscf") + import os + import pickle + import tempfile + + from pyscf import gto + + import quantui.scf_robust as scf_robust + from quantui.freq_raman_workers import ( + init_raman_worker, + run_displaced_polarizability, + ) + + 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) + coords = mol.atom_coords(unit="Bohr").flatten().tolist() + + seen_rescue = [] + _orig = scf_robust.run_scf_with_rescue + + def _spy(mf_arg, **kwargs): + seen_rescue.append(kwargs.get("rescue", True)) + return _orig(mf_arg, **kwargs) + + monkeypatch.setattr(scf_robust, "run_scf_with_rescue", _spy) + + tmp = tempfile.NamedTemporaryFile(delete=False, suffix=".pkl") + try: + pickle.dump(None, tmp) + tmp.close() + init_raman_worker(atom_str, "sto-3g", 0, 0, None, tmp.name, 1, False, False) + run_displaced_polarizability("d000_x_+", coords) + finally: + os.unlink(tmp.name) + + assert seen_rescue == [True] diff --git a/tests/test_ir_plot.py b/tests/test_ir_plot.py index 1e5b61a..a8853b4 100644 --- a/tests/test_ir_plot.py +++ b/tests/test_ir_plot.py @@ -123,6 +123,52 @@ def test_xaxis_low_to_high_broadened(self): # --------------------------------------------------------------------------- +class TestRangeCoversRealModes: + """AUDIT additional-concerns — the fixed 400-4000 cm⁻¹ default used to + be the ONLY range: a real water/STO-3G calculation has O-H stretches + at 4486.7/4788.3 cm⁻¹, silently clipped off the right edge in stick + mode and never entering the broadened kernel at all (evaluated only + on that fixed grid). The default must widen to cover every real mode. + """ + + def test_default_range_still_400_4000_when_all_modes_fit(self): + fig = plot_ir_spectrum(_SIMPLE_FREQS, _SIMPLE_INTS) + x_range = list(fig.layout.xaxis.range) + assert x_range[0] == 400 + assert x_range[1] == 4000 + + def test_stick_range_widens_for_a_high_frequency_mode(self): + # Real RHF/STO-3G water O-H stretch region. + freqs = [1785.6, 4486.7, 4788.3] + ints = [65.0, 5.0, 60.0] + fig = plot_ir_spectrum(freqs, ints, mode="stick") + x_range = list(fig.layout.xaxis.range) + assert x_range[1] > 4788.3, "x-axis must extend past the highest real mode" + x_data = [x for x in fig.data[0].x if x is not None] + assert 4788.3 in x_data, "the high-frequency stick must actually be plotted" + + def test_broadened_grid_widens_and_the_high_mode_is_broadened(self): + freqs = [1785.6, 4486.7, 4788.3] + ints = [65.0, 5.0, 60.0] + fig = plot_ir_spectrum(freqs, ints, mode="broadened", fwhm=20.0) + x = np.array(fig.data[0].x) + y = np.array(fig.data[0].y) + assert x.max() > 4788.3 + # The broadened trace must have real signal near 4788.3, not just + # a flat zero baseline past the old fixed grid's 4000 cm⁻¹ edge. + near_peak = y[(x > 4780) & (x < 4800)] + assert near_peak.max() > 1.0 + + def test_low_frequency_mode_below_400_is_not_clipped(self): + freqs = [150.0, 1500.0] + ints = [20.0, 40.0] + fig = plot_ir_spectrum(freqs, ints, mode="stick") + x_range = list(fig.layout.xaxis.range) + assert x_range[0] < 150.0 + x_data = [x for x in fig.data[0].x if x is not None] + assert 150.0 in x_data + + class TestEmptyInput: def test_empty_frequencies_no_exception(self): fig = plot_ir_spectrum([], []) diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py index 04e2e2f..161aa0f 100644 --- a/tests/test_optimizer.py +++ b/tests/test_optimizer.py @@ -34,16 +34,17 @@ except ImportError: pass -_ASE_PYSCF_AVAILABLE = False -try: - from ase.calculators.pyscf import PySCF as _c # noqa: F401 - - _ASE_PYSCF_AVAILABLE = True -except ImportError: - pass - +# AUDIT F22 — this skip used to also require `ase.calculators.pyscf.PySCF` +# (a module ASE itself does not ship), a check left over from an earlier +# design. quantui.optimizer implements and uses its own ASE Calculator, +# `_QuantUIPySCFCalc` (see optimizer.py's module docstring), and never +# imports ase.calculators.pyscf at all. Gating on it meant every test below +# skipped with a misleading "not installed" reason on any environment +# where ASE and PySCF were both genuinely installed and working — +# 16 tests, per the audit — because that unrelated module doesn't exist in +# any released ASE version. pyscf_only = pytest.mark.skipif( - not (ASE_AVAILABLE and _PYSCF_AVAILABLE and _ASE_PYSCF_AVAILABLE), + not (ASE_AVAILABLE and _PYSCF_AVAILABLE), reason="ase>=3.22 and pyscf not both installed (Linux/WSL only)", ) @@ -401,6 +402,48 @@ def test_max_steps_respected(self): assert result.n_steps <= 1 +class TestOptimizeGeometryRejectsUnconvergedScf: + """AUDIT F09 — BFGS must not accept forces from an unconverged SCF.""" + + # Not @pyscf_only: that marker gates on ase.calculators.pyscf.PySCF, + # which QuantUI's own _QuantUIPySCFCalc doesn't need (AUDIT F22) — this + # test only needs ASE + PySCF themselves, which are both available here. + @pytest.mark.skipif( + not (ASE_AVAILABLE and _PYSCF_AVAILABLE), + reason="ase and pyscf not both installed (Linux/WSL only)", + ) + @pytest.mark.slow + def test_unconverged_scf_raises_instead_of_optimizing(self, monkeypatch): + """Controlled reproduction: start H2 near its optimized geometry and + limit every real SCF to one cycle with rescue disabled (matching + the audit's reproduction). Before the fix, this silently reported + converged=True after a few steps despite every SCF evaluation being + unconverged; it must now raise instead. + """ + import pyscf.scf as pyscf_scf + + from quantui.optimizer import optimize_geometry + + _original_rhf = pyscf_scf.RHF + + def _one_cycle_rhf(mol): + mf = _original_rhf(mol) + mf.max_cycle = 1 + return mf + + monkeypatch.setattr(pyscf_scf, "RHF", _one_cycle_rhf) + + with pytest.raises(RuntimeError, match="did not converge"): + optimize_geometry( + _h2(0.74), + method="RHF", + basis="STO-3G", + fmax=0.05, + steps=10, + scf_rescue=False, + ) + + class TestOptimizeGeometryMetadataPreservation: """Charge and multiplicity survive through the optimization.""" @@ -553,6 +596,32 @@ def test_optimize_geometry_importable_from_quantui(self): from quantui import optimize_geometry # noqa: F401 +class TestPyscfOnlyGateNotObsolete: + """AUDIT F22 regression — the gate used to also require + ``ase.calculators.pyscf.PySCF``, a module ASE does not ship (this repo + implements its own ``_QuantUIPySCFCalc`` and never imports it). That + made every ``@pyscf_only`` test in this file skip with a misleading + "not installed" reason on any environment where ASE and PySCF were + both genuinely installed and working — 16 tests, per the audit. + """ + + def test_gate_does_not_require_ase_calculators_pyscf(self): + import importlib.util + + assert ( + importlib.util.find_spec("ase.calculators.pyscf") is None + ), "test environment assumption changed: ase.calculators.pyscf now exists" + # The gate must still evaluate to "available" here, since real ASE + # + PySCF ARE installed in this environment (see the module-level + # ASE_AVAILABLE / _PYSCF_AVAILABLE probes above) — it must not + # depend on the module asserted absent just above. + assert ASE_AVAILABLE and _PYSCF_AVAILABLE + assert pyscf_only.args[0] is False, ( + "pyscf_only should not skip when ASE and PySCF are both " + "installed, regardless of ase.calculators.pyscf" + ) + + # ============================================================================ # Run directly # ============================================================================ diff --git a/tests/test_orbital_visualization.py b/tests/test_orbital_visualization.py index 5bf7f17..6402b90 100644 --- a/tests/test_orbital_visualization.py +++ b/tests/test_orbital_visualization.py @@ -13,6 +13,7 @@ from quantui.orbital_visualization import ( HARTREE_TO_EV, + infer_charge_and_spin, load_orbital_info, orbital_info_from_arrays, orbital_summary_html, @@ -99,10 +100,73 @@ def minimal_cube_file(tmp_path): # --------------------------------------------------------------------------- -# OrbitalInfo construction +# infer_charge_and_spin — AUDIT F14 # --------------------------------------------------------------------------- +class TestInferChargeAndSpin: + """AUDIT F14 — a 1-D occupation array is not necessarily closed-shell + (ROHF is 1-D too), and a bare atomic-number sum overcounts an ECP + system's charge by the core electrons the ECP replaced. + """ + + def test_closed_shell_rhf_gives_zero_spin(self): + # Water RHF: 5 doubly-occupied MOs, no singly-occupied ones. + occ = [2.0, 2.0, 2.0, 2.0, 2.0] + mol_atom = [("O", [0, 0, 0]), ("H", [0, 0, 1]), ("H", [0, 1, 0])] + assert infer_charge_and_spin(mol_atom, occ) == (0, 0) + + def test_rohf_open_shell_1d_occupation_gives_nonzero_spin(self): + """The audit's exact reproduction: real OH doublet occupations + [2,2,2,2,1,0] used to infer (charge, spin)=(0,0) — impossible for + 9 electrons — because a 1-D array was assumed closed-shell. Must + now infer spin=1 (doublet) from the one singly-occupied orbital. + """ + occ = [2.0, 2.0, 2.0, 2.0, 1.0, 0.0] + mol_atom = [("O", [0, 0, 0]), ("H", [0, 0, 1])] + assert infer_charge_and_spin(mol_atom, occ) == (0, 1) + + def test_uhf_2d_occupation_still_works(self): + """The 2-D (UHF/UKS) path is unaffected by this fix.""" + occ = np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 0.0]]) + mol_atom = [("O", [0, 0, 0]), ("H", [0, 0, 1])] + charge, spin = infer_charge_and_spin(mol_atom, occ) + assert spin == 1 + assert charge == 8 + 1 - 5 # nuclear charge minus 5 electrons + + def test_ecp_system_without_basis_overcounts_charge(self): + """Without a basis to resolve the ECP, charge inference falls back + to the pre-fix all-electron behavior — documented, not silently + "fixed" by a guess it has no data to make.""" + occ = [2.0] # 2 explicit (valence) electrons + mol_atom = [("Na", [0, 0, 0]), ("H", [0, 0, 2.0])] + charge, _spin = infer_charge_and_spin(mol_atom, occ) + assert charge == 10 # 11 (Na) + 1 (H) - 2 electrons, all-electron count + + def test_ecp_system_with_basis_infers_correct_neutral_charge(self): + """The audit's exact reproduction: NaH/LANL2DZ used to infer +10 + instead of neutral, because 10 Na core electrons replaced by the + ECP were never subtracted from the atomic-number sum. + """ + pytest.importorskip("pyscf") + occ = [2.0] # 2 explicit electrons (correct for NaH/LANL2DZ) + mol_atom = [("Na", [0, 0, 0]), ("H", [0, 0, 2.0])] + charge, spin = infer_charge_and_spin(mol_atom, occ, basis="LANL2DZ") + assert charge == 0 + assert spin == 0 + + def test_all_electron_basis_unaffected_by_ecp_lookup(self): + pytest.importorskip("pyscf") + occ = [2.0, 2.0, 2.0, 2.0, 2.0] + mol_atom = [("O", [0, 0, 0]), ("H", [0, 0, 1]), ("H", [0, 1, 0])] + charge, spin = infer_charge_and_spin(mol_atom, occ, basis="STO-3G") + assert (charge, spin) == (0, 0) + + def test_none_inputs_return_zero_zero(self): + assert infer_charge_and_spin(None, [2.0]) == (0, 0) + assert infer_charge_and_spin([("H", [0, 0, 0])], None) == (0, 0) + + class TestLoadOrbitalInfo: def test_basic_load(self, simple_mo_data): diff --git a/tests/test_packaging.py b/tests/test_packaging.py new file mode 100644 index 0000000..706e926 --- /dev/null +++ b/tests/test_packaging.py @@ -0,0 +1,155 @@ +"""Wheel-content smoke tests (AUDIT F20). + +Every other test in this suite runs against an **editable** install +(``pip install -e .``, per CI and CLAUDE.md's setup instructions), which +never consults ``[tool.setuptools]``'s package list at all — the source +tree is used as-is. That is exactly why the bug this file guards against +went unnoticed: ``pyproject.toml`` explicitly listed +``packages = ["quantui", "quantui.backends"]``, silently omitting the real +``quantui.engines`` subpackage. A *built* wheel installed non-editably +contained no engine files at all, while CI — always editable — stayed +green. + +These tests build a real wheel via the PEP 517 frontend and inspect its +contents directly, so a future subpackage that isn't picked up by +discovery (or a future explicit list that forgets one) fails here instead +of only in an installed deployment. +""" + +from __future__ import annotations + +import subprocess +import sys +import zipfile +from pathlib import Path + +import pytest + +REPO_ROOT = Path(__file__).parent.parent + + +def _build_wheel(tmp_path: Path) -> Path: + """Build the project wheel into ``tmp_path`` and return its path.""" + subprocess.run( + [sys.executable, "-m", "build", "--wheel", "--outdir", str(tmp_path)], + cwd=REPO_ROOT, + check=True, + capture_output=True, + text=True, + ) + wheels = list(tmp_path.glob("*.whl")) + assert len(wheels) == 1, f"expected exactly one built wheel, found {wheels}" + return wheels[0] + + +@pytest.mark.slow +class TestWheelContents: + """Inspect the built wheel's file list directly (no install needed).""" + + def test_engines_subpackage_is_shipped(self, tmp_path): + """AUDIT F20 — quantui/engines/* must be present in the wheel. + + This is the exact regression: the package list omitted + quantui.engines, so a real (non-editable) install had no engine + module at all. + """ + pytest.importorskip("build") + wheel_path = _build_wheel(tmp_path) + with zipfile.ZipFile(wheel_path) as zf: + names = zf.namelist() + engine_files = [n for n in names if n.startswith("quantui/engines/")] + assert engine_files, f"no quantui/engines/* files in wheel: {names}" + assert "quantui/engines/__init__.py" in engine_files + assert "quantui/engines/base.py" in engine_files + assert "quantui/engines/pyscf_engine.py" in engine_files + assert "quantui/engines/pyfock_engine.py" in engine_files + + def test_backends_subpackage_is_shipped(self, tmp_path): + """Same check for quantui.backends, which the old explicit list did + include — a regression guard so a future edit can't drop it either.""" + pytest.importorskip("build") + wheel_path = _build_wheel(tmp_path) + with zipfile.ZipFile(wheel_path) as zf: + names = zf.namelist() + assert "quantui/backends/__init__.py" in names + assert "quantui/backends/worker.py" in names + + def test_tests_package_is_not_shipped(self, tmp_path): + """``tests/`` has its own __init__.py (it's a real Python package), + so package discovery must explicitly exclude it — otherwise the + wheel ships the whole test suite as an importable top-level + package, which is not what a distributed application wheel wants. + """ + pytest.importorskip("build") + wheel_path = _build_wheel(tmp_path) + with zipfile.ZipFile(wheel_path) as zf: + names = zf.namelist() + test_files = [n for n in names if n.startswith("tests/")] + assert not test_files, f"tests/* leaked into the wheel: {test_files}" + + def test_bundled_data_files_are_shipped(self, tmp_path): + """Package-data (the offline molecule library + vendored 3Dmol.js) + must survive alongside the package-discovery change — this is a + pre-existing guarantee (not part of F20), reasserted here since + this file is the one place a wheel is actually built and + inspected.""" + pytest.importorskip("build") + wheel_path = _build_wheel(tmp_path) + with zipfile.ZipFile(wheel_path) as zf: + names = zf.namelist() + assert "quantui/data/library/library.sqlite" in names + assert "quantui/data/js/3Dmol-min.js" in names + + +@pytest.mark.slow +class TestWheelIsolatedImport: + """Install the built wheel into a throwaway venv and import from + outside the checkout — the audit's own reproduction: "Importing + quantui.engines from the extracted wheel, with the editable source + finder removed from the probe process, raises ModuleNotFoundError." + """ + + def test_engines_importable_from_installed_wheel(self, tmp_path): + pytest.importorskip("build") + wheel_path = _build_wheel(tmp_path / "dist") + + venv_dir = tmp_path / "venv" + subprocess.run( + [sys.executable, "-m", "venv", str(venv_dir)], + check=True, + capture_output=True, + text=True, + ) + venv_python = venv_dir / "bin" / "python" + if not venv_python.exists(): # pragma: no cover - Windows layout + venv_python = venv_dir / "Scripts" / "python.exe" + + subprocess.run( + [str(venv_python), "-m", "pip", "install", "--quiet", str(wheel_path)], + check=True, + capture_output=True, + text=True, + ) + + # Run the import check from tmp_path, NOT the repo checkout, so + # Python can't fall back to the source tree on sys.path[0] and + # mask a genuinely broken installed package (the audit's own + # methodology: "with the editable source finder removed from the + # probe process"). + result = subprocess.run( + [ + str(venv_python), + "-c", + "import quantui.engines as e; " + "from quantui.engines import base, pyscf_engine, pyfock_engine; " + "print('OK', e.__file__)", + ], + cwd=str(tmp_path), + capture_output=True, + text=True, + ) + assert result.returncode == 0, ( + f"quantui.engines not importable from the installed wheel:\n" + f"stdout={result.stdout}\nstderr={result.stderr}" + ) + assert "OK" in result.stdout diff --git a/tests/test_pes_scan.py b/tests/test_pes_scan.py index 9fc7487..375a16b 100644 --- a/tests/test_pes_scan.py +++ b/tests/test_pes_scan.py @@ -400,6 +400,108 @@ def test_format_shows_converged_yes(self): ) +@_pyscf_available +@pytest.mark.slow +class TestRunPesScanRejectsUnconvergedScf: + """AUDIT F09 — a diatomic bond scan point hard-codes ``ok = True`` + (there's no relaxable DOF, so *geometric* convergence is trivial), but + that must not paper over the *electronic* SCF failing at that point. + """ + + @pytest.mark.slow + def test_unconverged_scf_marks_scan_point_as_failed(self, monkeypatch): + """Controlled reproduction, mirroring the optimizer's own F09 test: + force every SCF to one cycle with rescue disabled. Before the + optimizer-level F09 fix, this diatomic scan point's ``ok`` was + hard-set True regardless, so the failure was invisible. + """ + pytest.importorskip("pyscf") + import pyscf.scf as pyscf_scf + + from quantui.pes_scan import run_pes_scan + + _original_rhf = pyscf_scf.RHF + + def _one_cycle_rhf(mol): + mf = _original_rhf(mol) + mf.max_cycle = 1 + return mf + + monkeypatch.setattr(pyscf_scf, "RHF", _one_cycle_rhf) + + result = run_pes_scan( + _h2(), + method="RHF", + basis="STO-3G", + scan_type="bond", + atom_indices=[0, 1], + start=0.6, + stop=1.4, + steps=3, + scf_rescue=False, + ) + + assert result.converged_all is False + import math + + assert all(math.isnan(e) for e in result.energies_hartree) + + +class TestRunPesScanCheckpointCoordinateIdentity: + """AUDIT F10 — a cached checkpoint point must not be reused for a + different scan configuration just because its scalar "value" happens + to match the current target. + """ + + @pytest.mark.slow + def test_mismatched_scan_type_point_is_not_reused(self): + """Reproduces the audit's scenario: a banked point recorded under + a different scan_type at the same target value must be recomputed + for real, not silently returned as-is. + """ + pytest.importorskip("pyscf") + from quantui.checkpoint import CalcIdentity, Checkpoint + from quantui.pes_scan import run_pes_scan + + identity = CalcIdentity.from_molecule( + _h2(), calc_type="pes_scan", method="RHF", basis="STO-3G" + ) + ckpt = Checkpoint(identity) + ckpt.begin() + # A bogus point banked under a different scan_type, at the exact + # value the H-H bond scan's first grid point will target, with an + # energy nowhere near a real RHF/STO-3G H2 energy — if this were + # wrongly reused, the returned energy would be this exact value. + ckpt.append_point( + { + "index": 1, + "value": 0.6, + "scan_type": "angle", + "atom_indices": [0, 1, 0], + "energy_hartree": -999.0, + "ok": True, + "atoms": ["H", "H"], + "coordinates": [[0.0, 0.0, 0.0], [0.0, 0.0, 0.6]], + } + ) + + result = run_pes_scan( + _h2(), + method="RHF", + basis="STO-3G", + scan_type="bond", + atom_indices=[0, 1], + start=0.6, + stop=1.4, + steps=3, + checkpoint=ckpt, + resume=True, + ) + + assert result.energies_hartree[0] != pytest.approx(-999.0) + assert result.energies_hartree[0] < -0.5 # a real RHF/STO-3G H2 energy + + @_pyscf_available @pytest.mark.slow class TestRunPesScanIntegration: diff --git a/tests/test_raman_calc.py b/tests/test_raman_calc.py index 6b063d1..fc6cba1 100644 --- a/tests/test_raman_calc.py +++ b/tests/test_raman_calc.py @@ -34,6 +34,76 @@ def test_non_negative(self): assert _raman_invariants(da) >= 0.0 +class TestRamanUnitsRegression: + """AUDIT F02 — CPU Raman activities were ~45.54x too large because the + polarizability numerator (a0^3) was never converted to Angstrom^3, only + the displacement denominator (Bohr -> Angstrom) was. This reproduces the + audit's own method: an independent central difference of the real SCF + polarizability along the complete (mass-normalized) normal mode, using + none of raman_calc's atom-by-atom Jacobian/unit-conversion code. + """ + + @pytest.mark.slow + def test_h2_raman_activity_matches_independent_normal_mode_fd(self): + pytest.importorskip("pyscf") + pytest.importorskip("pyscf.prop.polarizability.rhf") + import numpy as _np + from pyscf import gto, scf + from pyscf.prop.polarizability import rhf as pol_mod + + from quantui.config import BOHR_TO_ANGSTROM as _BOHR_TO_ANG + from quantui.freq_calc import run_freq_calc + from quantui.molecule import Molecule + + bond_length = 0.74 # Angstrom + h2 = Molecule(["H", "H"], [[0.0, 0.0, 0.0], [0.0, 0.0, bond_length]]) + result = run_freq_calc(h2, method="RHF", basis="STO-3G") + assert result.raman_activities, "H2/RHF/STO-3G should produce a Raman activity" + assert result.displacements, "normal-mode displacements required for FD check" + + nm_flat = _np.asarray(result.displacements[0], dtype=float).reshape(-1) + + def _alpha_au(coords_bohr: _np.ndarray) -> _np.ndarray: + mol = gto.M( + atom=[ + ("H", tuple(coords_bohr[0])), + ("H", tuple(coords_bohr[1])), + ], + basis="STO-3G", + unit="Bohr", + verbose=0, + ) + mf = scf.RHF(mol) + mf.verbose = 0 + mf.kernel() + return _np.asarray( + pol_mod.polarizability(pol_mod.Polarizability(mf)), dtype=float + ) + + coords0_bohr = _np.array( + [[0.0, 0.0, 0.0], [0.0, 0.0, bond_length / _BOHR_TO_ANG]] + ) + # Scale the step so the largest per-atom Cartesian displacement is a + # modest 0.01 Bohr, regardless of how PySCF normalizes the mode. + eps = 0.01 / max(abs(nm_flat.max()), abs(nm_flat.min()), 1e-12) + disp_bohr = (eps * nm_flat).reshape(-1, 3) + + alpha_plus = _alpha_au(coords0_bohr + disp_bohr) + alpha_minus = _alpha_au(coords0_bohr - disp_bohr) + # d(alpha[a0^3]) / d(eps) == sum_k nm_k * d(alpha[a0^3])/d(x_k[Bohr]) + # by the chain rule — exactly the per-atom Jacobian raman_calc.py + # projects onto the same nm vector, just computed by directly + # perturbing along the mode instead of atom-by-atom. + dalpha_dq_au = (alpha_plus - alpha_minus) / (2.0 * eps) + dalpha_dq_ang = dalpha_dq_au * (_BOHR_TO_ANG**2) + + expected_activity = _raman_invariants(dalpha_dq_ang) + assert result.raman_activities[0] == pytest.approx(expected_activity, rel=0.05) + # The pre-fix bug deflated this by ~45.54x — well outside any + # plausible finite-difference discrepancy. + assert result.raman_activities[0] > 0.2 * expected_activity + + class TestRamanEnabled: def test_default_enabled(self, monkeypatch): monkeypatch.delenv("QUANTUI_RAMAN", raising=False) diff --git a/tests/test_raman_plot.py b/tests/test_raman_plot.py index 3a556bd..a8b7657 100644 --- a/tests/test_raman_plot.py +++ b/tests/test_raman_plot.py @@ -26,3 +26,25 @@ def test_skips_imaginary_frequencies(): fig = plot_raman_spectrum([-100.0, 1500.0], [5.0, 20.0], mode="stick") assert len(fig.data) == 2 assert fig.data[1].x == (1500.0,) + + +def test_range_widens_for_a_high_frequency_mode(): + """AUDIT additional-concerns — mirrors test_ir_plot.py: a real O-H + stretch above 4000 cm⁻¹ must not be clipped off the fixed default + range.""" + freqs = [1785.6, 4486.7, 4788.3] + acts = [2.0, 30.0, 15.0] + fig = plot_raman_spectrum(freqs, acts, mode="stick") + x_range = list(fig.layout.xaxis.range) + assert x_range[1] > 4788.3 + assert 4788.3 in fig.data[1].x + + +def test_broadened_grid_widens_for_a_high_frequency_mode(): + import numpy as np + + freqs = [1785.6, 4486.7, 4788.3] + acts = [2.0, 30.0, 15.0] + fig = plot_raman_spectrum(freqs, acts, mode="broadened", fwhm=20.0) + x = np.array(fig.data[0].x) + assert x.max() > 4788.3 diff --git a/tests/test_session_calc.py b/tests/test_session_calc.py index c69176d..6db3c74 100644 --- a/tests/test_session_calc.py +++ b/tests/test_session_calc.py @@ -595,6 +595,66 @@ def test_ccsd_t_water_runs_and_reports_triples(self): assert result.ccsd_t_correction_hartree is not None assert result.ccsd_t_correction_hartree < 0 + @pyscf_only + def test_ccsd_nonconvergence_is_not_reported_as_converged(self, monkeypatch): + """AUDIT F07 regression — controlled reproduction: restrict the + real CCSD solver to one iteration (a real, unpatched PySCF kernel; + only the iteration limit is forced). Before the fix, the returned + correlation energy was accepted unconditionally and + result.converged came solely from the HF reference, so this + reported converged=True despite _ccsd_obj.converged being False. + """ + import pyscf.cc as pyscf_cc + + from quantui.session_calc import run_in_session + + _original_ccsd = pyscf_cc.CCSD + + def _one_cycle_ccsd(mf): + obj = _original_ccsd(mf) + obj.max_cycle = 1 + return obj + + monkeypatch.setattr(pyscf_cc, "CCSD", _one_cycle_ccsd) + + result = run_in_session(molecule=_water(), method="CCSD", basis="STO-3G") + + assert result.cc_converged is False + assert result.converged is False + # The correlation energy is still surfaced (so the UI can show what + # happened) — just not stamped as a converged result. + assert result.ccsd_correlation_hartree is not None + + @pyscf_only + def test_ccsd_skipped_when_reference_scf_unconverged(self, monkeypatch): + """AUDIT F07 — post-HF work must not launch on an unconverged SCF. + + Runs the real SCF (so ``mf`` has genuine, valid orbitals), then + flips ``mf.converged`` to simulate a reference reported as + unconverged — session_calc.py reads exactly that attribute to + decide whether to launch CCSD, so this exercises the actual gate + without needing to engineer real SCF non-convergence. + """ + import quantui.scf_robust as scf_robust + from quantui.session_calc import run_in_session + + _real_run_scf_with_rescue = scf_robust.run_scf_with_rescue + + def _fake_run_scf_with_rescue(mf, *args, **kwargs): + energy = _real_run_scf_with_rescue(mf, *args, **kwargs) + mf.converged = False + return energy + + monkeypatch.setattr( + scf_robust, "run_scf_with_rescue", _fake_run_scf_with_rescue + ) + + result = run_in_session(molecule=_water(), method="CCSD", basis="STO-3G") + + assert result.converged is False + assert result.cc_converged is None + assert result.ccsd_correlation_hartree is None + # ============================================================================ # Run directly diff --git a/tests/test_slurm_ingest.py b/tests/test_slurm_ingest.py index 850b235..774ac63 100644 --- a/tests/test_slurm_ingest.py +++ b/tests/test_slurm_ingest.py @@ -30,6 +30,56 @@ def test_ingest_single_point(self, tmp_path, monkeypatch): assert data["calc_type"] == "single_point" assert (saved / "pyscf.log").read_text() == "log line\n" + def test_ingest_single_point_preserves_enriched_fields(self, tmp_path, monkeypatch): + """AUDIT F12 regression — a real water single-point round trip + (matching the audit's own reproduction) must not lose Mulliken + charges, dipole, atom symbols, SCF provenance, solvent/GPU/DF + metadata, or post-HF correlation fields between the worker's + staging JSON and the saved History result.json. + """ + patch_results_root(tmp_path, monkeypatch) + payload = { + "calc_type": "single_point", + "energy_hartree": -76.023190, + "homo_lumo_gap_ev": 10.0, + "converged": True, + "n_iterations": 6, + "method": "RHF", + "basis": "STO-3G", + "formula": "H2O", + "mulliken_charges": [-0.365510, 0.182755, 0.182755], + "dipole_moment_debye": 1.725515, + "dipole_vector_debye": [0.0, 1.725515, 0.0], + "atom_symbols": ["O", "H", "H"], + "scf_rescue_stage": "bootstrap", + "scf_variant": "RHF", + "mp2_correlation_hartree": -0.201, + "ccsd_correlation_hartree": -0.213, + "ccsd_t_correction_hartree": -0.004, + "cc_converged": True, + "dispersion_applied": False, + "solvent": "Water", + "gpu_used": True, + "gpu_name": "NVIDIA H200", + "density_fit": True, + } + record, _staging = make_staging_record(tmp_path, payload) + saved = ingest_staging_success(record) + data = json.loads((saved / "result.json").read_text()) + + assert data["mulliken_charges"] == [-0.365510, 0.182755, 0.182755] + assert data["dipole_moment_debye"] == 1.725515 + assert data["dipole_vector_debye"] == [0.0, 1.725515, 0.0] + assert data["atom_symbols"] == ["O", "H", "H"] + assert data["scf_variant"] == "RHF" + assert data["mp2_correlation_hartree"] == -0.201 + assert data["ccsd_correlation_hartree"] == -0.213 + assert data["ccsd_t_correction_hartree"] == -0.004 + assert data["solvent"] == "Water" + assert data["gpu_used"] is True + assert data["gpu_name"] == "NVIDIA H200" + assert data["density_fit"] is True + def test_ingest_geometry_opt_copies_trajectory(self, tmp_path, monkeypatch): patch_results_root(tmp_path, monkeypatch) traj = { diff --git a/tests/test_tddft_calc.py b/tests/test_tddft_calc.py index 63f754f..e7fb913 100644 --- a/tests/test_tddft_calc.py +++ b/tests/test_tddft_calc.py @@ -79,5 +79,64 @@ def test_scf_variant_reports_rks_for_closed_shell_dft(self): assert result.scf_variant == "RKS" +# ============================================================================ +# AUDIT F08 — per-root TD convergence must not be swallowed +# ============================================================================ + + +class TestTddftConvergence: + @pyscf_only + @pytest.mark.slow + def test_unconverged_roots_are_not_reported_as_converged(self, monkeypatch): + """Controlled reproduction from the audit: restrict the real TDHF + Davidson solve to one iteration (only max_cycle forced; the SCF and + TD kernels themselves are real PySCF). A real one-iteration-limited + TDHF/6-31G water solve gives converged=[False, False, False] and + excitations 9.804744, 11.993916, 12.420652 eV — the old code + reported these as a converged result because it never checked + td.converged at all. + """ + import pyscf.scf.hf as pyscf_hf + import pyscf.tdscf.rhf # noqa: F401 — import side effect registers RHF.TDHF + + from quantui.tddft_calc import run_tddft_calc + + # mf.TDHF() is registered via pyscf.lib.class_as_method, which binds + # the TDHF class into RHF.TDHF as a plain function at pyscf import + # time — monkeypatching pyscf.tdscf.rhf.TDHF afterward has no effect + # on that already-captured reference, so patch RHF.TDHF itself. + _original_tdhf_method = pyscf_hf.RHF.TDHF + + def _one_cycle_tdhf(self): + obj = _original_tdhf_method(self) + obj.max_cycle = 1 + return obj + + monkeypatch.setattr(pyscf_hf.RHF, "TDHF", _one_cycle_tdhf) + + result = run_tddft_calc(_water(), method="RHF", basis="6-31G", nstates=3) + + assert result.td_converged == [False, False, False] + assert result.n_converged_states == 0 + assert result.converged is False + # The excitations are still surfaced (so the UI can show what + # actually came out of the solver) — just not stamped converged. + assert len(result.excitation_energies_ev) == 3 + + @pyscf_only + @pytest.mark.slow + def test_converged_roots_report_full_convergence(self): + """Sanity check the happy path: a normal (unpatched) TDHF solve on + a small system converges every requested root.""" + from quantui.tddft_calc import run_tddft_calc + + result = run_tddft_calc(_water(), method="RHF", basis="STO-3G", nstates=2) + + assert result.td_converged is not None + assert all(result.td_converged) + assert result.n_converged_states == len(result.td_converged) + assert result.converged is True + + if __name__ == "__main__": pytest.main([__file__, "-v", "--tb=short"]) diff --git a/tests/test_worker_payload.py b/tests/test_worker_payload.py index 042b171..8d83c7b 100644 --- a/tests/test_worker_payload.py +++ b/tests/test_worker_payload.py @@ -19,6 +19,7 @@ session_result_payload, tddft_result_payload, ) +from quantui.freq_calc import ThermoData from quantui.molecule import Molecule @@ -37,6 +38,15 @@ def _session_result(**overrides) -> SimpleNamespace: atom_symbols=["Mn", "O", "O", "O", "O", "O", "O"], scf_rescue_stage="bootstrap", scf_variant="UKS", + mp2_correlation_hartree=-0.201, + ccsd_correlation_hartree=-0.213, + ccsd_t_correction_hartree=-0.004, + cc_converged=True, + dispersion_applied=False, + solvent="Water", + gpu_used=True, + gpu_name="NVIDIA H200", + density_fit=True, ) defaults.update(overrides) return SimpleNamespace(**defaults) @@ -102,6 +112,16 @@ def test_missing_from_an_older_sessionresult_defaults_to_none(self): assert payload["atom_symbols"] is None assert payload["scf_rescue_stage"] == "none" assert payload["scf_variant"] is None + # AUDIT F12 + assert payload["mp2_correlation_hartree"] is None + assert payload["ccsd_correlation_hartree"] is None + assert payload["ccsd_t_correction_hartree"] is None + assert payload["cc_converged"] is None + assert payload["dispersion_applied"] is None + assert payload["solvent"] is None + assert payload["gpu_used"] is False + assert payload["gpu_name"] is None + assert payload["density_fit"] is False def test_scf_rescue_stage_present(self): payload = session_result_payload(_session_result()) @@ -118,6 +138,21 @@ def test_calc_type_and_core_fields_unchanged(self): assert payload["converged"] is True assert payload["formula"] == "Mn(H2O)6" + def test_post_hf_and_solvent_gpu_df_fields_present(self): + """AUDIT F12 — these were computed onto SessionResult but never + serialized into staging JSON at all, distinct from (and upstream + of) _basic_result's own reconstruction gap in slurm_ingest.py.""" + payload = session_result_payload(_session_result()) + assert payload["mp2_correlation_hartree"] == -0.201 + assert payload["ccsd_correlation_hartree"] == -0.213 + assert payload["ccsd_t_correction_hartree"] == -0.004 + assert payload["cc_converged"] is True + assert payload["dispersion_applied"] is False + assert payload["solvent"] == "Water" + assert payload["gpu_used"] is True + assert payload["gpu_name"] == "NVIDIA H200" + assert payload["density_fit"] is True + class TestFreqTddftNmrResultPayloadScfVariant: """M-UX2 UXP2.10 — the same provenance field promoted to the other @@ -139,11 +174,99 @@ def test_freq_result_payload_carries_scf_variant(self): zpve_hartree=0.0, thermo=None, scf_variant="UKS", + density_fit=True, ) molecule = Molecule(["O", "H", "H"], [[0, 0, 0], [0.96, 0, 0], [0, 0.96, 0]]) payload = freq_result_payload(result, molecule) assert payload["scf_variant"] == "UKS" + def test_freq_result_payload_carries_density_fit(self): + """AUDIT F12 — density_fit was never serialized here at all, + though FreqResult carries it.""" + result = SimpleNamespace( + energy_hartree=-1600.0, + homo_lumo_gap_ev=None, + converged=True, + n_iterations=30, + method="B3LYP", + basis="def2-SVP", + formula="Fe(H2O)6", + displacements=None, + frequencies_cm1=[], + ir_intensities=[], + raman_activities=[], + zpve_hartree=0.0, + thermo=None, + scf_variant="UKS", + density_fit=True, + ) + molecule = Molecule(["O", "H", "H"], [[0, 0, 0], [0.96, 0, 0], [0, 0.96, 0]]) + payload = freq_result_payload(result, molecule) + assert payload["density_fit"] is True + + def test_freq_result_payload_carries_thermo(self): + """AUDIT F18 — FreqResult.thermo (H, S, G, ZPVE, temperature) was + computed by freq_calc.py but never made it into the staging JSON; + the batch save had only frequencies/intensities/activities/ + displacements/ZPVE, with thermochemistry silently discarded.""" + result = SimpleNamespace( + energy_hartree=-76.0, + homo_lumo_gap_ev=None, + converged=True, + n_iterations=12, + method="RHF", + basis="STO-3G", + formula="H2O", + displacements=None, + frequencies_cm1=[1600.0, 3700.0, 3800.0], + ir_intensities=[10.0, 5.0, 5.0], + raman_activities=[], + zpve_hartree=0.021, + thermo=ThermoData( + zpve_hartree=0.021, + H_hartree=-74.933498241, + S_jmol=188.538424, + G_hartree=-74.954908540, + temperature_k=298.15, + ), + scf_variant="RHF", + density_fit=False, + ) + molecule = Molecule(["O", "H", "H"], [[0, 0, 0], [0.96, 0, 0], [0, 0.96, 0]]) + payload = freq_result_payload(result, molecule) + thermo = payload["spectra"]["ir"]["thermo"] + assert thermo is not None + assert thermo["H_hartree"] == -74.933498241 + assert thermo["S_jmol"] == 188.538424 + assert thermo["G_hartree"] == -74.954908540 + assert thermo["temperature_k"] == 298.15 + assert thermo["pressure_atm"] == 1.0 + assert thermo["approximation"] == "ideal_gas_rigid_rotor_harmonic_oscillator" + + def test_freq_result_payload_thermo_none_when_missing(self): + """A Hessian-only run (or an older FreqResult) with no thermo object + must serialize a clean None, not raise.""" + result = SimpleNamespace( + energy_hartree=-76.0, + homo_lumo_gap_ev=None, + converged=True, + n_iterations=12, + method="RHF", + basis="STO-3G", + formula="H2O", + displacements=None, + frequencies_cm1=[], + ir_intensities=[], + raman_activities=[], + zpve_hartree=0.0, + thermo=None, + scf_variant="RHF", + density_fit=False, + ) + molecule = Molecule(["O", "H", "H"], [[0, 0, 0], [0.96, 0, 0], [0, 0.96, 0]]) + payload = freq_result_payload(result, molecule) + assert payload["spectra"]["ir"]["thermo"] is None + def test_tddft_result_payload_carries_scf_variant(self): result = SimpleNamespace( energy_hartree=-1600.0, @@ -161,6 +284,25 @@ def test_tddft_result_payload_carries_scf_variant(self): # tddft_calc.TDDFTResult sets scf_variant; a bare SimpleNamespace # without it must still serialize (None), not raise. assert payload["scf_variant"] is None + # AUDIT F12 — density_fit was never serialized here at all. + assert payload["density_fit"] is False + + def test_tddft_result_payload_carries_density_fit(self): + result = SimpleNamespace( + energy_hartree=-1600.0, + homo_lumo_gap_ev=None, + converged=True, + n_iterations=20, + method="B3LYP", + basis="def2-SVP", + formula="Co(H2O)6", + excitation_energies_ev=[], + oscillator_strengths=[], + wavelengths_nm=lambda: [], + density_fit=True, + ) + payload = tddft_result_payload(result) + assert payload["density_fit"] is True def test_nmr_result_payload_carries_scf_variant(self): result = SimpleNamespace( @@ -175,6 +317,9 @@ def test_nmr_result_payload_carries_scf_variant(self): reference_key="B3LYP/6-31G*", is_fallback_reference=False, scf_variant="RKS", + density_fit=True, ) payload = nmr_result_payload(result) assert payload["scf_variant"] == "RKS" + # AUDIT F12 — density_fit was never serialized here at all. + assert payload["density_fit"] is True diff --git a/tests/test_xc_resolution.py b/tests/test_xc_resolution.py index fe13fee..b56e961 100644 --- a/tests/test_xc_resolution.py +++ b/tests/test_xc_resolution.py @@ -2,11 +2,21 @@ The user's tier-3 calibration output showed ``H₂O wB97X-D/6-31G*`` erroring at 0.01 s — PySCF rejects ``mf.xc = "wb97x-d"`` because that composite -name is on the dftd3 black-list (pyscf/pyscf#2069). The fix: - -- Alias ``wB97X-D`` to bare ``wb97x``. -- Add ``wB97X-D`` to ``_NEEDS_D3`` so dispersion is applied via - ``pyscf.dftd3``, matching the UI label that already promises D3. +name is on the dftd3 black-list (pyscf/pyscf#2069). Session 55's original +fix aliased ``wB97X-D`` to bare ``wb97x`` and applied external Grimme D3 — +but that silently calculates a *different* functional (bare wb97x has +range-separation omega=0.3; the real wB97X-D has omega=0.2 and different +short-range exact exchange — AUDIT F03). The corrected fix: + +- Alias ``wB97X-D`` to its full LibXC name ``hyb_gga_xc_wb97x_d`` — the + actual Chai & Head-Gordon (2008) functional, whose own empirical + dispersion is baked into the fit. PySCF's short-alias black-list + (``pyscf.scf.dispersion.parse_dft``) intercepts "wb97x-d"/"wb97x_d" but + not the full LibXC name, so this avoids the original error without + substituting a different functional. +- ``wB97X-D`` does NOT go in ``_NEEDS_D3`` — wrapping it in + ``pyscf.dftd3`` would double-count dispersion under a method that + already includes its own. - Extract ``resolve_xc()`` + ``maybe_apply_d3()`` so every DFT entry point (session_calc / freq_calc / tddft_calc / optimizer / nmr_calc / the script-export template) shares the same resolution logic. Before @@ -21,6 +31,8 @@ import inspect +import pytest + from quantui.session_calc import ( _NEEDS_D3, _XC_ALIAS, @@ -35,15 +47,34 @@ class TestResolveXc: - def test_wb97x_d_resolves_to_bare_wb97x(self): - # The session-55 bug: PySCF rejects "wb97x-d". Bare wb97x is - # the right xc string; D3 dispersion is applied separately. - assert resolve_xc("wB97X-D") == "wb97x" + def test_wb97x_d_resolves_to_true_functional(self): + # AUDIT F03: PySCF rejects "wb97x-d" (short-alias black-list), but + # the fix must not substitute a different functional (bare wb97x) + # to work around that — it must resolve to the actual wB97X-D + # (Chai & Head-Gordon 2008) functional under its full LibXC name. + assert resolve_xc("wB97X-D") == "hyb_gga_xc_wb97x_d" + assert resolve_xc("wB97X-D") != "wb97x" def test_wb97x_d_case_insensitive(self): # Users sometimes type "WB97X-D" or "wb97x-d" — all should resolve. for spelling in ("wB97X-D", "WB97X-D", "wb97x-d", "Wb97x-D"): - assert resolve_xc(spelling) == "wb97x" + assert resolve_xc(spelling) == "hyb_gga_xc_wb97x_d" + + def test_wb97x_d_is_a_distinct_functional_from_bare_wb97x(self): + """AUDIT F03 numeric regression: resolve_xc("wB97X-D") must resolve + to a functional with different range-separation parameters than + bare wb97x, confirmed against PySCF/LibXC directly (independent of + the alias table's own claims).""" + pytest.importorskip("pyscf.dft") + from pyscf.dft import libxc + + resolved = resolve_xc("wB97X-D") + assert resolved != "wb97x" + wb97xd_omega = libxc.rsh_coeff(resolved)[0] + wb97x_omega = libxc.rsh_coeff("wb97x")[0] + assert wb97xd_omega == pytest.approx(0.2, abs=1e-6) + assert wb97x_omega == pytest.approx(0.3, abs=1e-6) + assert wb97xd_omega != wb97x_omega def test_pbe_d3_resolves_to_bare_pbe(self): # PBE-D3 is the long-standing pattern this fix mirrors. @@ -72,19 +103,21 @@ def test_unknown_method_passes_through(self): class TestNeedsD3: - def test_wb97x_d_needs_d3(self): - # The session-55 fix: wB97X-D now needs external D3. - assert needs_d3("wB97X-D") is True + def test_wb97x_d_does_not_need_external_d3(self): + # AUDIT F03: wB97X-D's dispersion is baked into the XC functional + # itself (hyb_gga_xc_wb97x_d) — wrapping it in pyscf.dftd3 would + # double-count dispersion, so it must NOT be in _NEEDS_D3. + assert needs_d3("wB97X-D") is False def test_pbe_d3_needs_d3(self): assert needs_d3("PBE-D3") is True def test_case_insensitive(self): - assert needs_d3("WB97X-D") is True + assert needs_d3("WB97X-D") is False assert needs_d3("pbe-d3") is True def test_dispersion_free_methods_dont_need_d3(self): - for method in ("RHF", "UHF", "B3LYP", "PBE0", "M06-2X", "HSE06"): + for method in ("RHF", "UHF", "B3LYP", "PBE0", "M06-2X", "HSE06", "wB97X-D"): assert needs_d3(method) is False def test_unknown_method_doesnt_need_d3(self): @@ -105,15 +138,32 @@ def __init__(self, label): class TestMaybeApplyD3: - def test_no_d3_method_returns_mf_unchanged(self): + """AUDIT F04 — maybe_apply_d3 now returns (mf, dispersion_applied) so + callers can record whether a D3-requiring result is actually missing + its dispersion correction, instead of silently keeping the original + method label on an uncorrected result.""" + + def test_no_d3_method_returns_mf_unchanged_and_none_flag(self): mf = _FakeMf("B3LYP") - result = maybe_apply_d3(mf, "B3LYP") - assert result is mf + result_mf, dispersion_applied = maybe_apply_d3(mf, "B3LYP") + assert result_mf is mf + assert dispersion_applied is None + + def test_wb97x_d_returns_mf_unchanged_and_none_flag(self): + # AUDIT F03: wB97X-D's dispersion is already in the XC functional — + # maybe_apply_d3 must be a no-op for it (never imports pyscf.dftd3), + # and dispersion_applied is None (not applicable), not False. + mf = _FakeMf("wB97X-D") + result_mf, dispersion_applied = maybe_apply_d3(mf, "wB97X-D") + assert result_mf is mf + assert dispersion_applied is None - def test_d3_method_with_missing_pyscf_returns_mf_unchanged(self, monkeypatch): + def test_d3_method_with_missing_pyscf_returns_mf_unchanged_and_false_flag( + self, monkeypatch + ): # Simulate pyscf.dftd3 being absent (typical on Windows where # PySCF isn't installable at all). The helper must return the - # original mf without raising. + # original mf, flagged dispersion_applied=False, without raising. import builtins original_import = builtins.__import__ @@ -125,10 +175,11 @@ def _fake_import(name, *args, **kwargs): monkeypatch.setattr(builtins, "__import__", _fake_import) - mf = _FakeMf("wB97X-D") + mf = _FakeMf("PBE-D3") # Without progress_stream — must not raise. - result = maybe_apply_d3(mf, "wB97X-D") - assert result is mf + result_mf, dispersion_applied = maybe_apply_d3(mf, "PBE-D3") + assert result_mf is mf + assert dispersion_applied is False def test_d3_warning_written_to_progress_stream(self, monkeypatch): import builtins @@ -144,11 +195,34 @@ def _fake_import(name, *args, **kwargs): monkeypatch.setattr(builtins, "__import__", _fake_import) stream = io.StringIO() - maybe_apply_d3(_FakeMf("wB97X-D"), "wB97X-D", progress_stream=stream) + maybe_apply_d3(_FakeMf("PBE-D3"), "PBE-D3", progress_stream=stream) out = stream.getvalue() # User must see the missing-dispersion warning. assert "dftd3 not available" in out - assert "wB97X-D" in out + assert "PBE-D3" in out + + def test_d3_warning_logged_even_without_progress_stream(self, monkeypatch, caplog): + # AUDIT F04: the optimizer path used to call maybe_apply_d3 with no + # progress_stream at all, so a missing pyscf.dftd3 gave NO warning + # anywhere. It must now always be logged, stream or not. + import builtins + import logging + + original_import = builtins.__import__ + + def _fake_import(name, *args, **kwargs): + if name == "pyscf.dftd3" or name.startswith("pyscf.dftd3"): + raise ImportError("simulated") + return original_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", _fake_import) + + with caplog.at_level(logging.WARNING, logger="quantui.session_calc"): + maybe_apply_d3(_FakeMf("PBE-D3"), "PBE-D3") + + assert any( + "dftd3 not available" in rec.message for rec in caplog.records + ), caplog.text # ===================================================================== @@ -170,7 +244,8 @@ def test_session_calc_uses_resolve_xc(self): src = inspect.getsource(session_calc) assert "resolve_xc(method)" in src - assert "maybe_apply_d3(mf, method" in src + assert "maybe_apply_d3(" in src + assert "mf, method, progress_stream=progress_stream" in src def test_freq_calc_uses_resolve_xc(self): from quantui import freq_calc @@ -210,13 +285,20 @@ def test_script_template_embeds_alias_resolution(self): # inlined. from quantui.config import PYSCF_SCRIPT_TEMPLATE - # The literal alias for wB97X-D in the template should be the - # bare functional (post-session-55 fix). Doubled-brace literals - # in the template appear as single braces in the output. - assert "'wB97X-D': 'wb97x'" in PYSCF_SCRIPT_TEMPLATE + # The literal alias for wB97X-D in the template should be the true + # functional's full LibXC name (AUDIT F03 fix), not bare wb97x. + # Doubled-brace literals in the template appear as single braces + # in the output. + assert "'wB97X-D': 'hyb_gga_xc_wb97x_d'" in PYSCF_SCRIPT_TEMPLATE assert "_NEEDS_D3" in PYSCF_SCRIPT_TEMPLATE - # The old (broken) "wb97x-d" string must NOT appear. + # Neither the black-listed short alias nor the wrong-functional + # bare-wb97x substitution should appear. assert "'wB97X-D': 'wb97x-d'" not in PYSCF_SCRIPT_TEMPLATE + assert "'wB97X-D': 'wb97x'" not in PYSCF_SCRIPT_TEMPLATE + # wB97X-D must not be wrapped in external D3 (its dispersion is + # already built into hyb_gga_xc_wb97x_d) — only PBE-D3 remains in + # _NEEDS_D3. Doubled braces are this template's literal-brace escape. + assert "_NEEDS_D3 = {{'PBE-D3'}}" in PYSCF_SCRIPT_TEMPLATE # =====================================================================