"""Phase G1-GDF -- analytic Γ-only periodic GDF (compcell) atomic gradient. The GDF two-electron energy at Γ with density D: E_2e = 1/2 tr(D . J) - 1/4 a_HF . tr(D . K) where J and K are built from the compcell cderi L via :func:`vibeqc.pbc_gdf._build_j_from_lpq` / :func:`_build_k_from_lpq`. **J gradient (Coulomb).** The DF-J gradient formula: dE_J/dR = S_P g̃_P S_muν D_muν dT[P,muν]/dR - 1/2 S_PQ g̃_P g̃_Q dM_PQ/dR where g̃ = M^+ r, r_P = S_muν T[P,muν] D_muν, and (M, T) are the compensated 2c metric and 3c tensor. ``M^+`` is the exact thresholded pseudoinverse used to build the SCF cderi. Its derivative includes the Fréchet response of the retained eigenspace, not only the full-rank ``-g̃ g̃^T`` term. The compensation matrix A maps the fused basis (modrho-aux + compensating charges) to the physical aux space: M = A.M_fused.A^T, T = A.T_fused. The gradient maps back onto the fused basis via g̃_fused = A^T.g̃. **K gradient (exchange).** The DF-K gradient: dE_K/dR = +a_HF S_PQ w_PQ dM_PQ/dR - 2 a_HF S_{P,muν} Y^P_{muν} d(P|muν)/dR with w = (η^P : η^Q), η = M^+ K_occ, K^P_occ,ij = C_occ^T T^P C_occ, Y^P = C_occ . η^P . C_occ^T. For pure DFT (a_HF = 0) this term is zero. When the fit drops metric modes, the J and K metric weights use the spectral divided differences of the thresholded inverse; the displayed outer-product forms are their full-rank limit. **Current scope:** Gamma-point closed-shell RHF with the compcell cderi and without the AFT correction. The J and K derivatives use the C++ lattice-summed gradient kernels below. The FT-based Ewald ``V_ne`` derivative includes its real-space, structure-factor, AO-centre, and ``G = 0`` terms. Production RSGDF, RKS/open-shell, multi-k, and optimizer wiring remain separate increments. - :func:`vibeqc._vibeqc_core.compute_2c_eri_lattice_gradient_weighted` - :func:`vibeqc._vibeqc_core.compute_3c_eri_lattice_gradient_weighted` """ from __future__ import annotations import hashlib from dataclasses import dataclass from typing import Optional, Tuple import numpy as np from ._vibeqc_core import ( BasisSet, CoulombMethod, LatticeSumOptions, PeriodicSystem, compute_2c_eri_lattice_gradient_weighted, compute_3c_eri_lattice_gradient_weighted, compute_overlap_lattice, kinetic_lattice_gradient_contribution, nuclear_repulsion_gradient_per_cell, overlap_lattice_gradient_contribution, ) from .aux_basis import ( _CompcellFitState, _basis_fingerprint, _build_lpq_compcell_state, ) __all__ = ["compute_gdf_gradient_rhf_gamma", "compute_gdf_gradient"] # --------------------------------------------------------------------------- # Helpers # --------------------------------------------------------------------------- def _gamma_density_lattice_set(template_lat, D: np.ndarray) -> np.ndarray: """Build a Γ-only lattice-resolved density. For the Γ-point, the real-space density is identical in every image cell (k=0 Bloch phase = 1).""" # Actually, for the GDF route at Γ, the density is the Γ-folded # single-cell density. The lattice-resolved D(g) = D_Γ for all g. # But for the gradient, the ERI lattice gradient uses a # LatticeMatrixSet. Since we compute the 2e gradient through the # DF path (not through the 4c ERI path), we don't need D(g). # The D passed to the 3c gradient kernel is the single-cell density. return np.asarray(D, dtype=np.float64) @dataclass class _CompcellGradientCache: """Cached intermediates from the compcell pipeline for gradient use. The compensation matrix A and the fused basis integrals M_fused and T_fused depend only on the geometry, not on the density. They are computed once and reused for every gradient evaluation. """ A: np.ndarray # (n_aux, n_fused) compensation matrix M_fused: np.ndarray # (n_fused, n_fused) pre-compensation 2c metric T_fused: np.ndarray # (n_fused, n_orb, n_orb) pre-compensation 3c tensor fused_basis: BasisSet # the fused basis (modrho-aux + chg) lat_opts_2c: LatticeSumOptions lat_opts_3c: LatticeSumOptions eigvals: np.ndarray eigvecs: np.ndarray keep_mask: np.ndarray inverse_eigvals: np.ndarray inverse_frechet: np.ndarray linear_dep_thr: float system_fingerprint: str ao_basis_fingerprint: str aux_basis_fingerprint: str compcell_eta: float lattice_cutoff_bohr: float nuclear_cutoff_bohr: float rcut_strategy: Optional[str] rcut_precision: float n_fused: int n_aux: int n_orb: int n_fit: int def _cache_from_compcell_fit_state( state: _CompcellFitState, *, system_fingerprint: str, ao_basis_fingerprint: str, aux_basis_fingerprint: str, compcell_eta: float, lattice_cutoff_bohr: float, nuclear_cutoff_bohr: float, rcut_strategy: Optional[str], rcut_precision: float, ) -> _CompcellGradientCache: """Prepare the exact thresholded inverse response used by the SCF fit.""" eigvals = np.asarray(state.eigvals, dtype=np.float64) eigvecs = np.asarray(state.eigvecs, dtype=np.float64) keep = np.asarray(state.keep_mask, dtype=bool) if eigvals.ndim != 1 or eigvecs.shape != (eigvals.size, eigvals.size): raise ValueError("compcell gradient: malformed metric eigensystem") if keep.shape != eigvals.shape or not np.any(keep): raise ValueError("compcell gradient: malformed retained-mode mask") max_eig = float(eigvals[-1]) threshold = float(state.linear_dep_thr) * max_eig eig_scale = max(abs(max_eig), 1.0) spectral_tol = 128.0 * np.finfo(np.float64).eps * eig_scale if np.any(np.abs(eigvals - threshold) <= spectral_tol): raise RuntimeError( "compcell gradient: an auxiliary metric eigenvalue lies on the " "linear-dependence threshold, so the fit derivative is undefined" ) inverse_eigvals = np.zeros_like(eigvals) inverse_eigvals[keep] = 1.0 / eigvals[keep] # Divided differences for the Fréchet derivative of the thresholded # inverse f(M), f(lambda)=1/lambda on kept modes and zero otherwise. inverse_frechet = np.zeros((eigvals.size, eigvals.size), dtype=np.float64) kept_pair = keep[:, None] & keep[None, :] inverse_frechet[kept_pair] = -( inverse_eigvals[:, None] * inverse_eigvals[None, :] )[kept_pair] mixed_pair = keep[:, None] ^ keep[None, :] eig_diff = eigvals[:, None] - eigvals[None, :] if np.any(np.abs(eig_diff[mixed_pair]) <= spectral_tol): raise RuntimeError( "compcell gradient: retained and dropped auxiliary modes are " "numerically degenerate, so the fit derivative is undefined" ) value_diff = inverse_eigvals[:, None] - inverse_eigvals[None, :] inverse_frechet[mixed_pair] = ( value_diff[mixed_pair] / eig_diff[mixed_pair] ) n_aux = int(state.A.shape[0]) n_fused = int(state.A.shape[1]) n_orb = int(state.T_fused.shape[1]) return _CompcellGradientCache( A=np.asarray(state.A, dtype=np.float64), M_fused=np.asarray(state.M_fused, dtype=np.float64), T_fused=np.asarray(state.T_fused, dtype=np.float64), fused_basis=state.fused_basis, lat_opts_2c=state.lat_opts_2c, lat_opts_3c=state.lat_opts_3c, eigvals=eigvals, eigvecs=eigvecs, keep_mask=keep, inverse_eigvals=inverse_eigvals, inverse_frechet=inverse_frechet, linear_dep_thr=float(state.linear_dep_thr), system_fingerprint=system_fingerprint, ao_basis_fingerprint=ao_basis_fingerprint, aux_basis_fingerprint=aux_basis_fingerprint, compcell_eta=float(compcell_eta), lattice_cutoff_bohr=float(lattice_cutoff_bohr), nuclear_cutoff_bohr=float(nuclear_cutoff_bohr), rcut_strategy=rcut_strategy, rcut_precision=float(rcut_precision), n_fused=n_fused, n_aux=n_aux, n_orb=n_orb, n_fit=int(np.count_nonzero(keep)), ) def _system_fingerprint(system: PeriodicSystem) -> str: """Return a stable fingerprint for gradient-relevant cell geometry.""" digest = hashlib.sha256(b"vibeqc-periodic-gradient-system-v1\0") metadata = (int(system.dim), int(system.charge), int(system.multiplicity)) digest.update(":".join(str(value) for value in metadata).encode("ascii")) digest.update(b"\0") digest.update( np.asarray(system.lattice, dtype=" Optional[str]: if strategy is None: return None return str(getattr(strategy, "value", strategy)) def _build_compcell_gradient_cache( system: PeriodicSystem, ao_basis: BasisSet, aux_basis: BasisSet, compcell_eta: float, lattice_opts: LatticeSumOptions, *, linear_dep_thr: float = 1e-9, rcut_strategy: Optional[object] = None, rcut_precision: float = 1e-8, fit_state: Optional[_CompcellFitState] = None, ) -> _CompcellGradientCache: """Build the geometry-dependent (density-independent) cache for compcell GDF gradient evaluation. Uses the same private fit builder as the SCF so the resolved 2c/3c cutoffs, compensated integrals, threshold, and retained eigenspace are identical. ``fit_state`` avoids rebuilding those intermediates when the gradient is requested directly by the SCF driver. """ if fit_state is None: fit_state = _build_lpq_compcell_state( system, ao_basis, aux_basis, molecule=system.unit_cell_molecule(), lat_opts=lattice_opts, linear_dep_thr=float(linear_dep_thr), compcell_eta=float(compcell_eta), apply_aft_correction=False, rcut_strategy=rcut_strategy, rcut_precision=float(rcut_precision), ) return _cache_from_compcell_fit_state( fit_state, system_fingerprint=_system_fingerprint(system), ao_basis_fingerprint=_basis_fingerprint(ao_basis), aux_basis_fingerprint=_basis_fingerprint(aux_basis), compcell_eta=float(compcell_eta), lattice_cutoff_bohr=float(lattice_opts.cutoff_bohr), nuclear_cutoff_bohr=float(lattice_opts.nuclear_cutoff_bohr), rcut_strategy=_rcut_strategy_value(rcut_strategy), rcut_precision=float(rcut_precision), ) def _validated_result_madelung( system: PeriodicSystem, result, *, caller: str, ) -> float: """Validate and return the exchange-divergence constant used by SCF.""" exxdiv = str(getattr(result, "exxdiv", "")) retained = float(getattr(result, "madelung_constant", float("nan"))) if exxdiv not in {"ewald", "none"} or not np.isfinite(retained): raise ValueError(f"{caller}: invalid exxdiv/Madelung provenance") if exxdiv == "none": if abs(retained) > 1e-15: raise ValueError( f"{caller}: exxdiv='none' requires a zero retained " "Madelung constant" ) return retained from .madelung import madelung_constant_for_cell expected = float(madelung_constant_for_cell(system)) if retained <= 0.0 or not np.isclose( retained, expected, rtol=1e-12, atol=1e-14, ): raise ValueError( f"{caller}: retained Ewald Madelung constant does not match " "the periodic cell" ) return retained def _validate_supported_gradient_result(result, *, caller: str) -> None: """Fail closed unless ``result`` is from the implemented RHF route.""" if not bool(getattr(result, "converged", False)): raise ValueError(f"{caller}: SCF result is not converged") backend = str(getattr(result, "backend", "")) functional = str(getattr(result, "functional", "")) if backend != "pbc-gdf-compcell" or functional: raise NotImplementedError( f"{caller}: only a Gamma RHF compcell result is supported; " f"got backend={backend!r}, functional={functional!r}." ) if bool(getattr(result, "apply_aft_correction", True)): raise NotImplementedError( f"{caller}: apply_aft_correction=False is required until the " "compcell AFT metric derivative is implemented." ) if str(getattr(result, "v_ne_backend", "")) != "analytic_ft": raise NotImplementedError( f"{caller}: the converged result must use the analytical FT " "Ewald V_ne backend." ) def _compute_j_gradient_compcell( system: PeriodicSystem, ao_basis: BasisSet, D: np.ndarray, cache: _CompcellGradientCache, ) -> np.ndarray: """Compute the DF-J (Coulomb) analytic gradient via the compcell path. Uses the standard DF gradient formula mapped onto the fused basis via the compensation matrix A. Parameters ---------- system, ao_basis, D Periodic system, orbital basis, and converged density (n_orb, n_orb). cache Pre-built geometry cache from :func:`_build_compcell_gradient_cache`. Returns ------- grad_J : np.ndarray of shape (n_atoms, 3) in Hartree/bohr. """ n_aux = cache.n_aux n_fused = cache.n_fused n_orb = cache.n_orb # Step 1: Compensate T. The metric eigensystem in the cache was built # from M_comp = A @ M_fused @ A.T by the exact SCF fit route. A = cache.A T_comp = np.einsum("iP,Pmn->imn", A, cache.T_fused, optimize=True) # Step 2: Contract T with D to get r. # r_P = S_muν T_comp[P, mu, ν] . D[mu, ν] T_flat = T_comp.reshape(n_aux, n_orb * n_orb) D_flat = np.asarray(D, dtype=np.float64).ravel() rho = T_flat @ D_flat # (n_aux,) # Step 3: Apply the exact thresholded inverse used by the SCF cderi. U = cache.eigvecs rho_eig = U.T @ rho gamma_tilde = U @ (cache.inverse_eigvals * rho_eig) # Step 4: Project g̃ back onto the fused basis. # g̃_fused = A^T . g̃ (n_fused,) gamma_tilde_fused = A.T @ gamma_tilde # (n_fused,) # Step 5: Differentiate the thresholded inverse. The mixed retained / # dropped blocks are the moving-projector response and cannot be written # as only -1/2 g g.T when the fit dropped a mode. omega_eig = 0.5 * cache.inverse_frechet * np.outer(rho_eig, rho_eig) omega_aux = U @ omega_eig @ U.T omega_fused = A.T @ omega_aux @ A # Step 6: Build the 3c weight W. # W[P, muν] = g̃_fused[P] . D[muν] (n_fused, n_orb x n_orb) row-major W = np.outer(gamma_tilde_fused, D_flat) # (n_fused, n_orb^2) W = np.ascontiguousarray(W) # ensure row-major for C++ # Step 7: Call the C++ lattice-summed gradient kernels. grad_2c = np.asarray( compute_2c_eri_lattice_gradient_weighted( cache.fused_basis, system, cache.lat_opts_2c, omega_fused ), dtype=np.float64, ) grad_3c = np.asarray( compute_3c_eri_lattice_gradient_weighted( ao_basis, cache.fused_basis, system, cache.lat_opts_3c, W ), dtype=np.float64, ) return grad_2c + grad_3c def _compute_k_gradient_compcell( system: PeriodicSystem, ao_basis: BasisSet, D: np.ndarray, C_occ: np.ndarray, alpha_hf: float, cache: _CompcellGradientCache, ) -> np.ndarray: """Compute the DF-K (exchange) analytic gradient via the compcell path. For pure DFT (alpha_hf = 0) returns zero. Parameters ---------- system, ao_basis, D, C_occ, alpha_hf Same conventions as the molecular DF-K gradient. cache Same as :func:`_compute_j_gradient_compcell`. Returns ------- grad_K : np.ndarray of shape (n_atoms, 3) in Hartree/bohr. """ if alpha_hf == 0.0: return np.zeros((len(system.unit_cell), 3), dtype=np.float64) n_aux = cache.n_aux n_fused = cache.n_fused n_orb = cache.n_orb n_occ = C_occ.shape[1] A = cache.A T_comp = np.einsum("iP,Pmn->imn", A, cache.T_fused, optimize=True) # Step 1: Build M^P_ij = C_occ^T . T^P . C_occ (K_occ) # M_occ[P, i, j] for each aux function P. M_occ = np.einsum("Pmn,mi,nj->Pij", T_comp, C_occ, C_occ, optimize=True) # M_occ shape: (n_aux, n_occ, n_occ) # Step 2: Apply the SCF's thresholded inverse column-wise. M_occ_flat = M_occ.reshape(n_aux, n_occ * n_occ) U = cache.eigvecs M_occ_eig = U.T @ M_occ_flat eta_flat = U @ (cache.inverse_eigvals[:, None] * M_occ_eig) eta = eta_flat.reshape(n_aux, n_occ, n_occ) # Step 3: Exact metric weight, including the retained/dropped projector # response. For a full-rank fit this reduces to +alpha * eta @ eta.T. rhs_gram_eig = M_occ_eig @ M_occ_eig.T omega_eig = -float(alpha_hf) * cache.inverse_frechet * rhs_gram_eig omega_aux = U @ omega_eig @ U.T omega_fused = A.T @ omega_aux @ A # Step 4: Build Y^P_{muν} = (C_occ . η^P . C_occ^T)_{muν}. Y = np.einsum("Pij,mi,nj->Pmn", eta, C_occ, C_occ, optimize=True) # Y shape: (n_aux, n_orb, n_orb) Y_fused = np.einsum("Pi,Pmn->imn", A, Y, optimize=True) # Y_fused shape: (n_fused, n_orb, n_orb) Y_flat = Y_fused.reshape(n_fused, n_orb * n_orb) Y_flat = np.ascontiguousarray(Y_flat) # Step 5: Gradient contract. # dE_K/dR = +a_HF S_PQ w_PQ dM_PQ/dR - 2 a_HF S_{P,muν} Y^P_{muν} d(P|muν)/dR grad_2c = np.asarray( compute_2c_eri_lattice_gradient_weighted( cache.fused_basis, system, cache.lat_opts_2c, omega_fused ), dtype=np.float64, ) grad_3c = np.asarray( compute_3c_eri_lattice_gradient_weighted( ao_basis, cache.fused_basis, system, cache.lat_opts_3c, Y_flat ), dtype=np.float64, ) return grad_2c - 2.0 * alpha_hf * grad_3c # --------------------------------------------------------------------------- # Public API # --------------------------------------------------------------------------- def compute_gdf_gradient_rhf_gamma( system: PeriodicSystem, basis: BasisSet, result, # PBCGDFResult *, aux_basis: BasisSet, compcell_eta: Optional[float] = None, alpha_hf: float = 1.0, lattice_opts: Optional[LatticeSumOptions] = None, cache: Optional[_CompcellGradientCache] = None, madelung: Optional[float] = None, gauge_lat_opts: Optional[LatticeSumOptions] = None, v_ne_ke_cutoff: Optional[float] = None, ) -> Tuple[np.ndarray, _CompcellGradientCache]: """Analytic Γ-only GDF (compcell) RHF atomic gradient. Parameters ---------- system, basis Periodic system and AO basis. result Converged :class:`PBCGDFResult` from :func:`vibeqc.run_pbc_gdf_rhf`. aux_basis Auxiliary basis (unscaled -- the gradient pipeline handles modrho rescaling internally, matching the SCF). compcell_eta Smooth-Gaussian exponent for the compensating charges. ``None`` uses the value retained on ``result``; an explicit value must match it. alpha_hf HF-exchange fraction. The supported RHF route requires ``1.0``. lattice_opts Lattice-sum options. Must match the SCF. If None, reconstructed from the scalar cutoffs retained on ``result``. cache Optional pre-built gradient cache from the identical system, AO and auxiliary bases, compensation exponent, and cutoff policy. If None, built fresh; incompatible provenance is rejected. madelung Ewald exchange-divergence constant. If None, uses the value retained on ``result``; an explicit value must match it. v_ne_ke_cutoff Reciprocal kinetic-energy cutoff used by the SCF's Ewald ``V_ne``. ``None`` uses the value retained on ``result``; an explicit value must match it. Returns ------- grad : np.ndarray of shape (n_atoms, 3) in Hartree/bohr. cache : _CompcellGradientCache The cache built (or passed in). Caller can reuse across multiple gradient evaluations at the same geometry. """ if int(system.dim) != 3: raise NotImplementedError( "compute_gdf_gradient_rhf_gamma: the analytic Ewald V_ne " "derivative currently supports 3D periodic systems only." ) _validate_supported_gradient_result( result, caller="compute_gdf_gradient_rhf_gamma", ) if not np.isclose(float(alpha_hf), 1.0, rtol=0.0, atol=0.0): raise ValueError( "compute_gdf_gradient_rhf_gamma: alpha_hf=1.0 is required for " "the supported RHF result" ) retained_madelung = _validated_result_madelung( system, result, caller="compute_gdf_gradient_rhf_gamma", ) if madelung is None: madelung = retained_madelung elif not np.isclose( float(madelung), retained_madelung, rtol=0.0, atol=1e-15 ): raise ValueError( "compute_gdf_gradient_rhf_gamma: madelung must match the " "converged SCF result" ) missing = object() result_rcut_strategy = getattr(result, "gdf_rcut_strategy", missing) if result_rcut_strategy is missing: raise ValueError( "compute_gdf_gradient_rhf_gamma: missing compcell rcut provenance" ) result_linear_dep_thr = float( getattr(result, "gdf_linear_dep_threshold", float("nan")) ) result_rcut_precision = float( getattr(result, "gdf_rcut_precision", float("nan")) ) result_fit_cutoff_2c = float( getattr(result, "gdf_fit_cutoff_2c", float("nan")) ) result_fit_cutoff_3c = float( getattr(result, "gdf_fit_cutoff_3c", float("nan")) ) result_aux_fingerprint = str( getattr(result, "aux_basis_fingerprint", "") ) result_compcell_eta = float( getattr(result, "compcell_eta", float("nan")) ) result_v_ne_ke_cutoff = float( getattr(result, "v_ne_ke_cutoff", float("nan")) ) result_lattice_cutoff = float( getattr(result, "gdf_lattice_cutoff_bohr", float("nan")) ) result_nuclear_cutoff = float( getattr(result, "gdf_nuclear_cutoff_bohr", float("nan")) ) if not ( np.isfinite(result_linear_dep_thr) and result_linear_dep_thr >= 0.0 and np.isfinite(result_rcut_precision) and result_rcut_precision > 0.0 and np.isfinite(result_fit_cutoff_2c) and result_fit_cutoff_2c > 0.0 and np.isfinite(result_fit_cutoff_3c) and result_fit_cutoff_3c > 0.0 and np.isfinite(result_lattice_cutoff) and result_lattice_cutoff > 0.0 and np.isfinite(result_nuclear_cutoff) and result_nuclear_cutoff > 0.0 and bool(result_aux_fingerprint) and np.isfinite(result_compcell_eta) and result_compcell_eta > 0.0 and np.isfinite(result_v_ne_ke_cutoff) and result_v_ne_ke_cutoff > 0.0 ): raise ValueError( "compute_gdf_gradient_rhf_gamma: invalid compcell fit provenance" ) if _basis_fingerprint(aux_basis) != result_aux_fingerprint: raise ValueError( "compute_gdf_gradient_rhf_gamma: auxiliary basis content does " "not match the converged SCF result" ) if compcell_eta is None: compcell_eta = result_compcell_eta elif not np.isclose( float(compcell_eta), result_compcell_eta, rtol=0.0, atol=0.0 ): raise ValueError( "compute_gdf_gradient_rhf_gamma: compcell_eta must match the " "converged SCF result" ) if v_ne_ke_cutoff is None: v_ne_ke_cutoff = result_v_ne_ke_cutoff elif not np.isclose( float(v_ne_ke_cutoff), result_v_ne_ke_cutoff, rtol=0.0, atol=0.0 ): raise ValueError( "compute_gdf_gradient_rhf_gamma: v_ne_ke_cutoff must match the " "converged SCF result" ) if lattice_opts is None: lattice_opts = LatticeSumOptions() lattice_opts.cutoff_bohr = result_lattice_cutoff lattice_opts.nuclear_cutoff_bohr = result_nuclear_cutoff elif not ( np.isclose( float(lattice_opts.cutoff_bohr), result_lattice_cutoff, rtol=0.0, atol=1e-14, ) and np.isclose( float(lattice_opts.nuclear_cutoff_bohr), result_nuclear_cutoff, rtol=0.0, atol=1e-14, ) ): raise ValueError( "compute_gdf_gradient_rhf_gamma: lattice_opts cutoffs must " "match the converged SCF result" ) if gauge_lat_opts is None: from .pbc_gdf import _gauge_lat_opts_ewald_3d gauge_lat_opts = _gauge_lat_opts_ewald_3d(lattice_opts, system) else: if gauge_lat_opts.coulomb_method != CoulombMethod.EWALD_3D: raise ValueError( "compute_gdf_gradient_rhf_gamma: gauge_lat_opts must use " "EWALD_3D" ) if not ( np.isclose( float(gauge_lat_opts.cutoff_bohr), result_lattice_cutoff, rtol=0.0, atol=1e-14, ) and np.isclose( float(gauge_lat_opts.nuclear_cutoff_bohr), result_nuclear_cutoff, rtol=0.0, atol=1e-14, ) ): raise ValueError( "compute_gdf_gradient_rhf_gamma: gauge_lat_opts cutoffs " "must match the converged SCF result" ) n_elec = system.n_electrons() if n_elec % 2 != 0: raise ValueError( "compute_gdf_gradient_rhf_gamma: closed-shell only " f"(got {n_elec} electrons)" ) n_occ = n_elec // 2 D = np.asarray(result.density, dtype=np.float64) C_occ = np.asarray(result.mo_coeffs[:, :n_occ], dtype=np.float64) n_atoms = len(system.unit_cell) n_orb = basis.nbasis # Build or reuse the geometry-dependent (density-independent) cache. if cache is None: cache = _build_compcell_gradient_cache( system, basis, aux_basis, compcell_eta=float(compcell_eta), lattice_opts=lattice_opts, linear_dep_thr=result_linear_dep_thr, rcut_strategy=result_rcut_strategy, rcut_precision=result_rcut_precision, ) cache_provenance_matches = ( str(getattr(cache, "system_fingerprint", "")) == _system_fingerprint(system) and str(getattr(cache, "ao_basis_fingerprint", "")) == _basis_fingerprint(basis) and str(getattr(cache, "aux_basis_fingerprint", "")) == result_aux_fingerprint and np.isclose( float(getattr(cache, "compcell_eta", float("nan"))), result_compcell_eta, rtol=0.0, atol=0.0, ) and np.isclose( float(getattr(cache, "lattice_cutoff_bohr", float("nan"))), result_lattice_cutoff, rtol=0.0, atol=0.0, ) and np.isclose( float(getattr(cache, "nuclear_cutoff_bohr", float("nan"))), result_nuclear_cutoff, rtol=0.0, atol=0.0, ) and getattr(cache, "rcut_strategy", missing) == _rcut_strategy_value(result_rcut_strategy) and np.isclose( float(getattr(cache, "rcut_precision", float("nan"))), result_rcut_precision, rtol=0.0, atol=0.0, ) ) if ( not cache_provenance_matches or cache.n_fit != int(getattr(result, "n_fit", -1)) or cache.linear_dep_thr != result_linear_dep_thr or not np.isclose( float(cache.lat_opts_2c.cutoff_bohr), result_fit_cutoff_2c, rtol=0.0, atol=1e-12, ) or not np.isclose( float(cache.lat_opts_3c.cutoff_bohr), result_fit_cutoff_3c, rtol=0.0, atol=1e-12, ) ): raise ValueError( "compute_gdf_gradient_rhf_gamma: compcell cache provenance " "does not match the converged SCF fit" ) # ---- 1-electron Pulay + overlap + nuclear terms --------------------- # These are the same as the existing periodic gradient (G1a) for the # Ewald path. We build the lattice-resolved density and W matrices # using the GDF SCF's hcore (not rebuilt molecule Fock, since at # Γ the GDF Fock doesn't have the G=0 gauge issue of EWALD_3D). # Build the overlap lattice for the cell list. S_lat = compute_overlap_lattice(basis, system, lattice_opts) # Lattice-resolved density: at Γ, D(g) = D_Γ for all images. D_gamma = np.asarray(D, dtype=np.float64) D_set = compute_overlap_lattice(basis, system, lattice_opts) # Homogeneous Γ density: D(g) = D_Γ for all image cells. # At Γ, k=0 Bloch phase is 1 in every cell, so the real-space # density is identical in all images. The SCF energy is built from # this homogeneous density; the gradient must use it too. home_cell_only = all( tuple(int(v) for v in np.asarray(c.index).reshape(3)) == (0, 0, 0) for c in D_set.cells ) zero_block = np.zeros_like(D_gamma) for c_idx in range(len(D_set.cells)): idx = tuple(int(v) for v in np.asarray(D_set.cells[c_idx].index).reshape(3)) D_set.set_block( c_idx, D_gamma if (idx == (0, 0, 0) or not home_cell_only) else zero_block ) # Energy-weighted density W for the overlap-Lagrangian. # Use the converged Fock eigenvalues directly -- the GDF route doesn't # have the EWALD_3D G=0 gauge shift, so e from the SCF is fine. eps = np.asarray(result.mo_energies, dtype=np.float64) C = np.asarray(result.mo_coeffs, dtype=np.float64) W_gamma = 2.0 * (C[:, :n_occ] * eps[:n_occ][None, :]) @ C[:, :n_occ].T W_set = compute_overlap_lattice(basis, system, lattice_opts) for c_idx in range(len(W_set.cells)): idx = tuple(int(v) for v in np.asarray(W_set.cells[c_idx].index).reshape(3)) W_set.set_block( c_idx, W_gamma if (idx == (0, 0, 0) or not home_cell_only) else zero_block, ) grad = np.zeros((n_atoms, 3), dtype=np.float64) # Nuclear repulsion. grad += np.asarray(nuclear_repulsion_gradient_per_cell(system, gauge_lat_opts)) # Overlap Lagrangian: -tr(W dS/dR). grad += np.asarray( overlap_lattice_gradient_contribution(basis, system, W_set, lattice_opts) ) # Kinetic + nuclear-attraction Pulay: tr(D d(T+V)/dR). grad += np.asarray( kinetic_lattice_gradient_contribution(basis, system, D_set, lattice_opts) ) # V_ne Hellmann-Feynman and AO-centre response in the exact # FT-based Ewald gauge used by the SCF. from .periodic_v_ne_gradient import ( compute_v_ne_ewald_3d_ft_gamma_gradient, ) grad += compute_v_ne_ewald_3d_ft_gamma_gradient( basis, system, gauge_lat_opts, D, ke_cutoff=float(v_ne_ke_cutoff), ) # ---- 2-electron DF gradient ----------------------------------------- grad_J = _compute_j_gradient_compcell( system, basis, D, cache, ) if alpha_hf > 0.0: grad_K = _compute_k_gradient_compcell( system, basis, D, C_occ, float(alpha_hf), cache, ) # E_2e = E_J - 1/4 a_HF tr(D.K_shifted) # = E_J - 1/4 a_HF tr(D.K_raw) + 1/4 a_HF ξ tr(D.S.D.S) # dE_2e/dR = dE_J/dR - 1/4 a_HF dE_K_raw/dR + 1/2 a_HF ξ tr(D.S.D.dS/dR) grad += grad_J + grad_K # Exxdiv shift gradient: K_shifted = K_raw + ξ.S.D.S (ADDED in # apply_exxdiv_ewald_to_K). So E_K = -1/4 Tr(D.K_raw) - 1/4 ξ Tr(D.S.D.S). # d(-1/4 ξ Tr(D.S.D.S))/dR at fixed D = -1/2 ξ Tr(D.S.D.dS/dR). # overlap_lattice_gradient_contribution computes -Tr(W.dS/dR), # so pass W_exx = +1/2 ξ D.S.D. if abs(madelung) > 0.0: S = np.asarray(result.overlap, dtype=np.float64) DSD = D @ S @ D W_exx_gamma = 0.5 * madelung * DSD W_exx_set = compute_overlap_lattice(basis, system, lattice_opts) for c_idx in range(len(W_exx_set.cells)): idx = tuple( int(v) for v in np.asarray(W_exx_set.cells[c_idx].index).reshape(3) ) W_exx_set.set_block( c_idx, W_exx_gamma if (idx == (0, 0, 0) or not home_cell_only) else np.zeros_like(W_exx_gamma), ) grad += np.asarray( overlap_lattice_gradient_contribution( basis, system, W_exx_set, lattice_opts ) ) else: grad += grad_J return grad, cache def compute_gdf_gradient( system: PeriodicSystem, basis: BasisSet, result, # PBCGDFResult *, aux_basis: Optional[BasisSet] = None, aux_basis_name: Optional[str] = None, compcell_eta: Optional[float] = None, lattice_opts: Optional[LatticeSumOptions] = None, v_ne_ke_cutoff: Optional[float] = None, _fit_state: Optional[_CompcellFitState] = None, ) -> np.ndarray: """Compute the GDF analytic gradient from a converged PBCGDFResult. Convenience wrapper around :func:`compute_gdf_gradient_rhf_gamma`. Handles gauge setup, gradient cache construction, and validation of the exchange-divergence constant retained by the converged SCF result. Parameters ---------- system, basis Periodic system and AO basis. result Converged :class:`PBCGDFResult` from :func:`vibeqc.run_pbc_gdf_rhf`. aux_basis Auxiliary basis. If None, built from ``aux_basis_name``. aux_basis_name Aux basis name for auto-construction. ``None`` uses the name retained on ``result``; an explicit name must match the converged SCF value. compcell_eta Compensation exponent. ``None`` uses the value retained on ``result``; an explicit value must match the converged SCF value. lattice_opts Lattice-sum options. If None, reconstructs the base and nuclear cutoffs retained by the converged SCF result. v_ne_ke_cutoff Ewald ``V_ne`` reciprocal cutoff. ``None`` uses the cutoff retained on ``result``; an explicit value must match that converged SCF value. Returns ------- grad : (n_atoms, 3) ndarray in Hartree/bohr. """ if int(system.dim) != 3: raise NotImplementedError( "compute_gdf_gradient: the analytic Ewald V_ne derivative " "currently supports 3D periodic systems only." ) _validate_supported_gradient_result(result, caller="compute_gdf_gradient") result_v_ne_ke_cutoff = float( getattr(result, "v_ne_ke_cutoff", float("nan")) ) if not np.isfinite(result_v_ne_ke_cutoff): raise ValueError( "compute_gdf_gradient: the converged result does not retain " "its Ewald V_ne reciprocal cutoff." ) if v_ne_ke_cutoff is None: v_ne_ke_cutoff = result_v_ne_ke_cutoff elif float(v_ne_ke_cutoff) != result_v_ne_ke_cutoff: raise ValueError( "compute_gdf_gradient: v_ne_ke_cutoff must match the converged " f"SCF value ({result_v_ne_ke_cutoff:g} Ha)." ) result_aux_basis_name = str(getattr(result, "aux_basis_name", "")) if not result_aux_basis_name: raise ValueError( "compute_gdf_gradient: the converged result does not retain " "its auxiliary-basis name." ) if aux_basis_name is None: aux_basis_name = result_aux_basis_name elif str(aux_basis_name) != result_aux_basis_name: raise ValueError( "compute_gdf_gradient: aux_basis_name must match the converged " f"SCF value ({result_aux_basis_name!r})." ) if aux_basis is not None: supplied_aux_basis_name = str(getattr(aux_basis, "name", "")) if ( supplied_aux_basis_name and supplied_aux_basis_name != result_aux_basis_name ): raise ValueError( "compute_gdf_gradient: aux_basis must match the converged " f"SCF value ({result_aux_basis_name!r})." ) result_aux_fingerprint = str( getattr(result, "aux_basis_fingerprint", "") ) if not result_aux_fingerprint: raise ValueError( "compute_gdf_gradient: the converged result does not retain " "its auxiliary-basis fingerprint." ) result_compcell_eta = float( getattr(result, "compcell_eta", float("nan")) ) if not np.isfinite(result_compcell_eta): raise ValueError( "compute_gdf_gradient: the converged result does not retain " "its compcell_eta value." ) if compcell_eta is None: compcell_eta = result_compcell_eta elif float(compcell_eta) != result_compcell_eta: raise ValueError( "compute_gdf_gradient: compcell_eta must match the converged " f"SCF value ({result_compcell_eta:g})." ) missing = object() result_rcut_strategy = getattr(result, "gdf_rcut_strategy", missing) if result_rcut_strategy is missing: raise ValueError( "compute_gdf_gradient: the converged result does not retain " "its compcell rcut strategy." ) if result_rcut_strategy is not None: result_rcut_strategy = str(result_rcut_strategy) result_rcut_precision = float( getattr(result, "gdf_rcut_precision", float("nan")) ) result_linear_dep_thr = float( getattr(result, "gdf_linear_dep_threshold", float("nan")) ) if not np.isfinite(result_rcut_precision) or result_rcut_precision <= 0.0: raise ValueError( "compute_gdf_gradient: the converged result does not retain " "a valid compcell rcut precision." ) if not np.isfinite(result_linear_dep_thr) or result_linear_dep_thr < 0.0: raise ValueError( "compute_gdf_gradient: the converged result does not retain " "a valid GDF linear-dependence threshold." ) result_lattice_cutoff = float( getattr(result, "gdf_lattice_cutoff_bohr", float("nan")) ) result_nuclear_cutoff = float( getattr(result, "gdf_nuclear_cutoff_bohr", float("nan")) ) result_fit_cutoff_2c = float( getattr(result, "gdf_fit_cutoff_2c", float("nan")) ) result_fit_cutoff_3c = float( getattr(result, "gdf_fit_cutoff_3c", float("nan")) ) retained_cutoffs = ( result_lattice_cutoff, result_nuclear_cutoff, result_fit_cutoff_2c, result_fit_cutoff_3c, ) if not all(np.isfinite(value) and value > 0.0 for value in retained_cutoffs): raise ValueError( "compute_gdf_gradient: the converged result does not retain " "valid base and resolved compcell lattice cutoffs." ) if lattice_opts is None: lattice_opts = LatticeSumOptions() lattice_opts.cutoff_bohr = result_lattice_cutoff lattice_opts.nuclear_cutoff_bohr = result_nuclear_cutoff elif not ( np.isclose( float(lattice_opts.cutoff_bohr), result_lattice_cutoff, rtol=0.0, atol=1e-14, ) and np.isclose( float(lattice_opts.nuclear_cutoff_bohr), result_nuclear_cutoff, rtol=0.0, atol=1e-14, ) ): raise ValueError( "compute_gdf_gradient: lattice_opts cutoffs must match the " "converged SCF values." ) result_madelung = _validated_result_madelung( system, result, caller="compute_gdf_gradient", ) from .pbc_gdf import _gauge_lat_opts_ewald_3d from .aux_basis import make_aux_basis_set gauge_opts = _gauge_lat_opts_ewald_3d(lattice_opts, system) mol = system.unit_cell_molecule() if aux_basis is None: aux_basis = make_aux_basis_set(mol, aux_name=aux_basis_name) if _basis_fingerprint(aux_basis) != result_aux_fingerprint: raise ValueError( "compute_gdf_gradient: auxiliary basis content must match the " "converged SCF result." ) # Build gradient cache (expensive, geometry-dependent, density-independent) cache = _build_compcell_gradient_cache( system, basis, aux_basis, compcell_eta=float(compcell_eta), lattice_opts=lattice_opts, linear_dep_thr=result_linear_dep_thr, rcut_strategy=result_rcut_strategy, rcut_precision=result_rcut_precision, fit_state=_fit_state, ) result_n_fit = int(getattr(result, "n_fit", -1)) if cache.n_fit != result_n_fit: raise ValueError( "compute_gdf_gradient: rebuilt compcell fit rank does not match " f"the converged SCF result ({cache.n_fit} != {result_n_fit})." ) if not np.isclose( cache.linear_dep_thr, result_linear_dep_thr, rtol=0.0, atol=0.0, ): raise ValueError( "compute_gdf_gradient: rebuilt compcell fit threshold does not " "match the converged SCF result." ) if not np.isclose( float(cache.lat_opts_2c.cutoff_bohr), result_fit_cutoff_2c, rtol=0.0, atol=1e-12, ) or not np.isclose( float(cache.lat_opts_3c.cutoff_bohr), result_fit_cutoff_3c, rtol=0.0, atol=1e-12, ): raise ValueError( "compute_gdf_gradient: rebuilt compcell fit cutoffs do not match " "the converged SCF result." ) grad, _ = compute_gdf_gradient_rhf_gamma( system, basis, result, aux_basis=aux_basis, compcell_eta=float(compcell_eta), alpha_hf=1.0, lattice_opts=lattice_opts, cache=cache, madelung=result_madelung, gauge_lat_opts=gauge_opts, v_ne_ke_cutoff=v_ne_ke_cutoff, ) return grad