"""High-level "run a periodic job" — companion to :func:`run_job`. Mirrors the molecular ``run_job`` API for periodic SCFs: - ``output.out`` text log (banner, system, basis, SCF trace, energies, properties) - ``output.system`` runtime manifest (CPU, OS, libs, wall-time) - ``output.molden`` Γ-point MOs in MOLDEN format (when MOs exist) - ``output.xsf`` SCF density on a primitive-cell grid (when ``write_density=True``) Usage:: from vibeqc import PeriodicSystem, BasisSet from vibeqc.periodic_runner import run_periodic_job sys_p = PeriodicSystem(...) basis = BasisSet(sys_p.unit_cell_molecule(), "sto-3g") run_periodic_job( sys_p, basis, method="RHF", output="output", ) Status: native Γ-only RHF/RKS GDF is available; multi-k GDF/FFTDF is under construction. The earlier PySCF-backed GDF spike is retired because PySCF and CRYSTAL are external reference programs, not in-process vibe-qc backends. Explicit ``jk_method="fft_poisson"`` still dispatches the legacy native EWALD_3D debug path where available; ``AUTO`` now resolves to native Γ-only GDF for RHF jobs. """ from __future__ import annotations import os import time from pathlib import Path from typing import TYPE_CHECKING, Optional, Union import numpy as np if TYPE_CHECKING: from .bands import BandStructure from ._vibeqc_core import ( Atom, BasisSet, Crystal, InitialGuess, PeriodicKSOptions, PeriodicRHFOptions, PeriodicSystem, SpaceGroup, attach_symmetry, to_primitive, ) from .banner import banner, library_versions from .occupations import ( hartree_to_kelvin_temperature, resolve_smearing_temperature, ) from .output import ( ManifestUpdater, OutputPlan, dry_run_manifest, is_dry_run_requested, write_extended_xyz, ) from .output import ( write_cif as _write_cif, ) from .output import ( write_population as _write_population, ) from .output import ( write_poscar as _write_poscar, ) from .output.citations import ( format_references_block, load_default_database, write_bibtex, write_references, ) from .periodic_jk_method import ( PeriodicJKMethod, describe_jk_method, pick_jk_method, validate_jk_method, ) from .periodic_k_gdf import ( run_krhf_periodic_gdf, run_krks_periodic_gdf, ) from .periodic_rhf_gdf import ( PeriodicRHFGDFResult, run_rhf_periodic_gamma_gdf, ) from .progress import ProgressLogger, resolve_progress __all__ = ["run_periodic_job"] # ============================================================ # Helpers — symmetry reduction # ============================================================ def _reduce_system_to_primitive( system: PeriodicSystem, *, symprec: float = 1e-4, ) -> tuple[PeriodicSystem, SpaceGroup]: """Reduce a ``PeriodicSystem`` to its primitive cell via spglib. Returns ``(primitive_system, space_group)`` where ``space_group`` is the symmetry analysis of the **original** cell. The primitive system carries its own symmetry analysis as well. Raises ``ValueError`` if the primitive cell is identical to the input (i.e. the input is already primitive), since no reduction is possible. """ # Attach symmetry to the original cell first. attach_symmetry(system, symprec=symprec) sg_original = system.symmetry if sg_original is None: raise RuntimeError("attach_symmetry did not populate system.symmetry") # Build a Crystal from the PeriodicSystem so we can call to_primitive. L = np.asarray(system.lattice, dtype=float, order="F") n_atoms_in = len(system.unit_cell) # Fractional coordinates: r_frac = L^{-1} · r_cart inv_L = np.linalg.inv(L) frac_coords = np.empty((3, n_atoms_in), dtype=float, order="F") species = [] for i, atom in enumerate(system.unit_cell): r_cart = np.array(atom.xyz, dtype=float) frac = inv_L @ r_cart frac_coords[:, i] = frac species.append(int(atom.Z)) crystal_in = Crystal(L, frac_coords, species) crystal_prim = to_primitive(crystal_in, symprec=symprec) n_atoms_out = crystal_prim.n_atoms if n_atoms_out == n_atoms_in: raise ValueError( "reduce_to_primitive: input cell is already primitive " f"({n_atoms_in} atoms, space group " f"{sg_original.international_symbol} " f"(No. {sg_original.number})). " "Set reduce_to_primitive=False to skip reduction." ) # Build a new PeriodicSystem from the primitive Crystal. prim_lattice = np.asarray(crystal_prim.lattice, dtype=float, order="F") prim_atoms = [] for col in range(n_atoms_out): r_cart = prim_lattice @ np.asarray( crystal_prim.fractional_coords[:, col], dtype=float ) prim_atoms.append(Atom(int(crystal_prim.species[col]), r_cart.tolist())) system_prim = PeriodicSystem( dim=system.dim, lattice=prim_lattice, unit_cell=prim_atoms, charge=system.charge, multiplicity=system.multiplicity, ) # Attach symmetry to the primitive cell too. attach_symmetry(system_prim, symprec=symprec) return system_prim, sg_original def _build_primitive_summary( system_in: PeriodicSystem, system_prim: PeriodicSystem, sg: SpaceGroup, ) -> str: """Human-readable summary of the cell reduction.""" n_in = len(system_in.unit_cell) n_out = len(system_prim.unit_cell) ratio = n_in / n_out if n_out > 0 else 1.0 L_in = np.asarray(system_in.lattice, dtype=float) L_out = np.asarray(system_prim.lattice, dtype=float) vol_in = float(abs(np.linalg.det(L_in))) vol_out = float(abs(np.linalg.det(L_out))) lines = [ " Cell reduction (reduce_to_primitive=True)", " " + "-" * 56, f" space group = {sg.international_symbol} " f"(No. {sg.number}), point group {sg.point_group}", f" symmetry order = {sg.order}", f" input cell = {n_in} atoms, volume = {vol_in:.4f} bohr³", f" primitive cell = {n_out} atoms, volume = {vol_out:.4f} bohr³", f" reduction factor = {ratio:.1f}× fewer atoms", ] if sg.equivalent_atoms: lines.append( f" inequivalent = {len(set(sg.equivalent_atoms))} atom type(s)" ) lines.append("") return "\n".join(lines) # ============================================================ # Helpers — text output sections # ============================================================ def _system_summary(system: PeriodicSystem) -> str: """Print lattice + atoms in bohr.""" L = np.asarray(system.lattice, dtype=float) dim = int(system.dim) if dim == 1: measure_label = "periodic length" measure_unit = "bohr" measure = float(np.linalg.norm(L[:, 0])) elif dim == 2: measure_label = "periodic area" measure_unit = "bohr^2" measure = float(np.linalg.norm(np.cross(L[:, 0], L[:, 1]))) else: measure_label = "cell volume" measure_unit = "bohr^3" measure = float(abs(np.linalg.det(L))) out = [] out.append(" Periodicity") out.append(" " + "-" * 56) out.append(f" dimensionality = {dim}D") active_axes = ", ".join(f"a{i + 1}" for i in range(dim)) out.append(f" active axes = {active_axes}") out.append(f" {measure_label:<14s} = {measure:.4f} {measure_unit}") out.append("") out.append(" Lattice (bohr)") out.append(" " + "-" * 56) for i in range(3): out.append(f" a{i + 1} = {L[i, 0]:14.8f} {L[i, 1]:14.8f} {L[i, 2]:14.8f}") if dim < 3: supercell_volume = float(abs(np.linalg.det(L))) out.append(f" embedding volume = {supercell_volume:.4f} bohr^3") out.append("") out.append(f" Atoms (bohr) — {len(system.unit_cell)} in unit cell") out.append(" " + "-" * 56) for i, atom in enumerate(system.unit_cell, start=1): x, y, z = atom.xyz out.append(f" {i:4d} Z={atom.Z:3d} {x:14.8f} {y:14.8f} {z:14.8f}") out.append( f" n_electrons = {system.n_electrons()} " f"multiplicity = {system.multiplicity}" ) out.append("") return "\n".join(out) def _basis_summary(basis: BasisSet) -> str: return ( " Basis\n" " " + "-" * 56 + "\n" f" name = {basis.name}\n" f" nbasis = {basis.nbasis}\n" f" nshells = {basis.nshells}\n" "\n" ) def _scf_trace_text(result, *, with_decomp: bool = False) -> str: """Render the SCF iteration trace from result.scf_trace.""" lines = [ " SCF iteration trace", " " + "-" * 70, f" {'iter':>4s} {'energy (Ha)':>18s} {'dE':>14s} " f"{'||[F,DS]||':>14s} {'DIIS':>5s}", ] for step in result.scf_trace: de = f"{step.delta_e:+.4e}" if step.iter > 1 else " --" lines.append( f" {step.iter:4d} {step.energy:18.10f} {de:>14s} " f"{step.grad_norm:14.4e} {step.diis_subspace:5d}" ) lines.append("") return "\n".join(lines) def _energy_summary(result, *, label: str) -> str: text = ( " Energy summary\n" " " + "-" * 56 + "\n" f" {'method':>20s} = {label}\n" f" {'converged':>20s} = {bool(result.converged)}\n" f" {'n_iter':>20s} = {int(result.n_iter)}\n" f" {'E_electronic (Ha)':>20s} = {float(result.e_electronic):20.10f}\n" f" {'E_nuclear (Ha)':>20s} = {float(result.e_nuclear):20.10f}\n" ) if hasattr(result, "e_coulomb"): text += f" {'E_coulomb (Ha)':>20s} = {float(result.e_coulomb):20.10f}\n" if hasattr(result, "e_xc") and float(getattr(result, "e_xc", 0.0)) != 0.0: text += f" {'E_xc (Ha)':>20s} = {float(result.e_xc):20.10f}\n" if hasattr(result, "e_hf_exchange"): text += ( f" {'E_HF_exchange (Ha)':>20s} = {float(result.e_hf_exchange):20.10f}\n" ) text += f" {'E_total (Ha)':>20s} = {float(result.energy):20.10f}\n" smearing_T = float(getattr(result, "smearing_temperature", 0.0)) if smearing_T > 0.0: text += ( f" {'kBT_smearing (Ha)':>20s} = {smearing_T:20.10f}\n" f" {'T_elec (K)':>20s} = " f"{hartree_to_kelvin_temperature(smearing_T):20.3f}\n" f" {'entropy S/kB':>20s} = " f"{float(getattr(result, 'entropy', 0.0)):20.10f}\n" f" {'free_energy (Ha)':>20s} = " f"{float(getattr(result, 'free_energy', result.energy)):20.10f}\n" f" {'fermi_level (Ha)':>20s} = " f"{float(getattr(result, 'fermi_level', 0.0)):20.10f}\n" ) text += "\n" return text def _mo_summary(result, n_show: int = 20) -> str: """Print HOMO/LUMO + nearest occupied/virtual MO energies.""" mo_energies = result.mo_energies # Multi-k: list of per-k arrays; use the first k-point's energies. if isinstance(mo_energies, (list, tuple)): if len(mo_energies) == 0: return "" eps = np.asarray(mo_energies[0]) occ_list = getattr(result, "occupations", None) occ = ( np.asarray(occ_list[0], dtype=float) if isinstance(occ_list, (list, tuple)) and len(occ_list) > 0 else np.empty(0) ) is_multi_k = True else: eps = np.asarray(mo_energies) occ = np.asarray(getattr(result, "occupations", np.empty(0)), dtype=float) is_multi_k = False if occ.shape == eps.shape and occ.size: occ_mask = occ > 1e-8 virt_mask = occ < 2.0 - 1e-8 homo_idx = int(np.where(occ_mask)[0][-1]) if np.any(occ_mask) else -1 lumo_idx = int(np.where(virt_mask)[0][0]) if np.any(virt_mask) else eps.size centre = max(homo_idx + 1, 0) if homo_idx >= 0 else min(lumo_idx, eps.size) else: # Fallback for older result objects without occupation metadata. n_occ = int(np.sum(eps < 0)) homo_idx = n_occ - 1 lumo_idx = n_occ centre = n_occ lines = [" Molecular orbital energies (Ha) — sorted low → high", " " + "-" * 56] if is_multi_k: n_k = len(mo_energies) lines.append(f" (showing k-point 0 of {n_k}; per-k MOs in output file)") half = n_show // 2 lo = max(0, centre - half) hi = min(eps.size, centre + half) for i in range(lo, hi): marker = " HOMO" if i == homo_idx else " LUMO" if i == lumo_idx else "" if occ.shape == eps.shape and occ.size: lines.append( f" {i + 1:4d} {eps[i]:18.10f} occ={occ[i]:8.5f}{marker}" ) else: lines.append(f" {i + 1:4d} {eps[i]:18.10f}{marker}") lines.append("") return "\n".join(lines) # ============================================================ # Main entry point # ============================================================ def _has_valid_mo_coeffs(result) -> bool: """Check whether result.mo_coeffs is non-empty (array or list).""" mc = getattr(result, "mo_coeffs", None) if mc is None: return False if isinstance(mc, (list, tuple)): return len(mc) > 0 and hasattr(mc[0], "size") and mc[0].size > 0 return hasattr(mc, "size") and mc.size > 0 def _gamma_proxy_for_multi_k(result) -> _GammaProxy: """Wrap a multi-k result to expose Γ-point (k=0) MOs for molden/etc.""" return _GammaProxy( mo_coeffs=result.mo_coeffs[0] if isinstance(result.mo_coeffs, (list, tuple)) else result.mo_coeffs, mo_energies=result.mo_energies[0] if isinstance(result.mo_energies, (list, tuple)) else result.mo_energies, occupations=result.occupations[0] if isinstance(getattr(result, "occupations", None), (list, tuple)) else getattr(result, "occupations", None), density=result.density[0] if isinstance(result.density, (list, tuple)) else result.density, overlap=getattr(result, "overlap", None), ) class _GammaProxy: """Duck-typed result wrapping Γ-point (k=0) of a multi-k result.""" __slots__ = ("mo_coeffs", "mo_energies", "occupations", "density", "overlap") def __init__(self, *, mo_coeffs, mo_energies, occupations, density, overlap): self.mo_coeffs = mo_coeffs self.mo_energies = mo_energies self.occupations = occupations self.density = density self.overlap = overlap def run_periodic_job( system: PeriodicSystem, basis: BasisSet, *, method: str = "RHF", functional: Optional[str] = None, jk_method: Union[str, "PeriodicJKMethod"] = "auto", aux_basis: Optional[str] = None, output: Union[str, os.PathLike] = "output", use_diis: bool = True, damping: float = 0.0, fmixing_percent: Optional[float] = None, smearing_temperature: Union[float, str, None] = 0.0, smearing_unit: str = "hartree", smearing_method: str = "fermi-dirac", smearing_metallic: Optional[bool] = None, smearing_band_gap_hartree: Optional[float] = None, diis_start_iter: int = 2, diis_subspace_size: int = 8, max_iter: int = 80, conv_tol_energy: float = 1e-7, initial_guess: str = "SAD", write_molden_file: bool = True, write_density: bool = False, density_spacing_bohr: float = 0.2, write_xyz_file: bool = True, write_poscar_file: bool = True, write_xsf_structure_file: bool = True, write_cif_file: bool = True, write_population_file: bool = True, citations: bool = True, dry_run: bool = False, record_hostname: bool = True, progress: Union[bool, ProgressLogger, None] = None, verbose: Optional[int] = None, # --- BIPOLE-specific --------------------------------------------- ewald_omega: Optional[float] = None, ewald_precision: float = 1e-8, use_oda: bool = False, oda_trust_lambda_max: float = 1.0, use_mom: bool = False, level_shift: float = 0.0, # --- Multi-k (GDF and BIPOLE) ----------------------------------- kpoints: Optional[Union[Tuple[int, int, int], List[int]]] = None, # --- Cell reduction ----------------------------------------------- reduce_to_primitive: bool = False, symmetry_precision: float = 1e-4, symmetry: Union[bool, str] = False, symmetry_stabilize: bool = False, # --- BIPOLE multipole far-field ---------------------------------- use_multipole_far_field: bool = False, multipole_l_max: int = 2, # --- Geometry optimization ---------------------------------------- optimize: bool = False, optimize_max_iter: int = 30, optimize_conv_tol_grad: float = 1e-4, optimize_cell: bool = False, # QVF visualisation archive (v1). output_qvf: bool = False, # Harmonic vibrational frequencies via finite-difference Hessian. hessian: bool = False, # Pre-computed band structure to include in QVF (vibe-view bands plot). band_structure: Optional["BandStructure"] = None, ): """Run a periodic SCF job and write the standard output files. Mirrors :func:`vibeqc.run_job` but for periodic systems. Γ-only RHF/RKS GDF is the default when no k-mesh is specified. Multi-k RHF/RKS GDF and UHF/UKS BIPOLE are available via ``kpoints``. Parameters ---------- system, basis Periodic system + AO basis. method ``"RHF"`` (only currently supported method). functional XC functional for ``method="RKS"``. output Path stem; produces ``{output}.out``, ``{output}.system``, ``{output}.molden``, ``{output}.xsf`` (when ``write_density``). band_structure Optional :class:`BandStructure` pre-computed by :func:`vibeqc.band_structure` (or ``_hcore``). When given together with ``output_qvf=True``, the band structure is embedded in the QVF archive so vibe-view can render an interactive Plotly band-structure plot. Compute it before calling this function — the same workflow used for matplotlib plotting with :func:`vibeqc.plot.band_structure_figure`. aux_basis Optional auxiliary basis for ``jk_method="gdf"``. If omitted, vibe-qc chooses the current native-GDF default for ``basis.name``. use_diis, damping, fmixing_percent, diis_start_iter, diis_subspace_size, max_iter, conv_tol_energy SCF controls forwarded to the periodic driver. ``fmixing_percent`` mirrors CRYSTAL's ``FMIXING`` keyword: the percentage of the previous Fock/KS matrix mixed into the matrix diagonalised on the next cycle. It is separate from density damping. ``smearing_temperature`` may be a numeric electronic ``k_B T`` (interpreted via ``smearing_unit``), ``"auto"``, ``"metal"``, ``"small-gap"``, ``"debug"``, ``"none"`` / ``"off"``, or ``None``. Only ``smearing_method="fermi-dirac"`` is implemented today. ``smearing_metallic`` and ``smearing_band_gap_hartree`` guide the conservative ``"auto"`` guess. initial_guess ``"SAD"`` (default) or ``"HCORE"``. write_molden_file Emit ``{output}.molden`` of the Γ-point MOs (using the unit-cell molecule + basis as the molecular target). write_density Emit ``{output}.xsf`` with the SCF density on a primitive-cell grid (XSF works for any lattice; cube is orthorhombic-only). density_spacing_bohr Grid spacing for the XSF density. Default 0.2 bohr. hessian When True, compute harmonic vibrational frequencies for the unit-cell molecule via finite-difference Hessian. Frequencies and IR intensities printed to .out and embedded in QVF for vibe-view. Default False. Cost: ~6N SCF evaluations for the unit cell. """ method_upper = method.upper() if method_upper not in ("RHF", "RKS", "UHF", "UKS"): raise NotImplementedError( f"run_periodic_job: method={method!r} not supported. " f"Supported: RHF, RKS, UHF, UKS." ) if method_upper in ("RHF", "UHF") and functional is not None: raise ValueError( f"functional={functional!r} given but method={method_upper}; " f"use method='RKS' or 'UKS' for KS calculations" ) if method_upper in ("RKS", "UKS") and functional is None: raise ValueError(f"method={method_upper!r} requires functional=...") # Resolve the J/K method (AUTO → concrete pick) and validate the # combination of method × lattice × basis. This is intentionally # done before any expensive setup so user errors fire early. resolved_jk = pick_jk_method( jk_method, lattice=np.asarray(system.lattice, dtype=float), basis_name=basis.name, n_atoms=len(system.unit_cell), scf_method=method_upper, ) validate_jk_method( resolved_jk, lattice=np.asarray(system.lattice, dtype=float), basis_name=basis.name, ) # --- Cell reduction (symmetry) ------------------------------------ system_original_info: Optional[str] = None # Resolve the 'symmetry' kwarg into a reduction flag. _sym_val = str(symmetry).lower() if isinstance(symmetry, str) else "" _do_reduce = ( reduce_to_primitive or symmetry is True or _sym_val in ("auto", "reduce") ) _do_attach = _do_reduce or _sym_val == "attach" if _do_attach and system.symmetry is None: attach_symmetry(system, symprec=symmetry_precision) if _do_reduce: try: system_prim, sg_input = _reduce_system_to_primitive( system, symprec=symmetry_precision, ) except ValueError: # Already primitive — just attach symmetry and continue if system.symmetry is None: attach_symmetry(system, symprec=symmetry_precision) else: if system.symmetry is None: attach_symmetry(system, symprec=symmetry_precision) system_original_info = _build_primitive_summary( system, system_prim, sg_input, ) basis = BasisSet(system_prim.unit_cell_molecule(), basis.name) system = system_prim fock_mixing = 0.0 if fmixing_percent is not None: fock_mixing = float(fmixing_percent) / 100.0 if not (0.0 <= fock_mixing < 1.0): raise ValueError( "run_periodic_job: fmixing_percent must be in [0, 100); " f"got {fmixing_percent}" ) output_stem = Path(os.fspath(output)) out_path = output_stem.with_suffix(".out") molden_path = output_stem.with_suffix(".molden") xsf_path = output_stem.with_suffix(".xsf") out_path.parent.mkdir(parents=True, exist_ok=True) # Dry-run short-circuit (Phase O5). Mirrors the molecular run_job # path: build the OutputPlan from current kwargs, write a one-shot # ``{output}.system`` with ``[outputs].status = "dry_run"``, print # the declared-artefacts summary, and return None without running # the SCF. Honours both the ``dry_run=True`` kwarg and the # ``VIBEQC_DRY_RUN=1`` env var that vq's submit pre-flight sets. if dry_run or is_dry_run_requested(): dry_plan = OutputPlan.from_run_job_kwargs( output=output_stem, method=method_upper, basis=basis.name, functional=functional, write_molden_file=write_molden_file, write_xyz=write_xyz_file, write_poscar=write_poscar_file, write_xsf_structure=write_xsf_structure_file, write_cif=write_cif_file, output_qvf=output_qvf, write_density_xsf=write_density, write_population=write_population_file, citations=citations, crash_dump=False, job_kind="periodic_scf", ) dry_run_manifest( dry_plan, record_hostname=record_hostname, ) return None plog = resolve_progress(progress, verbose=verbose) # Build the option object expected by the selected native driver. opts = ( PeriodicKSOptions() if method_upper in ("RKS", "UKS") else PeriodicRHFOptions() ) if method_upper in ("RKS", "UKS"): opts.functional = str(functional) opts.use_diis = bool(use_diis) opts.damping = float(damping) opts.fock_mixing = float(fock_mixing) smearing = resolve_smearing_temperature( smearing_temperature, unit=smearing_unit, method=smearing_method, metallic=smearing_metallic, band_gap_hartree=smearing_band_gap_hartree, n_electrons=system.n_electrons(), ) opts.smearing_temperature = float(smearing.temperature) opts.diis_start_iter = int(diis_start_iter) opts.diis_subspace_size = int(diis_subspace_size) opts.max_iter = int(max_iter) opts.conv_tol_energy = float(conv_tol_energy) # pybind11 enums expose ``__members__`` rather than supporting # ``Enum[name]`` subscripting directly. try: opts.initial_guess = InitialGuess.__members__[initial_guess.upper()] except KeyError as exc: valid = ", ".join(InitialGuess.__members__.keys()) raise ValueError( f"unknown initial_guess={initial_guess!r} (valid: {valid})" ) from exc label = f"{method_upper}" if functional: label = f"{label} / {functional}" t_job_start = time.perf_counter() hessian_result = None # populated by hessian=True; consumed by QVF writer with open(out_path, "w", encoding="utf-8", buffering=1) as f: # --- Banner --------------------------------------------------- f.write(banner() + "\n\n") libs = library_versions() f.write(f" Job: PERIODIC {label} basis={basis.name}\n") f.write(f" J/K method: {describe_jk_method(resolved_jk)}\n") if jk_method != "auto" and jk_method != PeriodicJKMethod.AUTO: f.write(f" (user-requested: {jk_method!r})\n") else: f.write(f" (resolved from AUTO)\n") f.write("\n") f.write(_system_summary(system)) if system_original_info is not None: f.write(system_original_info) elif system.symmetry is not None: sg = system.symmetry f.write( f" Symmetry: {sg.international_symbol} (No. {sg.number}), " f"point group {sg.point_group}, order {sg.order}\n\n" ) f.write(_basis_summary(basis)) f.write(" SCF options\n") f.write(" " + "-" * 56 + "\n") f.write(f" use_diis = {opts.use_diis}\n") f.write(f" damping = {opts.damping}\n") if fmixing_percent is not None: f.write(f" fmixing_percent = {float(fmixing_percent)}\n") if opts.smearing_temperature > 0.0 or smearing.source != "explicit": f.write(f" smearing_method = {smearing.method}\n") f.write(f" smearing_source = {smearing.source}\n") if smearing.reason: f.write(f" smearing_reason = {smearing.reason}\n") f.write(f" smearing_temperature = {opts.smearing_temperature}\n") if opts.smearing_temperature > 0.0: f.write( " smearing_temperature_K = " f"{hartree_to_kelvin_temperature(opts.smearing_temperature)}\n" ) f.write(f" diis_start_iter = {opts.diis_start_iter}\n") f.write(f" diis_subspace_size = {opts.diis_subspace_size}\n") f.write(f" max_iter = {opts.max_iter}\n") f.write(f" conv_tol_energy = {opts.conv_tol_energy}\n") f.write(f" initial_guess = {initial_guess.upper()}\n") if resolved_jk == PeriodicJKMethod.GDF: f.write(f" aux_basis = {aux_basis or ''}\n") f.write("\n") plog.banner(f"run_periodic_job PERIODIC {label} basis={basis.name}") plog.info(f"Output file: {out_path}") # --- SCF (dispatch on resolved jk_method) -------------------- t0 = time.perf_counter() if resolved_jk == PeriodicJKMethod.GDF: if method_upper in ("UHF", "UKS"): raise NotImplementedError( f"GDF path does not support {method_upper}. " f"Use jk_method='bipole' for open-shell periodic jobs." ) # Multi-k GDF dispatch when a k-mesh is requested. if kpoints is not None: plog.info(f"kmesh = {kpoints}") f.write(f" kpoints = {kpoints}\n") if method_upper == "RKS": result = run_krks_periodic_gdf( system, basis, kpoints, opts, functional=functional, aux_basis=aux_basis, fock_mixing=fock_mixing, progress=plog, ) else: result = run_krhf_periodic_gdf( system, basis, kpoints, opts, functional=None, aux_basis=aux_basis, fock_mixing=fock_mixing, progress=plog, ) else: result = run_rhf_periodic_gamma_gdf( system, basis, opts, functional=(functional if method_upper == "RKS" else None), aux_basis=aux_basis, fock_mixing=fock_mixing, symmetry_stabilize=symmetry_stabilize, progress=plog, ) elif resolved_jk == PeriodicJKMethod.FFT_POISSON: # Legacy EWALD_3D path. Known broken on ionic crystals # (warning already emitted by validate_jk_method). if method_upper != "RHF": raise NotImplementedError( f"FFT_POISSON path only wired for RHF in v0.7.1; got " f"method={method!r}. Native FFTDF/GDF is pending " "for KS." ) from . import run_rhf_periodic_gamma_ewald3d result = run_rhf_periodic_gamma_ewald3d( system, basis, opts, omega=0.5, progress=plog, ) elif resolved_jk == PeriodicJKMethod.DIRECT: # Dispatch DIRECT through the periodic_jk_direct wrapper. # AUTO never picks DIRECT (it diverges on tight ionic # crystals); the user has explicitly opted in here. # NOTE: this is a placeholder dispatch — there's no DIRECT # SCF driver yet that builds J + K via jk_via_direct each # iter. It's a wrapper one would call directly outside the # periodic_runner SCF loop. Track the SCF-driver port in # docs/design_native_gdf.md (DIRECT loop variant). raise NotImplementedError( "DIRECT SCF driver is not wired into run_periodic_job " "yet. Use vibeqc.periodic_jk_direct.jk_via_direct(...) " "for one-shot J/K builds. AUTO picks GDF for SCF; opt " "in to DIRECT only for vacuum-padded debug studies." ) elif resolved_jk == PeriodicJKMethod.BIPOLE: # CRYSTAL-gauge Ewald J-split — all four method flavours. from ._vibeqc_core import monkhorst_pack as _mp # Use user-provided k-mesh if given, else default to Γ-only if kpoints is not None: kp = ( list(kpoints) if isinstance(kpoints, (list, tuple)) else [kpoints, kpoints, kpoints] ) kmesh = _mp(system, list(kp)) plog.info(f" BIPOLE multi-k: mesh = {kp}") else: kmesh = _mp(system, [1, 1, 1]) if level_shift != 0.0: opts.level_shift = float(level_shift) if method_upper == "RHF": from .pbc_bipole import run_pbc_bipole_rhf result = run_pbc_bipole_rhf( system, basis, kmesh, opts, linear_dep_threshold=1e-7, use_ewald_j_split=True, ewald_omega=ewald_omega, ewald_precision=ewald_precision, use_oda=use_oda, oda_trust_lambda_max=oda_trust_lambda_max, use_mom=use_mom, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, progress=plog, ) elif method_upper == "UHF": from .pbc_bipole_uhf import run_pbc_bipole_uhf result = run_pbc_bipole_uhf( system, basis, kmesh, opts, linear_dep_threshold=1e-7, use_ewald_j_split=True, ewald_omega=ewald_omega, ewald_precision=ewald_precision, use_oda=use_oda, oda_trust_lambda_max=oda_trust_lambda_max, use_mom=use_mom, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, progress=plog, ) elif method_upper == "RKS": from .pbc_bipole_rks import run_pbc_bipole_rks result = run_pbc_bipole_rks( system, basis, kmesh, opts, functional=functional, linear_dep_threshold=1e-7, use_ewald_j_split=True, ewald_omega=ewald_omega, ewald_precision=ewald_precision, use_oda=use_oda, oda_trust_lambda_max=oda_trust_lambda_max, use_mom=use_mom, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, progress=plog, ) elif method_upper == "UKS": from .pbc_bipole_uks import run_pbc_bipole_uks result = run_pbc_bipole_uks( system, basis, kmesh, opts, functional=functional, linear_dep_threshold=1e-7, use_ewald_j_split=True, ewald_omega=ewald_omega, ewald_precision=ewald_precision, use_oda=use_oda, oda_trust_lambda_max=oda_trust_lambda_max, use_mom=use_mom, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, progress=plog, ) else: raise NotImplementedError( f"Periodic JK method {resolved_jk.value!r} dispatch " f"is not yet implemented in v0.7.1-spike." ) t_scf = time.perf_counter() - t0 # --- Write SCF trace + energies ------------------------------ f.write(_scf_trace_text(result)) f.write(_energy_summary(result, label=label)) f.write(_mo_summary(result)) # --- Molden --------------------------------------------------- if write_molden_file and _has_valid_mo_coeffs(result): try: from ._vibeqc_core import Molecule as _Mol from .io.molden import write_molden from .runner import _DispersionAugmented # noqa: F401 mol = system.unit_cell_molecule() # Multi-k: use Γ-point (k=0) MOs for molden output. molden_result = result if isinstance(result.mo_coeffs, (list, tuple)): molden_result = _gamma_proxy_for_multi_k(result) write_molden( molden_path, mol, basis, molden_result, title=output_stem.name ) f.write(f" Molecular orbitals written to {molden_path.name}\n") except Exception as exc: f.write( f" WARNING: molden write failed: {type(exc).__name__}: {exc}\n" ) # --- Population dump (Phase O6, periodic) --------------------- _qvf_pop = None if write_population_file and _has_valid_mo_coeffs(result): try: from ._vibeqc_core import Molecule as _Mol # noqa: F401 mol_p = system.unit_cell_molecule() pop_result = result if isinstance(result.mo_coeffs, (list, tuple)): pop_result = _gamma_proxy_for_multi_k(result) _write_population( output_stem, pop_result, basis, mol_p, ) f.write( f" Population dump written to " f"{output_stem.name}.population.txt + " f"{output_stem.name}.population.json\n" ) # Also compute summary for QVF atom_properties section. if output_qvf: try: from vibeqc.output.formats.population import ( compute_population_summary, ) _qvf_pop = compute_population_summary( result, basis, mol_p, ) except Exception: pass except Exception as exc: f.write( f" WARNING: population dump failed: {type(exc).__name__}: {exc}\n" ) # --- Geometry siblings (Phase O5) ----------------------------- # Extended-XYZ (lattice in comment line), POSCAR, and the # structure-only XSF — every viewer / chem-toolkit / fellow # QC code can consume at least one of these. Default-on with # opt-out via the matching kwarg; failures are best-effort # warnings so a finished SCF never gets dragged down by a # geometry writer. if write_xyz_file: xyz_path = output_stem.with_suffix(".xyz") try: energy_ha = float(getattr(result, "energy", float("nan"))) if energy_ha != energy_ha: # NaN sentinel energy_ha = None write_extended_xyz( output_stem, system, energy_ha=energy_ha, ) f.write( f" Final geometry written to {xyz_path.name} " f"(extended XYZ with lattice)\n" ) except Exception as exc: f.write( f" WARNING: extended-xyz write failed: " f"{type(exc).__name__}: {exc}\n" ) if write_poscar_file: poscar_path = output_stem.with_suffix(".POSCAR") try: _write_poscar( output_stem, system, comment=f"vibe-qc periodic {label} basis={basis.name}", ) f.write(f" Structure written to {poscar_path.name} (VASP-5 POSCAR)\n") except Exception as exc: f.write( f" WARNING: POSCAR write failed: {type(exc).__name__}: {exc}\n" ) if write_cif_file: cif_path = output_stem.with_suffix(".cif") try: _write_cif( output_stem, system, comment=f"vibe-qc periodic {label} basis={basis.name}", ) f.write(f" Structure written to {cif_path.name} (CIF, P 1)\n") except Exception as exc: f.write(f" WARNING: CIF write failed: {type(exc).__name__}: {exc}\n") # --- Citations (Phase O5b) ------------------------------------ # Assemble the citation list for this periodic job (software + # libint + libxc-if-DFT + basis-set + functional + DIIS + # spglib + ECP-if-used) and emit {stem}.bibtex / .references # siblings plus a "## References" block in the .out so the # text log is self-contained. Failures are non-fatal — the # SCF result is the load-bearing artefact, citations are # observability. cite_block_text: str | None = None if citations: try: _db = load_default_database() _refs = _db.assemble( method=method_upper, basis=basis.name, functional=functional, periodic=True, # ⇒ spglib fires uses_ecp=False, # ECP path lands in O5c ) write_bibtex(output_stem, _refs) write_references(output_stem, _refs) cite_block_text = format_references_block(_refs) bibtex_path = output_stem.with_suffix(".bibtex") f.write( f" Citations written to {bibtex_path.name} + " f"{output_stem.with_suffix('.references').name}\n" ) except Exception as exc: f.write( f" WARNING: citation emission failed: " f"{type(exc).__name__}: {exc}\n" ) if write_xsf_structure_file: xsf_struct_path = output_stem.with_suffix(".xsf") try: from .xsf import write_xsf_structure # When write_density=True we'd collide with the # density XSF below; route the structure-only XSF to # a different suffix in that case. if write_density: xsf_struct_path = output_stem.with_suffix( ".structure.xsf", ) write_xsf_structure(xsf_struct_path, system) f.write( f" Structure written to {xsf_struct_path.name} " f"(XSF crystal block)\n" ) except Exception as exc: f.write( f" WARNING: XSF structure write failed: " f"{type(exc).__name__}: {exc}\n" ) # --- Density XSF + grid for QVF ------------------------------- # Compute the primitive-cell density grid once; reuse it for # both the XSF writer and the QVF archive (when requested). # The grid origin is (0,0,0) and the span is the lattice- # vector transpose (rows = spanning vectors), matching the XSF # convention. All units are bohr — write_qvf stores them # as-is in the grid descriptor (see _grid_descriptor). _qvf_volume_data = None _needs_density_grid = write_density or output_qvf if _needs_density_grid: try: from ._vibeqc_core import ( CoulombMethod, LatticeMatrixSet, LatticeSumOptions, compute_overlap_lattice, ) from .periodic_density import evaluate_periodic_density_on_grid from .xsf import write_xsf_volume # Wrap the Γ-only D as a degenerate LatticeMatrixSet # (block 0 = D, others zero). Match the convention # used by the periodic XC builder. # Build the k-summed real-space density for multi-k results. lat_opts = LatticeSumOptions() lat_opts.coulomb_method = CoulombMethod.EWALD_3D lat_opts.cutoff_bohr = 18.0 D_set = compute_overlap_lattice(basis, system, lat_opts) # Multi-k: k-summed density = Σ_k w_k D(k). if isinstance(result.density, (list, tuple)): weights = getattr(result, "kpoint_weights", None) if weights is not None and len(weights) == len(result.density): D = sum( float(w) * np.real(np.asarray(d)) for w, d in zip(weights, result.density) ) else: D = sum(np.real(np.asarray(d)) for d in result.density) / len( result.density ) else: D = np.asarray(result.density) zero = np.zeros_like(D) for i in range(len(D_set)): D_set.set_block(i, D if i == 0 else zero) # Resolve grid shape (same rule as write_xsf_density). L_bohr = np.asarray(system.lattice, dtype=float) shape = tuple( max( 1, int( np.ceil(np.linalg.norm(L_bohr[:, i]) / density_spacing_bohr) ), ) for i in range(3) ) rho, _ = evaluate_periodic_density_on_grid( basis, system, D_set, grid_shape=shape, spacing_bohr=density_spacing_bohr, ) # --- XSF output --- if write_density: write_xsf_volume( xsf_path, system, data=rho, name=f"{output_stem.name}_density", ) f.write(f" Density written to {xsf_path.name}\n") # --- Package for QVF ---------------------------------- # QVF writes grid in bohr. Origin is (0,0,0) for a # primitive-cell grid. voxel_vectors are per-voxel # step vectors: lattice column / shape_i. if output_qvf: import numpy as np per_voxel = L_bohr.T / np.array(shape, dtype=float) _qvf_volume_data = { "Electron density": ( rho, np.zeros(3, dtype=float), per_voxel, ), } except Exception as exc: f.write( f" WARNING: density grid evaluation failed: " f"{type(exc).__name__}: {exc}\n" ) # --- Timing --------------------------------------------------- t_total = time.perf_counter() - t_job_start n_iter = int(getattr(result, "n_iter", 0)) iter_avg = (t_scf / n_iter) if n_iter > 0 else float("nan") f.write("\n Timings (wall clock, seconds)\n") f.write(" " + "-" * 56 + "\n") f.write(f" {'SCF total':<28s} {t_scf:12.3f}\n") f.write(f" {'SCF avg per iter':<28s} {iter_avg:12.3f} ({n_iter} iters)\n") f.write(f" {'Job total':<28s} {t_total:12.3f}\n") f.flush() # --- References block (Phase O5b) ----------------------------- # Embed the assembled citation list in the .out so the text # log is self-contained — a user reading the .out doesn't have # to chase the .bibtex / .references siblings to know what to # cite. Same content that lives in the .references file, hard- # wrapped to match the SCF-trace layout. if cite_block_text: f.write("\n" + cite_block_text + "\n") f.flush() # --- Harmonic vibrational analysis (finite-difference Hessian) --- if hessian: _scf_converged = bool(getattr(result, "converged", False)) if not _scf_converged: f.write( "\n ## Vibrational Frequencies\n" " " + "-" * 52 + "\n" " SKIPPED -- SCF did not converge.\n" ) f.flush() else: try: from ._vibeqc_core import ( RHFOptions, RKSOptions, UHFOptions, UKSOptions, ) from .hessian import ( HessianFDOptions, compute_hessian_fd, ir_intensities, ) if method_upper == "RHF": _hess_scf_opts = RHFOptions() elif method_upper == "UHF": _hess_scf_opts = UHFOptions() elif method_upper == "RKS": _hess_scf_opts = RKSOptions() _hess_scf_opts.functional = str(functional) elif method_upper == "UKS": _hess_scf_opts = UKSOptions() _hess_scf_opts.functional = str(functional) else: raise ValueError( f"Hessian not available for method={method_upper}" ) _hess_opts = HessianFDOptions( include_dipole_derivatives=True, ) _uc_mol = system.unit_cell_molecule() hessian_result = compute_hessian_fd( _uc_mol, basis.name, method=method_upper, scf_options=_hess_scf_opts, hessian_options=_hess_opts, ) f.write( f"\n ## Vibrational Frequencies\n" f" " + "-" * 52 + "\n" f" Finite-difference Hessian (unit cell)" f" (step = {_hess_opts.step_bohr:.3f} bohr," f" {hessian_result.n_displacements} displacements)\n" f" Imaginary modes: {hessian_result.imaginary_count}\n" f" Linear molecule: {hessian_result.is_linear}\n\n" ) _n_skip = 5 if hessian_result.is_linear else 6 _freqs = hessian_result.frequencies_cm1 _ir = None try: _ir = ir_intensities(hessian_result) except Exception as _ir_exc: f.write( f" (warning: IR intensities not available: " f"{type(_ir_exc).__name__}: {_ir_exc})\n" ) _header = " {:<6s} {:>10s}" + ( " {:>12s}" if _ir is not None else "" ) if _ir is not None: f.write( _header.format("Mode", "Freq/cm\u207b\u00b9", "IR/(km/mol)") + "\n" ) f.write(" " + "-" * 52 + "\n") for k in range(_n_skip, len(_freqs)): _label = f"{k - _n_skip + 1}" _freq = _freqs[k] _iri = _ir[k] if _freq < 0: _freq_str = f"{abs(_freq):.1f}i" else: _freq_str = f"{_freq:.1f}" f.write( f" {_label:<6s} {_freq_str:>10s} {_iri:>12.2f}\n" ) else: f.write(_header.format("Mode", "Freq/cm\u207b\u00b9") + "\n") f.write(" " + "-" * 52 + "\n") for k in range(_n_skip, len(_freqs)): _label = f"{k - _n_skip + 1}" _freq = _freqs[k] if _freq < 0: _freq_str = f"{abs(_freq):.1f}i" else: _freq_str = f"{_freq:.1f}" f.write(f" {_label:<6s} {_freq_str:>10s}\n") f.write("\n") f.flush() except Exception as _hess_exc: f.write( f"\n ## Vibrational Frequencies\n" f" " + "-" * 52 + "\n" f" FAILED: {type(_hess_exc).__name__}: {_hess_exc}\n" ) f.flush() hessian_result = None # --- System manifest with [plan] + [outputs] (Phase O5) ------------- # Replaces the bare write_system_manifest call with the # OutputPlan / ManifestUpdater pipeline so the periodic .system # carries the same declarative-plan + status-flag surface as the # molecular path. vq's submit pre-flight + status polling reads # the same fields whether the job is molecular or periodic. real_plan = OutputPlan.from_run_job_kwargs( output=output_stem, method=method_upper, basis=basis.name, functional=functional, write_molden_file=write_molden_file, write_xyz=write_xyz_file, write_poscar=write_poscar_file, write_xsf_structure=write_xsf_structure_file, write_density_xsf=write_density, write_population=write_population_file, citations=citations, crash_dump=False, output_qvf=output_qvf, job_kind="periodic_scf", ) manifest = ManifestUpdater( real_plan, record_hostname=record_hostname, wall_seconds=t_total, ) # --- QVF visualisation archive (v1) ---------------------------------- if output_qvf: try: from vibeqc.output.formats.qvf import write_qvf as _write_qvf _qvf_path = _write_qvf( output_stem, real_plan, system=system, result=result, method=method_upper, basis=basis.name, functional=functional, wall_seconds=t_total, volume_data=_qvf_volume_data, hessian_result=hessian_result, band_structure=band_structure, population_summary=_qvf_pop, ) try: manifest.mark_written(_qvf_path) except Exception: pass except Exception as _qvf_exc: import warnings warnings.warn(f"QVF write failed: {type(_qvf_exc).__name__}: {_qvf_exc}") # Record actually-written artefacts so [[outputs.files]] reflects # what landed on disk (best-effort: a writer that failed silently # leaves its row with written=false, which is the right signal # for vq + downstream tooling). for declared in real_plan.files: if declared.path.is_file(): try: manifest.mark_written(declared.path) except Exception: # mark_written hashes the file; a brittle write # shouldn't trip the wrap-up step. pass manifest.finish(wall_seconds=t_total) # ---- Geometry optimization --------------------------------------- opt_result = None if optimize and result.converged: from .bipole_optimize import relax_atoms, relax_full plog.info("") plog.banner("Geometry optimization") basis_name = basis.name # Reconstruct kmesh from the options used during SCF from ._vibeqc_core import monkhorst_pack as _mp if kpoints is not None: kp = ( list(kpoints) if isinstance(kpoints, (list, tuple)) else [kpoints, kpoints, kpoints] ) km = _mp(system, list(kp)) else: km = _mp(system, [1, 1, 1]) if optimize_cell and system.dim == 3: from .bipole_optimize import relax_cell_gradient # First relax atoms, then cell (gradient-based), then atoms again opt_result = relax_atoms( system, basis_name, km, method.upper(), functional=functional, max_iter=optimize_max_iter, conv_tol_grad=optimize_conv_tol_grad, cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)), ewald_precision=ewald_precision, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, ) opt_result = relax_cell_gradient( opt_result.system, basis_name, km, method.upper(), functional=functional, max_iter=10, cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)), ewald_precision=ewald_precision, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, ) opt_result = relax_atoms( opt_result.system, basis_name, km, method.upper(), functional=functional, max_iter=optimize_max_iter, conv_tol_grad=optimize_conv_tol_grad, cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)), ewald_precision=ewald_precision, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, ) else: opt_result = relax_atoms( system, basis_name, km, method.upper(), functional=functional, max_iter=optimize_max_iter, conv_tol_grad=optimize_conv_tol_grad, cutoff_bohr=float(getattr(opts.lattice_opts, "cutoff_bohr", 8.0)), ewald_precision=ewald_precision, use_multipole_far_field=use_multipole_far_field, multipole_l_max=multipole_l_max, ) if opt_result is not None: opt_sys = opt_result.system with open(out_path, "a", encoding="utf-8", buffering=1) as f: f.write("\n Optimized geometry\n") f.write(" " + "-" * 56 + "\n") f.write(f" E_final = {opt_result.energy:+.10f} Ha\n") f.write(f" n_iter = {opt_result.n_iter}\n") f.write(f" converged = {opt_result.converged}\n") f.write("\n") f.write(_system_summary(opt_sys)) f.flush() # Write optimized geometry files if write_xyz_file: try: write_extended_xyz( output_stem.with_suffix(".opt"), opt_sys, energy_ha=opt_result.energy, ) except Exception: pass if write_poscar_file: try: _write_poscar( output_stem.with_suffix(".opt"), opt_sys, comment=f"vibe-qc optimized {label}", ) except Exception: pass result = opt_result # Return optimization result return result