"""User-selectable J/K (Coulomb + exchange) method for periodic SCF, plus an intelligent AUTO picker. vibe-qc supports several alternatives to FFT-Poisson for the periodic two-electron Hartree-Fock / KS build. Each has a different sweet spot in (basis quality, lattice shape, system size, accuracy target). This module provides: * :class:`PeriodicJKMethod` -- enum naming each method. * :func:`pick_jk_method` -- heuristic that turns ``AUTO`` into a concrete choice given the system + basis + cell shape. * :func:`validate_jk_method` -- hard error / warning on illegal or risky combinations. * :func:`describe_jk_method` -- short human-readable label for logs. Status by method: Method Status --------------------- --------------------------------------------- GDF Native Gamma and multi-k RHF/RKS/UHF/UKS implementation. The earlier PySCF-backed spike is retired: PySCF and CRYSTAL are external reference programs, not in-process vibe-qc backends. AICCM2026DEV_B Independent finite Born--von Karman torus RHF/RKS/UHF/UKS. The finite-character construction is evaluated on its unreduced, Gamma-centred mesh; a real-Gamma form is Fourier-equivalent only for the same specified chi Hamiltonian. Experimental. BIPOLE ✓ multi-k RHF/UHF/RKS/UKS via CRYSTAL-gauge Ewald J-split: shared Ewald a, analytic J^LR, direct J^SR + K (HF) or V_xc (DFT), experimental opt-in multipole far-field branch off by default. DIRECT Limited scope. Plumbed against the existing C++ ``build_fock_2e_real_space`` which does the proper triple-cell-sum + Schwarz screening on each pair displacement. BUT -- with omega=0 (full Coulomb), the truncated real-space sum DIVERGES for tight ionic crystals (verified on MgO sto-3g: max|J-J_py| ≈ 363 Ha at cutoff 18 bohr). DIRECT is only valid for vacuum-padded (molecular-limit) cells where the cutoff actually clips the density. AUTO will refuse to pick DIRECT for non-vacuum-padded crystals. Native GDF/FFTDF is the production target. FFT_POISSON ✓ native EWALD_3D path. The public name is legacy; the default Hartree J backend now uses the analytical AO-pair FT with G=0 omitted. The historical FFT-Poisson grid backend is env-gated for parity tests. GPW ✓ Γ-only RHF/RKS/UHF/UKS via FFT-Poisson on a smooth real-space grid; multi-k pure-DFT RKS (M3e); forces, Hessian, DFT+U, smearing, D3-BJ dispersion. Experimental -- emits GAPWExperimentalWarning. See docs/design_periodic_gapw.md. GAPW ✓ All-electron Gaussian-augmented plane-wave (Lippert & Hutter 1999, doi:10.1007/978-3- 540-48272-4_3; Krack & Parrinello 2000, doi:10.1039/B001167N). Adds per-atom radial augmentation on top of GPW. Experimental -- guarded by GAPWExperimentalWarning. RSGDF not yet implemented (range-separated GDF; Ye & Berkelbach, J. Chem. Phys. 154, 131104 (2021), DOI 10.1063/5.0046617). RIJCOSX Γ-only RHF on the dedicated periodic COSX driver; true multi-k RHF/RKS/UHF/UKS via the native GDF loop with k_exchange='cosx'. CFMM not yet implemented (continuous fast-multipole; White et al. CPL 230, 8). The AUTO picker chooses based on: - **Lattice shape**: all implemented methods accept general 3D lattice vectors at the API level; method quality is still basis- and system-dependent. - **Basis compactness**: a "compact" basis (sto-3g, minimal) lets DIRECT converge; diffuse bases (cc-pVTZ, def2-TZVP) need GDF or RSGDF. - **System size**: huge supercells favor CFMM (when implemented). The default policy does not silently fall back to a method that is known to be scientifically unsafe for tight ionic crystals. Until native FFTDF/GDF lands, ``AUTO`` raises with an explicit explanation. """ from __future__ import annotations import enum import warnings from typing import Optional import numpy as np __all__ = [ "PeriodicJKMethod", "pick_jk_method", "validate_jk_method", "describe_jk_method", "is_orthorhombic", "resolve_jk_method_string", ] class PeriodicJKMethod(enum.Enum): """User-selectable periodic J/K builder. Use ``AUTO`` to let vibe-qc choose; pass any other value to force a specific method. The resolved method is logged in the output file's banner so reproducibility doesn't depend on the AUTO heuristic version. """ AUTO = "auto" GDF = "gdf" # native Gamma and multi-k Gaussian density fitting # --- CCM approaches and a retained same-H representation control ------- # aiccm2026dev-a / gamma / gamma-ccm: union-and-weight/Wigner--Seitz four-center # Γ-CCM construction. # real-gamma: k-free real-Γ representation of the neutral fitted-torus # control Hamiltonian. It is not evidence for the Γ-CCM construction. # chi: finite-translation-group χ-CCM construction. AICCM2026DEV_A = "aiccm2026dev-a" # Γ-CCM union-and-weight construction AICCM2026DEV_A_REAL_GAMMA = "aiccm2026dev-a-real-gamma" # neutral control AICCM2026DEV_B = "aiccm2026dev-b" # finite-torus RHF/RKS/UHF/UKS (χ-CCM) BIPOLE = "bipole" # CRYSTAL-gauge Ewald J-split RHF/UHF/RKS/UKS (multi-k); # multi-k; shared Ewald a; hybrids OK DIRECT = "direct" # partial -- see module docstring FFT_POISSON = "fft_poisson" # implemented; native FFT metric RIJCOSX = "rijcosx" # RI-J + periodic COSX-K; Gamma RHF + multi-k all spin cases # GPW / GAPW: v0.11.x periodic route (Lippert & Hutter 1997, # doi:10.1080/00268979709482119; Krack & Parrinello 2000, # doi:10.1039/B001167N). GPW = FFT-Poisson Hartree-J on a smooth # real-space grid; GAPW = all-electron via per-atom radial # augmentation on top of GPW. Both are experimental (guarded by # GAPWExperimentalWarning). See docs/design_periodic_gapw.md. GPW = "gpw" # Γ-only RHF/RKS/UHF/UKS; multi-k RKS GAPW = "gapw" # all-electron via per-atom radial augmentation # Vacuum-free 2D (slab) Coulomb -- the rigorous Parry / de Leeuw-Perram # gauge routed via lattice_opts.coulomb_method=SLAB_EWALD_2D. AUTO for # dim=2 closed-/open-shell (RHF/RKS/UKS). Not a bulk (dim=3) method. SLAB_EWALD_2D = "slab_ewald_2d" RSGDF = "rsgdf" # not yet implemented CFMM = "cfmm" # not yet implemented # These two selectors do not yet compute through ``run_periodic_job``. The # first names the union-and-weight Γ-CCM construction; the second is retained # for the neutral fitted-torus real-Gamma control and must not be interpreted as # another Γ-CCM construction route. Both remain deliberately outside # ``_IMPLEMENTED`` and fail closed with precise library-entry pointers. _UNWIRED_CCM_ROUTES: frozenset[PeriodicJKMethod] = frozenset( {PeriodicJKMethod.AICCM2026DEV_A, PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA} ) # Short / spoken aliases accepted anywhere a jk_method *string* is resolved, # normalised to the canonical enum value before lookup. ``gamma`` and # ``gamma-ccm`` name the Γ-CCM construction. ``neutral-bloch`` / ``gdf`` / # ``aiccm-ri`` name the neutral fitted-torus Bloch control, while # ``real-gamma`` names its Fourier-related real-supercell control. Keys are # lower-cased. _JK_STRING_ALIASES: dict[str, str] = { "real-gamma": PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA.value, "real_gamma": PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA.value, "gamma": PeriodicJKMethod.AICCM2026DEV_A.value, "neutral-bloch": PeriodicJKMethod.GDF.value, "neutral_bloch": PeriodicJKMethod.GDF.value, "bloch-control": PeriodicJKMethod.GDF.value, "gdf-control": PeriodicJKMethod.GDF.value, "aiccm-ri": PeriodicJKMethod.GDF.value, "gamma-ccm": PeriodicJKMethod.AICCM2026DEV_A.value, "gamma_ccm": PeriodicJKMethod.AICCM2026DEV_A.value, "chi": PeriodicJKMethod.AICCM2026DEV_B.value, "chi-ccm": PeriodicJKMethod.AICCM2026DEV_B.value, "slab": PeriodicJKMethod.SLAB_EWALD_2D.value, "ewald_2d": PeriodicJKMethod.SLAB_EWALD_2D.value, "ewald-2d": PeriodicJKMethod.SLAB_EWALD_2D.value, "slab-ewald": PeriodicJKMethod.SLAB_EWALD_2D.value, } _REJECTED_A_PREFIXED_CONTROL_STRINGS = frozenset( { PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA.value, "aiccm2026dev-a-direct", } ) def resolve_jk_method_string(name: str) -> "PeriodicJKMethod": """Public alias of the jk_method string resolver (aliases applied). ``"real-gamma"`` → :attr:`PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA`, ``"gamma"`` / ``"gamma-ccm"`` → :attr:`PeriodicJKMethod.AICCM2026DEV_A`, ``"neutral-bloch"`` / ``"aiccm-ri"`` → :attr:`PeriodicJKMethod.GDF`, and ``"chi"`` → :attr:`PeriodicJKMethod.AICCM2026DEV_B`. A-prefix spellings for the real-Gamma control are rejected because A names the Γ-CCM construction. """ return _resolve_jk_string(name) def _resolve_jk_string(name: str) -> "PeriodicJKMethod": """Resolve a jk_method *string* (with aliases) to a :class:`PeriodicJKMethod`. Applies :data:`_JK_STRING_ALIASES` before the enum lookup. In particular, ``"gamma"`` / ``"gamma-ccm"`` select the union-and-weight construction, while ``"neutral-bloch"`` and ``"aiccm-ri"`` select ordinary GDF as the neutral fitted-torus Bloch control. Raises :class:`ValueError` with the full valid list on an unknown name. """ key = name.strip().lower() if key in _REJECTED_A_PREFIXED_CONTROL_STRINGS: raise ValueError( f"Periodic JK method {name!r} is rejected because the " "'aiccm2026dev-a' prefix denotes the union-and-weight Γ-CCM " "construction, while real-Gamma is a neutral fitted-torus " "representation control. Use 'aiccm2026dev-a' for Γ-CCM or " "'real-gamma' for the control." ) key = _JK_STRING_ALIASES.get(key, key) try: return PeriodicJKMethod(key) except ValueError: valid = ", ".join( m.value for m in PeriodicJKMethod if m != PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA ) aliases = ", ".join(sorted(_JK_STRING_ALIASES)) raise ValueError( f"Unknown periodic JK method: {name!r}. Valid: {valid}. " f"Accepted aliases: {aliases}." ) # Methods that are actually wired up natively. _IMPLEMENTED: frozenset[PeriodicJKMethod] = frozenset( { PeriodicJKMethod.GDF, # Gamma and multi-k in periodic_runner PeriodicJKMethod.AICCM2026DEV_B, # finite BvK torus RHF/RKS/UHF/UKS PeriodicJKMethod.BIPOLE, # multi-k RHF/UHF/RKS/UKS via pbc_bipole PeriodicJKMethod.DIRECT, # only valid for vacuum-padded cells # PeriodicJKMethod.FFT_POISSON -- RETIRED v0.13.0 as a user route # (validate_jk_method raises). Internal Γ-only EWALD_3D drivers # remain for dilute periodic + mechanics, fail-closed on dense cells. PeriodicJKMethod.GPW, # M2-full / M3a / M3b -- Γ-only RHF PeriodicJKMethod.GAPW, # M3 -- adds Gaussian augmentation PeriodicJKMethod.RIJCOSX, # Gamma RHF + multi-k GDF/COSX PeriodicJKMethod.SLAB_EWALD_2D, # dim=2 vacuum-free slab (Γ + multi-k) # PeriodicJKMethod.RSGDF, # PeriodicJKMethod.CFMM, } ) # ============================================================ # Helpers # ============================================================ def is_orthorhombic(lattice: np.ndarray, tol: float = 1e-10) -> bool: """True iff ``lattice`` is diagonal (axis-aligned orthorhombic). Compatibility helper for older routing code and docs. ``lattice`` is the (3, 3) Cartesian matrix whose columns are a₁, a₂, a₃. """ L = np.asarray(lattice, dtype=float) off_diag = L - np.diag(np.diag(L)) diag_norm = max(np.linalg.norm(np.diag(L)), 1.0) return float(np.max(np.abs(off_diag))) < tol * diag_norm def _basis_is_compact(basis_name: str) -> bool: """Is this a compact (minimal / sto-3g-like) basis where DIRECT can converge with manageable cell-cutoff? Heuristic: name starts with "sto" (any zeta count) or contains "minimal". Could be extended once DIRECT is back online. """ n = (basis_name or "").lower() return n.startswith("sto") or "minimal" in n # ============================================================ # AUTO picker # ============================================================ def pick_jk_method( method: PeriodicJKMethod | str, *, lattice: np.ndarray, basis_name: str, n_atoms: int, scf_method: str = "RHF", dim: int = 3, ) -> PeriodicJKMethod: """Resolve an ``AUTO`` choice into a concrete :class:`PeriodicJKMethod`. Returns ``method`` unchanged if it is already concrete (and implemented). Raises if the resolved method is not implemented or is incompatible with the lattice / basis / dimensionality. ``dim`` is the system's periodic dimensionality (``system.dim``). For ``dim == 2`` (a slab), AUTO resolves to the vacuum-free :attr:`~PeriodicJKMethod.SLAB_EWALD_2D` gauge for RHF/RKS/UKS; every bulk (dim=3) J/K builder fails closed on a slab because running it would treat the slab as a 3-D crystal of sheets ``a3`` apart, giving ``a3``-/k-mesh- dependent energies (CLAUDE.md Sec. 7). For ``dim in (1, 3)`` AUTO keeps the established heuristic: GDF for closed-shell (RHF/RKS), BIPOLE for open-shell (UHF/UKS); explicit open-shell GDF is reachable via ``jk_method='gdf'``. """ if isinstance(method, str): method = _resolve_jk_string(method) dim = int(dim) spin = scf_method.upper() # --- dim=2 (slab) routing: the vacuum-free SLAB_EWALD_2D gauge only ------ if dim == 2: if method in (PeriodicJKMethod.AUTO, PeriodicJKMethod.SLAB_EWALD_2D): if spin in ("RHF", "RKS", "UKS"): return PeriodicJKMethod.SLAB_EWALD_2D raise NotImplementedError( "dim=2 (slab) SCF is available for RHF/RKS/UKS via the " "vacuum-free SLAB_EWALD_2D gauge; " f"scf_method={scf_method!r} is not yet supported for slabs " "(open-shell HF / UHF on slabs is a follow-on)." ) # An explicit bulk (dim=3) builder on a slab would silently run the # slab as a 3-D crystal of sheets a3 apart -> a3-/k-mesh-dependent # garbage. Fail closed with a pointer (CLAUDE.md Sec. 7). raise NotImplementedError( f"jk_method={method.value!r} is a bulk (dim=3) Coulomb builder and " "cannot run a dim=2 slab correctly: it would treat the slab as a " "3-D crystal of sheets a3 apart, giving a3-/k-mesh-dependent " "energies. Use jk_method='auto' (or 'slab_ewald_2d') for the " "rigorous vacuum-free 2D Coulomb." ) # SLAB_EWALD_2D is a dim=2-only gauge; reject it on bulk / polymer cells. if method == PeriodicJKMethod.SLAB_EWALD_2D and dim != 2: raise NotImplementedError( f"jk_method='slab_ewald_2d' requires a dim=2 slab; got dim={dim}." ) if method != PeriodicJKMethod.AUTO: return method # dim in (1, 3): established AUTO heuristic (dim=1 hits downstream raises). if spin in ("UHF", "UKS"): return PeriodicJKMethod.BIPOLE return PeriodicJKMethod.GDF # ============================================================ # Validation # ============================================================ def validate_jk_method( method: PeriodicJKMethod, *, lattice: np.ndarray, basis_name: str, ) -> None: """Raise / warn on illegal or risky combinations.""" if method == PeriodicJKMethod.FFT_POISSON: # Retired as a user-facing route (v0.13.0). The Γ-only EWALD_3D # path returns a ~2 Ha-wrong energy on dense ionic crystals (its # molecular-limit density convention D(g!=0)=0 is invalid once # periodic images overlap) and is superseded everywhere by GDF / # BIPOLE / GPW. The internal drivers (run_r{h,k}f_periodic_gamma_ # ewald3d) remain for dilute periodic + mechanics tests, fail-closed # on dense cells (CLAUDE.md Sec.7). raise ValueError( "jk_method='fft_poisson' (Γ-only EWALD_3D) is retired as of " "v0.13.0: it returns a ~2 Ha-wrong energy on dense ionic " "crystals and is superseded everywhere by GDF / BIPOLE / GPW. " "Use jk_method='gdf' (default), 'bipole', or 'gpw'." ) if method in _UNWIRED_CCM_ROUTES: # Fail closed with a route-specific pointer. Keep the Γ-CCM # construction separate from the neutral fitted-torus control. if method == PeriodicJKMethod.AICCM2026DEV_A: raise NotImplementedError( "jk_method='aiccm2026dev-a' (Γ-CCM) is the recognized " "union-and-weight/Wigner--Seitz four-center construction " "but is not yet wired into run_periodic_job (the FR-1 " "bundle follow-up). Call " "vibeqc.periodic.ccm.scf.run_ccm_rhf(ccm, " "method='aiccm2026dev-a') / " "vibeqc.periodic.ccm.dft.run_ccm_rks, or " "vibeqc.periodic.ccm.run_ccm_scf(ccm, route='four-center') " "(the default). The neutral-bloch/GDF and real-gamma routes are " "neutral fitted-torus representation controls, not Γ-CCM " "construction evidence." ) raise NotImplementedError( "jk_method='aiccm2026dev-a-real-gamma' is the retained k-free " "real-Γ representation control for the neutral fitted-torus " "Hamiltonian; it is not a Γ-CCM construction selector and is not " "yet wired into run_periodic_job " "(the FR-1 bundle follow-up; its per-supercell-Γ result needs the " "per-unit-cell result adapter). Run it today via " "vibeqc.periodic.ccm.direct.run_ccm_rhf_direct(ccm) or " "vibeqc.periodic.ccm.run_ccm_scf(ccm, route='real-gamma'); " "benchmark harness route 'aiccm-hf-direct'." ) # GPW / GAPW are wired through `run_periodic_job` as of M3b/M3c: # the GAPW dispatch uses the per-atom augmentation correction # for all-electron accuracy. Both are Γ-only RHF/UHF/RKS/UKS. if method == PeriodicJKMethod.GAPW: # GAPW is implemented -- see periodic_gapw_augment.py pass if method not in _IMPLEMENTED: raise NotImplementedError( f"Periodic JK method {method.value!r} is not yet wired up " f"natively. Currently available: " f"{ {m.value for m in _IMPLEMENTED} }" ) if method == PeriodicJKMethod.GDF: return if method == PeriodicJKMethod.AICCM2026DEV_B: return if method == PeriodicJKMethod.BIPOLE: # BIPOLE requires 3D (the Ewald gauge is 3D-only). # 1D/2D fall back to DIRECT_TRUNCATED unless the user # opts into another method. L = np.asarray(lattice, dtype=float) if np.linalg.matrix_rank(L) < 3: raise ValueError( "PeriodicJKMethod.BIPOLE requires a 3D lattice. " "For 1D/2D systems use GDF or DIRECT." ) return if method == PeriodicJKMethod.DIRECT: # build_fock_2e_real_space at omega=0 is real-space-truncated, # NOT Ewald-split. The truncated 1/r lattice sum diverges for # tight ionic crystals. Verified on MgO sto-3g cutoff=18: # max|J-J_pyscf| ≈ 363 Ha. DIRECT is only valid for cells # padded enough that the Coulomb tail outside cutoff is small. warnings.warn( "PeriodicJKMethod.DIRECT (build_fock_2e_real_space, omega=0) " "is only valid for VACUUM-PADDED (molecular-limit) cells. " "On tight ionic crystals the truncated real-space Coulomb " "lattice sum diverges (e.g., MgO sto-3g shows a historical " "hundreds-of-Hartree J error at cutoff 18 bohr). Native " "GDF/FFTDF is the target for tight crystals; DIRECT is " "appropriate for vacuum-padded molecules in PBC " "boxes (Makov-Payne regime). Future Method 2 augmentation " "would add a reciprocal-space LR contribution to make it " "valid for tight crystals (= Ewald composition).", stacklevel=2, ) # ============================================================ # Description # ============================================================ _DESCRIPTIONS = { PeriodicJKMethod.AUTO: "AUTO -- pick at runtime", PeriodicJKMethod.GDF: "GDF -- Gaussian density fitting (native Gamma and " "multi-k RHF/RKS/UHF/UKS; CCM aliases 'neutral-bloch'/'aiccm-ri' denote " "the neutral fitted-torus Bloch control, not the Γ-CCM construction; " "PySCF/CRYSTAL are external references only)", PeriodicJKMethod.AICCM2026DEV_A: "AICCM2026DEV-A / Γ-CCM -- " "union-and-weight/Wigner--Seitz four-center construction. Library entry: " "run_ccm_rhf(method='aiccm2026dev-a') / run_ccm_scf(route='gamma')", PeriodicJKMethod.AICCM2026DEV_A_REAL_GAMMA: "AICCM2026DEV-A real-gamma " "compatibility selector -- neutral fitted-torus same-H real-Γ representation " "control (not Γ-CCM construction evidence), with one real generalized " "eigenproblem and the derived exchange-q0 seam. Library entry: " "run_ccm_rhf_direct / run_ccm_scf(route='real-gamma')", PeriodicJKMethod.AICCM2026DEV_B: "AICCM2026DEV-B / χ-CCM (chi) -- independent finite " "translation-group-character BvK-torus RHF/RKS/UHF/UKS " "(4-center, RI, or RIJCOSX; " "exact Gamma-supercell / Gamma-centred character-mesh equivalence only for " "the same specified χ Hamiltonian; experimental)", PeriodicJKMethod.BIPOLE: "BIPOLE -- CRYSTAL-gauge Ewald J-split (multi-k; " "shared Ewald a across V_ne/E_nn/J^LR; " "all methods RHF/UHF/RKS/UKS; hybrids OK)", PeriodicJKMethod.DIRECT: "DIRECT -- 4-center periodic lattice sum + Schwarz screening " "(partial: missing-s_q bug)", PeriodicJKMethod.FFT_POISSON: "FFT_POISSON -- RETIRED v0.13.0 (Γ-only EWALD_3D; " "wrong on dense ionic crystals). Use GDF (default) / BIPOLE / GPW", PeriodicJKMethod.RIJCOSX: "RIJCOSX -- RI-J (density fitting) + periodic COSX-K " "(Gamma RHF plus true multi-k RHF/RKS/UHF/UKS; experimental -- periodic " "COSX is novel)", PeriodicJKMethod.GPW: "GPW -- Gaussian + plane-wave J via FFT-Poisson on a smooth " "real-space grid (Γ-only RHF/RKS/UHF/UKS; multi-k RKS; " "forces, Hessian, DFT+U, smearing, D3-BJ; experimental)", PeriodicJKMethod.GAPW: "GAPW -- Gaussian-augmented plane-wave J (Lippert-Hutter 1997 + " "Krack-Parrinello 2000); smooth plane-wave grid + per-atom " "radial augmentation for all-electron accuracy (experimental)", PeriodicJKMethod.RSGDF: "RSGDF -- range-separated GDF (not yet implemented)", PeriodicJKMethod.CFMM: "CFMM -- continuous fast multipole (not yet implemented)", } def describe_jk_method(method: PeriodicJKMethod) -> str: """Short human-readable description for output logs.""" return _DESCRIPTIONS.get(method, str(method))