diff --git a/.gitignore b/.gitignore index 36b13f1..505a3b1 100644 --- a/.gitignore +++ b/.gitignore @@ -1,176 +1,10 @@ -# ---> Python -# Byte-compiled / optimized / DLL files +# Python-generated files __pycache__/ -*.py[cod] -*$py.class - -# C extensions -*.so - -# Distribution / packaging -.Python +*.py[oc] build/ -develop-eggs/ dist/ -downloads/ -eggs/ -.eggs/ -lib/ -lib64/ -parts/ -sdist/ -var/ wheels/ -share/python-wheels/ -*.egg-info/ -.installed.cfg -*.egg -MANIFEST +*.egg-info -# PyInstaller -# Usually these files are written by a python script from a template -# before PyInstaller builds the exe, so as to inject date/other infos into it. -*.manifest -*.spec - -# Installer logs -pip-log.txt -pip-delete-this-directory.txt - -# Unit test / coverage reports -htmlcov/ -.tox/ -.nox/ -.coverage -.coverage.* -.cache -nosetests.xml -coverage.xml -*.cover -*.py,cover -.hypothesis/ -.pytest_cache/ -cover/ - -# Translations -*.mo -*.pot - -# Django stuff: -*.log -local_settings.py -db.sqlite3 -db.sqlite3-journal - -# Flask stuff: -instance/ -.webassets-cache - -# Scrapy stuff: -.scrapy - -# Sphinx documentation -docs/_build/ - -# PyBuilder -.pybuilder/ -target/ - -# Jupyter Notebook -.ipynb_checkpoints - -# IPython -profile_default/ -ipython_config.py - -# pyenv -# For a library or package, you might want to ignore these files since the code is -# intended to run in multiple environments; otherwise, check them in: -# .python-version - -# pipenv -# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. -# However, in case of collaboration, if having platform-specific dependencies or dependencies -# having no cross-platform support, pipenv may install dependencies that don't work, or not -# install all needed dependencies. -#Pipfile.lock - -# UV -# Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control. -# This is especially recommended for binary packages to ensure reproducibility, and is more -# commonly ignored for libraries. -#uv.lock - -# poetry -# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. -# This is especially recommended for binary packages to ensure reproducibility, and is more -# commonly ignored for libraries. -# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control -#poetry.lock - -# pdm -# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. -#pdm.lock -# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it -# in version control. -# https://pdm.fming.dev/latest/usage/project/#working-with-version-control -.pdm.toml -.pdm-python -.pdm-build/ - -# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm -__pypackages__/ - -# Celery stuff -celerybeat-schedule -celerybeat.pid - -# SageMath parsed files -*.sage.py - -# Environments -.env +# Virtual environments .venv -env/ -venv/ -ENV/ -env.bak/ -venv.bak/ - -# Spyder project settings -.spyderproject -.spyproject - -# Rope project settings -.ropeproject - -# mkdocs documentation -/site - -# mypy -.mypy_cache/ -.dmypy.json -dmypy.json - -# Pyre type checker -.pyre/ - -# pytype static type analyzer -.pytype/ - -# Cython debug symbols -cython_debug/ - -# PyCharm -# JetBrains specific template is maintained in a separate JetBrains.gitignore that can -# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore -# and can be added to the global gitignore or merged into this file. For a more nuclear -# option (not recommended) you can uncomment the following to ignore the entire idea folder. -#.idea/ - -# Ruff stuff: -.ruff_cache/ - -# PyPI configuration file -.pypirc - diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..9e4cc63 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "libxc"] + path = libxc + url = https://gitlab.com/libxc/libxc diff --git a/.python-version b/.python-version new file mode 100644 index 0000000..24ee5b1 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.13 diff --git a/README.md b/README.md index 17c553e..14a6ce0 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,69 @@ # funxc -drop-in replacement to libxc functionals in jax \ No newline at end of file +Differentiable exchange–correlation functionals in JAX: a clean-room +reimplementation of [libxc](https://libxc.gitlab.io/), aiming for numeric +compatibility with `libxc` 7.0.0 while being composable, autodifferentiable, +and hardware-agnostic. + +## Why + +1. `libxc` defines its mathematics in Maple (non-FOSS) and ships generated C; + `funxc` defines each functional *family* once as a pure JAX function. +2. JAX provides a hardware-agnostic acceleration layer (CPU/GPU/TPU) and + exact derivatives to any order via autodiff, rather than generated derivative code. +3. Differentiable functionals open the path from SCF-on-a-BO-surface toward + end-to-end gradient methods. + +Most of libxc's 600+ functionals are a small set of analytic kernels with +different parameters. `funxc` implements the kernels as code and the +parameters as a TOML database: adding a parameter-only variant is one TOML +entry. + +## Usage + +JAX-native: + +```python +import funxc +import jax, jax.numpy as jnp + +f = funxc.functional("GGA_X_PBE") # polarized by default +zk = f.exc(rho, sigma) # rho (N,2), sigma (N,3) +out = f.exc_vxc(rho, sigma) # zk, vrho, vsigma via autodiff + +# per-point functions compose with any JAX transform +grad_e = jax.grad(f.energy_density, argnums=(0, 1)) +``` + +pylibxc drop-in shim: + +```python +from funxc.libxc_compat import LibXCFunctional +f = LibXCFunctional("gga_x_pbe", "polarized") +out = f.compute({"rho": rho, "sigma": sigma}) # zk, vrho, vsigma arrays +``` + +## Compatibility & testing + +Every registered functional is validated against libxc 7.0.0's own +regression data (energies and first derivatives, polarized and unpolarized) +to a tolerance of `2e-9 + 1e-7·|ref|` for a given `|ref|` per functional; +typical agreement is machine precision. Kernels are guarded so that values +and gradients (including second derivatives) stay finite at ρ→0 tails, σ=0, +ζ=±1, s→∞ and τ→0. + +``` +uv sync +uv run pytest +``` + +## Status + +Stage 1 complete (39 functionals): Slater exchange, PW92 + VWN local +correlation, the PBE rational-exchange family (PBE, PBEsol, revPBE, xPBE, +APBE, ...), RPBE, B88, the PBE H(t) correlation family, LYP, and the B97 +power-series family (B97-D, B97-3c, B97-GGA1, the HCTH functionals, HLE16). +Meta-GGA plumbing (τ and ∇²ρ inputs, vtau/vlapl via autodiff) is in place +and validated with the LTA family. Stage 2 kernels (TPSS, SCAN/r2SCAN) and +hybrid/range-separated recipes are next; see `mapOfFunctionals.md` for the +roadmap. diff --git a/libxc b/libxc new file mode 160000 index 0000000..7bd5bb4 --- /dev/null +++ b/libxc @@ -0,0 +1 @@ +Subproject commit 7bd5bb41415968db94c499a2f093309c9a2dcf53 diff --git a/mapOfFunctionals.md b/mapOfFunctionals.md new file mode 100644 index 0000000..30e2b02 --- /dev/null +++ b/mapOfFunctionals.md @@ -0,0 +1,292 @@ +# Libxc Tree of Knowledge: A Categorical Index of Exchange–Correlation Functionals by Shared Analytic Kernel + +## TL;DR +- Libxc's 600+ functional entries collapse to roughly **50–70 distinct analytic enhancement-factor kernels**; the great majority of "different" functionals are the *same* math kernel with a different parameter tuple (e.g. ~30+ PBE-type exchange functionals share one rational `Fx(s)` form), so a clean-room reimplementation needs only the small kernel set as code plus a parameter database. +- The tree has four rungs (LDA → GGA → meta-GGA, plus kinetic-energy **K** functionals), each split into **X / C / XC**, with hybrids and range-separation treated as a thin recipe layer of mixing coefficients (α, β, ω) over the semilocal kernels — exact exchange is host-supplied, never a kernel to implement. +- The irreducible "long tail" is small: genuinely unique forms include Becke–Roussel (BR89) exchange-hole inversion, LYP correlation, the VWN/PW92 Padé-and-log local-correlation parametrizations, OPTX, AM05, the B97 power series, and the SCAN/TPSS iso-orbital meta-GGA interpolations. + +## Key Findings + +1. **The compression ratio is dramatic.** MolSSI's official page states: *"The library currently implements over 600 density functional approximations, and is used by over 30 electronic structure programs ranging from atomic structure calculations to large quantum chemistry and materials modeling programs."* The original 2012 implementation (Marques, Oliveira & Burnus, *Comput. Phys. Commun.* 183, 2272) contained "around 180" functionals, while Marques' own estimate of the total number of *known* approximations at that time was "of the order of 250–300" (quoted in Lehtola's GPAW-2021 talk). Today's 600+ entries reduce to ~50–70 distinct analytic kernels because parameter-only variants dominate. The single largest family — PBE-type GGA exchange — alone subsumes dozens of functionals through the (κ, μ) tuple. + +2. **The PBE rational enhancement `Fx(s) = 1 + κ − κ/(1 + μ·s²/κ)` is the workhorse kernel.** It originated as Becke's B86 form and is shared verbatim (in libxc's `maple/gga_x_pbe.mpl`, where it reads `pbe_f0 := s -> 1 + kappa*(1 - kappa/(kappa + mu*s^2))`) by PBE, PBEsol, revPBE, xPBE, APBE, PBEmol, PBEfe, PBEint, PBE_TCA, and the lambda family. The TPSS and SCAN meta-GGAs reuse the same rational shell with a more complex argument `x` replacing `μ·s²`. + +3. **Local correlation is dominated by two parametrization forms:** the VWN Padé-in-√rs log/arctan form and the Perdew–Wang 1992 `G(rs) = −2A(1+α₁·rs)·ln[1 + 1/(2A(β₁·rs^½+…))]` form. PW92 (`lda_c_pw`) is the LDA base of nearly all GGA/mGGA correlation (PBE, PW91, TPSS, SCAN). + +4. **The B97 power series `u(x) = γx²/(1+γx²)` with `Σ cᵢ·uⁱ` is a second universal kernel**, spanning B97, B97-1, B97-2, HCTH (all variants), ωB97, ωB97X, B97-D and the meta-GGA B97M extensions — they differ only in the cᵢ coefficient vectors and the γ values. + +5. **Meta-GGA exchange splits into three kernel lineages:** (a) TPSS/revTPSS/SCAN/rSCAN/r2SCAN — the PBE rational shell with an iso-orbital argument built from `z = τW/τ` or `α = (τ−τW)/τunif`; (b) Minnesota (M06-L, M05, M06, M08, M11, MN15) — PBE/PW86 exchange times a kinetic-energy power series `f(w)` plus a VS98-style term; (c) Becke–Roussel BR89 exchange-hole inversion (used by BR89, B00, and underlying the MVS, TM, and mBJ potentials). + +--- + +## Details — The Tree + +### RUNG 0 — Variable dependencies (top-level split) +- **LDA**: ε depends on ρ (and ζ = spin polarization) only. +- **GGA**: adds σ = |∇ρ|². Reduced gradient `s = |∇ρ|/(2(3π²)^(1/3)ρ^(4/3))`, or `x = |∇ρ|/ρ^(4/3)`; note `s = x/(2(3π²)^(1/3))`. +- **meta-GGA**: adds τ (kinetic energy density) and/or ∇²ρ (Laplacian). Key dimensionless ingredients: `z = τW/τ`, `α = (τ−τW)/τunif`, where `τW = |∇ρ|²/(8ρ)`, `τunif = (3/10)(3π²)^(2/3)ρ^(5/3)`. +- Each rung × {X exchange, C correlation, XC combined, K kinetic}. + +--- + +### LDA RUNG + +#### Kernel L-X1: Dirac–Slater exchange +- **Variables:** ρ. **Reference:** Dirac/Slater (`LDA_X`). +- **Form:** εx = −(3/4)(3/π)^(1/3) ρ^(1/3) per particle; spin via εx[ρ↑,ρ↓] = ½(εx[2ρ↑]+εx[2ρ↓]). +- **Members / parameters:** `LDA_X` (canonical Cx), `LDA_X_REL` (relativistic correction factor), `LDA_X_2D`, `LDA_X_1D_SOFT`/`_EXPONENTIAL`, `XALPHA` (Slater Xα, adjustable α default 0.7), `LDA_X_RAE`, `LDA_X_SLOC` (a=1.67, b=0.3), `LDA_X_YUKAWA`/`_ERF` (screened). Variation = the prefactor Cx and the dimension. +- **Citation:** Dirac, *Proc. Camb. Phil. Soc.* 26, 376 (1930); Slater, *Phys. Rev.* 81, 385 (1951). + +#### Kernel L-C1: VWN Padé log/arctan local correlation +- **Variables:** ρ, ζ. **Reference:** Vosko–Wilk–Nusair 1980 (`LDA_C_VWN`). +- **Form:** `ec(rs,ζ) = ec^P + [ec^F − ec^P]·fc(ζ)ζ⁴ + (fc(ζ)/fc''(0))(1−ζ⁴)αc(rs)`, each of ec^P, ec^F, αc parametrized by + `G(y) = A{ln(y²/X(y)) + (2b/Q)arctan(Q/(2y+b)) − (b·y0/X(y0))[ln((y−y0)²/X(y)) + (2(b+2y0)/Q)arctan(Q/(2y+b))]}`, + with `y = √rs`, `X(y) = y²+by+c`, `Q = √(4c−b²)`, `fc''(0) = 4/[9(2^(1/3)−1)]`. Paramagnetic channel constants: A = 0.0310907, b = 3.72744, c = 12.9352, x0 = −0.10498. +- **Members:** `LDA_C_VWN` (=VWN5), `LDA_C_VWN_RPA`, `LDA_C_VWN_1`, `_2`, `_3`, `_4` — variants differ in (A,b,c,x0) tuples and whether the RPA or Ceperley–Alder fit is used. VWN3 and VWN5 are the historically confused pair (VWN3 = RPA-fit body parametrization; VWN5 = Table-5 Monte-Carlo parametrization). +- **Citation:** Vosko, Wilk & Nusair, *Can. J. Phys.* 58, 1200 (1980), doi:10.1139/p80-159. + +#### Kernel L-C2: Perdew–Wang 1992 G(rs) local correlation +- **Variables:** ρ, ζ. **Reference:** PW92 (`LDA_C_PW`). +- **Form:** `G(rs) = −2A(1+α₁·rs)·ln[1 + 1/(2A(β₁·rs^(1/2) + β₂·rs + β₃·rs^(3/2) + β₄·rs^(p+1)))]`, spin-interpolated with the von Barth–Hedin fc(ζ). +- **Parameter table (three channels: ec0 paramagnetic, ec1 ferromagnetic, αc spin-stiffness):** + + | param | ec0 | ec1 | αc | + |---|---|---|---| + | A | 0.031091 | 0.015545 | 0.016887 | + | α₁ | 0.21370 | 0.20548 | 0.11125 | + | β₁ | 7.5957 | 14.1189 | 10.357 | + | β₂ | 3.5876 | 6.1977 | 3.6231 | + | β₃ | 1.6382 | 3.3662 | 0.88026 | + | β₄ | 0.49294 | 0.62517 | 0.49671 | + + p = 1 for `LDA_C_PW`; p = 3/4 for `LDA_C_PW_RPA`; `LDA_C_PW_MOD` uses higher-precision A (0.0310907…). +- **Members:** `LDA_C_PW`, `LDA_C_PW_RPA`, `LDA_C_PW_MOD`, `LDA_C_OB_PW` (Ortiz–Ballone refit). +- **Citation:** Perdew & Wang, *Phys. Rev. B* 45, 13244 (1992). + +#### Kernel L-C3: PZ81 piecewise (high/low-density) correlation +- **Form:** for rs ≥ 1, `ec = γ/(1+β₁√rs+β₂·rs)`; for rs < 1, `ec = A·ln rs + B + C·rs·ln rs + D·rs`. **Members:** `LDA_C_PZ`, `LDA_C_PZ_MOD`, `LDA_C_OB_PZ`. **Citation:** Perdew & Zunger, *Phys. Rev. B* 23, 5048 (1981). **Stability note:** second-derivative discontinuity at rs = 1. + +#### Kernel L-C4: Chachiyo single-formula correlation +- **Form:** `ec = a·ln(1 + b/rs + b/rs²)`-type closed two-parameter expression. **Members:** `LDA_C_CHACHIYO`, `LDA_C_CHACHIYO_MOD`, `GGA_C_CHACHIYO` (gradient extension). **Citation:** Chachiyo, *J. Chem. Phys.* 145, 021101 (2016). + +#### Kernel L-C5: Wigner / small-rational one-offs +- `LDA_C_WIGNER` (ec = −a/(rs+b)), `LDA_C_RPA`, `LDA_C_HL` (Hedin–Lundqvist), `LDA_C_GL` (Gunnarsson–Lundqvist), `LDA_C_vBH`, `LDA_C_RC04`, `LDA_C_W20`, `LDA_C_GOMBAS` — each a distinct small rational/log form; mostly irreducible. +- **Warm-dense-matter / finite-T XC:** `LDA_XC_KSDT`, `LDA_XC_GDSMFB`, `LDA_XC_CORRKSDT` — temperature-dependent Padé parametrizations (shared form, different fits). +- **Electron–proton (NEO):** `LDA_C_EPC17`, `_17_2`, `_18_1`, `_18_2` — implemented via the spin-channel trick (electron density as "up", proton density as "down"). + +#### Kernel L-K1: Thomas–Fermi kinetic +- **Form:** `τ = CF·ρ^(5/3)`, `CF = (3/10)(3π²)^(2/3)`. **Members:** `LDA_K_TF`, `LDA_K_LP` (Lee–Parr Gaussian), `LDA_K_ZLP`. **Citation:** Thomas (1927), Fermi (1928). + +--- + +### GGA RUNG — EXCHANGE + +#### Kernel G-X1: PBE-type rational enhancement (THE big family) +- **Variables:** ρ, σ. **Reference:** PBE (`GGA_X_PBE`); form first from Becke B86. +- **Form (libxc `maple/gga_x_pbe.mpl`):** `Fx(s) = 1 + κ − κ/(1 + μ·s²/κ)`, εx = εx^LDA·Fx(s). Gradient expansion `Fx = 1 + μ·s² + O(s⁴)`. +- **Member / parameter table:** + + | Functional | κ | μ | + |---|---|---| + | `GGA_X_PBE` | 0.804 | 0.2195149727645171 | + | `GGA_X_PBE_SOL` (PBEsol) | 0.804 | 10/81 = 0.12345679 | + | `GGA_X_PBE_R` (revPBE) | 1.245 | 0.2195149727645171 | + | `GGA_X_XPBE` (Xu–Goddard) | 0.91954 | 0.23889 | + | `GGA_X_APBE` | 0.804 | 0.260 | + | `GGA_X_PBE_MOL` | 0.804 | ≈0.27583 | + | `GGA_X_PBEFE` | 0.437 | 0.346 | + | `GGA_X_PBE_TCA` | 1.227 | 0.2195149727645171 | + | `GGA_X_PBEINT` | 0.804 | μ(s) interpolating 10/81→0.21951, α=0.197 | + | `GGA_X_LAMBDA_{LO,CH,OC2}_N` | N-dependent κ (tightened Lieb–Oxford) | 0.2195149727645171 | + + Also share this rational shell (refit params): `GGA_X_PBE_JSJR`, `PBE_GAUSSIAN`, `mPBE` (Adamo–Barone), and vdW reparametrizations (`PBEK1_VDW`, `OPTPBE_VDW`). The VASP wiki documents the PBE→PBEsol move concisely: PBE's "μ=0.21951 in exchange and β=0.066725 in correlation are changed to μ=10/81≈0.12345679 and β=0.046 to get the PBEsol functional." +- **Citation:** Perdew, Burke & Ernzerhof, *PRL* 77, 3865 (1996); Becke, *J. Chem. Phys.* 84, 4524 (1986). + +#### Kernel G-X2: RPBE exponential enhancement +- **Form:** `Fx(s) = 1 + κ(1 − exp(−μ·s²/κ))`. Same (κ, μ) = (0.804, 0.21951) as PBE but a distinct functional form (exponential, not rational), hence a separate kernel. +- **Members:** `GGA_X_RPBE`; related `GGA_X_PBEpow`. **Citation:** Hammer, Hansen & Nørskov, *Phys. Rev. B* 59, 7413 (1999). + +#### Kernel G-X3: B88 arcsinh exchange +- **Variables:** ρ, σ. **Form:** `Fx = 1 + (β/Cx)·x²/(1 + γ·β·x·arcsinh(x))`, `x = |∇ρ|/ρ^(4/3)`, γ = 6. +- **Member / parameters:** `GGA_X_B88` (β = 0.0042), `GGA_X_B88_6311G` (β = 0.0051), `GGA_X_OPTB88_VDW`, `GGA_X_MB88`, `GGA_X_EB88` (excogitated β≈0.0050), `GGA_X_B88M`, LRC variants (β = 0.323). All differ only in β (and γ for some). +- **Citation:** Becke, *Phys. Rev. A* 38, 3098 (1988). **Stability note:** arcsinh is well-behaved; the asymptotic exchange potential is unbounded (the −1/(2r) energy-density decay is designed-in). + +#### Kernel G-X4: B86 / B86b family +- **Form:** `Fx = 1 + (μ s²)/(1 + γ s²)^p` with p ∈ {1, 4/5}. **Members:** `GGA_X_B86`, `GGA_X_B86_MGC`, `GGA_X_B86_R`, `GGA_X_OPTB86B_VDW`. (A close cousin of the PBE rational shell with the added power p.) +- **Citation:** Becke, *J. Chem. Phys.* 84, 4524 (1986); 85, 7184 (1986). + +#### Kernel G-X5: PW91 exchange (sinh-based rational) +- **Form:** `Fx = (1 + a₁·s·arcsinh(a₂·s) + (a₃ + a₄·exp(−b·s²))s²)/(1 + a₁·s·arcsinh(a₂·s) + a₅·s⁴)`. **Members:** `GGA_X_PW91`, `GGA_X_MPW91` (Adamo–Barone), `GGA_X_PW91_MOD`. **Citation:** Perdew et al., *Phys. Rev. B* 46, 6671 (1992). + +#### Kernel G-X6: OPTX (sum-of-rational-squared) +- **Form:** `Fx = a₁ + a₂·(γx²/(1+γx²))²`, `x = |∇ρ|/ρ^(4/3)`. Parameters: a₁ = 1.05151, a₂ = 1.43169, γ = 0.006. +- **Members:** `GGA_X_OPTX`, `GGA_X_ITYH_OPTX` (short-range erf recipe). **Citation:** Handy & Cohen, *Mol. Phys.* 99, 403 (2001). + +#### Kernel G-X7: G96 (Gill 1996) +- **Form:** Fx with an x^(3/2) gradient term, simpler than B88. **Member:** `GGA_X_G96`. **Citation:** Gill, *Mol. Phys.* 89, 433 (1996). + +#### Kernel G-X8: AM05 (Airy-gas / subsystem interpolation) +- **Form:** interpolation between LDA and Airy-gas (Laplacian-level) limits with a Langreth–Vosko-style index. **Members:** `GGA_X_AM05`, `GGA_X_AIRY`, `GGA_X_LAG`. **Citation:** Armiento & Mattsson, *Phys. Rev. B* 72, 085108 (2005). Largely irreducible. + +#### Kernel G-X9: Asymptotically-corrected potentials (LB94/LBα) +- **Form:** modeled on the *potential* (not energy): `vx = −β x²/(1+3β x·arcsinh(x))` added to the LDA potential. **Members:** `GGA_X_LB94`, `GGA_X_LBM`, `GGA_X_FD_LB94`. Potential-only — flag as a special case (no consistent εx). + +#### Kernel G-X10: HJS screened exchange-hole models +- **Form:** error-function-screened exchange-hole integral, parametrized per base functional. **Members:** `GGA_X_HJS_PBE`, `_PBE_SOL`, `_B88`, `_B97X`, `_B88_V2` — same HJS machinery, parameters per base. **Citation:** Henderson, Janesko & Scuseria, *J. Chem. Phys.* 128, 194105 (2008). Underlies HSE / range-separated builds. + +#### Misc one-off GGA exchange (irreducible long tail) +- `GGA_X_WC` (Wu–Cohen; PBE rational with `x = 10/81·s² + (μ−10/81)s²·e^(−s²) + ln(1+c·s⁴)`), `GGA_X_SOGGA`/`SOGGA11` (two-term second-order GGA), `GGA_X_C09X`, `GGA_X_HTBS` (switches between RPBE & PBEsol), `GGA_X_EV93` (Engel–Vosko rational poly in s²), `GGA_X_AK13`, `GGA_X_KT1`, `GGA_X_FT97_A/B`, `GGA_X_BAYESIAN`, `GGA_X_BCGP`, the `GGA_X_2D_*` 2D forms, and the Tozer–Handy XC forms. + +--- + +### GGA RUNG — CORRELATION + +#### Kernel G-C1: PBE/PW91 correlation H(t) gradient correction +- **Variables:** ρ, σ, ζ. **Reference:** PBE correlation (`GGA_C_PBE`). +- **Form:** `εc = εc^LDA(PW92) + H(rs,ζ,t)`, with `H = γφ³·ln[1 + (β/γ)t²·(1+At²)/(1+At²+A²t⁴)]`, `A = (β/γ)[exp(−εc^LDA/(γφ³))−1]^(−1)`, `φ = ((1+ζ)^(2/3)+(1−ζ)^(2/3))/2`, t = scaled gradient. +- **Member / parameters:** + + | Functional | β | γ | + |---|---|---| + | `GGA_C_PBE` | 0.06672455060314922 | (1−ln2)/π² ≈ 0.0310906909 | + | `GGA_C_PBE_SOL` | 0.046 | (1−ln2)/π² | + | `GGA_C_PBEFE` | 0.043 | (1−ln2)/π² | + | `GGA_C_PBE_JRGX` | 3μ^GE/π² ≈ 0.0375 | (1−ln2)/π² | + | `GGA_C_PBELOC`, `GGA_C_SG4`, `GGA_C_PBEINT` | β = β(rs,t) density-dependent | — | + + Members (all H(t) form, parameter-only): `GGA_C_PBE`, `PBE_SOL`, `XPBE`, `PBEFE`, `PBE_MOL`, `PBE_JRGX`, `APBE`, `PBEINT`, `PBELOC`, `SG4`, `REGTPSS`, `ZPBESOL`. +- **Citation:** Perdew, Burke & Ernzerhof, *PRL* 77, 3865 (1996); PW91: Perdew et al., *PRB* 46, 6671 (1992). + +#### Kernel G-C2: LYP correlation (Colle–Salvetti-derived) +- **Variables:** ρ, σ (originally ∇²ρ, eliminated by parts). **Form (Miehlich):** `εc = −a·(4/(1+d·ρ^(−1/3)))·(ρ↑ρ↓/ρ) − a·b·ω(ρ)·[...gradient and ρ^(8/3) terms...]`, `ω(ρ) = exp(−c·ρ^(−1/3))/(1+d·ρ^(−1/3))·ρ^(−11/3)`. +- **Parameters (libxc `gga_c_lyp_set_params`):** a = 0.04918, b = 0.132, c = 0.2533, d = 0.349. +- **Members:** `GGA_C_LYP`, plus refits used in OLYP, XLYP, MOHLYP, PBELYP1W. **Citation:** Lee, Yang & Parr, *Phys. Rev. B* 37, 785 (1988); Miehlich et al., *Chem. Phys. Lett.* 157, 200 (1989). Irreducible kernel. + +#### Kernel G-C3: P86 (Perdew 1986) gradient correlation +- **Form:** `εc = εc^LDA + e^(−Φ)·C(rs)·|∇ρ|²/ρ^(4/3)` over a spin factor d(ζ); C(rs) is a Rasolt–Geldart rational. **Members:** `GGA_C_P86`, `GGA_C_P86VWN`, `GGA_C_P86VWN_FT`, `GGA_C_P86_FT`. Differ by underlying LDA (PZ vs VWN5) and the ftilde constant. **Citation:** Perdew, *Phys. Rev. B* 33, 8822 (1986). + +#### Kernel G-C4: Wilson–Levy type +- **Form:** `εc = √(1−ζ²)·(a+b·x)/(c+d(x↑+x↓)+rs)`. **Members:** `GGA_C_WL` (a=−0.74860, b=0.06001, c=3.60073, d=0.90000), `GGA_C_WI`, `GGA_C_WI0`. **Citation:** Wilson & Levy, *Phys. Rev. B* 41, 12930 (1990). + +#### Misc GGA correlation one-offs +- `GGA_C_AM05` (paired with AM05-X), `GGA_C_OPTC` (Cohen–Handy), `GGA_C_W94` (Wilson 94), `GGA_C_CS1`, `GGA_C_BMK`, `GGA_C_CCDF`, `GGA_C_SOGGA11`, `GGA_C_GAPC`, `GGA_C_TCA`. + +--- + +### GGA RUNG — KINETIC ENERGY (orbital-free) + +#### Kernel G-K1: TF–λ–vW linear combination +- **Form:** `T = λ·TW + γ·TTF`. **Members (by λ):** `GGA_K_VW` (pure von Weizsäcker), `GGA_K_GE2` (λ=1/9), `GGA_K_GOLDEN` (λ=13/45), `GGA_K_YT65` (λ=1/5), `GGA_K_BALTIN` (λ=5/9), `GGA_K_TFVW`. Parameter-only (the λ coefficient). **Citation:** von Weizsäcker, *Z. Phys.* 96, 431 (1935). + +#### Kernel G-K2: PBE-type kinetic enhancement (reuses the G-X1 shell) +- **Form:** `Ft(s) = 1 + κ − κ/(1+μs²/κ)` applied to τTF. **Members:** `GGA_K_APBE`, `GGA_K_APBEINT`, `GGA_K_REVAPBE`, `GGA_K_LC94`, `GGA_K_PERDEW`, `GGA_K_TW1–4` (Tran–Wesolowski), `GGA_K_LLP`, `GGA_K_THAKKAR`, `GGA_K_OL1/OL2`, `GGA_K_FR_PW86`, `GGA_K_DK` (DePristo–Kress rational), `GGA_K_VSK`, `GGA_K_VJKS`, `GGA_K_ERNZERHOF`. Most are PBE/B88-shell forms refit for the kinetic energy. **Citation:** Constantin et al., *PRL* 106, 186406 (2011) (APBE). + +--- + +### META-GGA RUNG — EXCHANGE + +#### Kernel M-X1: TPSS-type rational shell with iso-orbital argument +- **Variables:** ρ, σ, τ. **Reference:** TPSS (`MGGA_X_TPSS`). +- **Form:** `Fx = 1 + κ − κ²/(κ + x(s,z,α))`, where x is a long polynomial in s, `z = τW/τ`, `α = (τ−τW)/τunif`, and the reduced-Laplacian mimic `q̃b = (9/20)(α−1)/√(1+b·α(α−1)) + (2/3)s²`. +- **Member / parameters:** `MGGA_X_TPSS` (κ=0.804, μGE=10/81, μPBE=0.21951, b=0.40, c=1.59096, e=1.537), `MGGA_X_REVTPSS` (μ=0.14, c=2.35203946, e=2.16769874, f=3), `MGGA_X_TPSSLOC`, `MGGA_X_BLOC`, `MGGA_X_REGTPSS`, `MGGA_X_PKZB` (precursor), `MGGA_X_MBEEF`, `MGGA_X_TM` (Tao–Mo), `MGGA_X_REVTM`. **Citation:** Tao, Perdew, Staroverov & Scuseria, *PRL* 91, 146401 (2003). + +#### Kernel M-X2: SCAN-type iso-orbital interpolation +- **Variables:** ρ, σ, τ (via α). **Reference:** SCAN (`MGGA_X_SCAN`). +- **Form:** `Fx(s,α) = {h1x(s,α) + fx(α)[h0x − h1x]}·gx(s)`, with `h0x = 1.174 = 1+κ0`, `h1x = 1 + k1 − k1/(1 + x/k1)`, `gx(s) = 1 − exp(−a1·s^(−1/2))`, `fx(α) = exp(−c1x·α/(1−α))·θ(1−α) − dx·exp(c2x/(1−α))·θ(α−1)`. +- **Parameters:** SCAN: μGE=10/81, k1=0.065, c1x=0.667, c2x=0.8, dx=1.24, a1=4.9479, with b1, b2≈0.12083, b3=0.5, b4 fixed by GE4. revSCAN: c1x=0.607, c2x=0.7, dx=1.37. r2SCAN regularizes via `α̃` with floor η=10^(−3): per Furness et al. (2020), *"the constants Cη = 20/27 + 5η/3, depending on the α̅ regularization parameter η = 10⁻³, and C2x ≈ −0.162742 eliminate erroneous contributions from dfx(α̅)/dα̅ at α̅ → 1."* The damping constant is **dp2 = 0.361** — note the original paper's main text printed 0.316, corrected to 0.361 in the published Correction (*J. Phys. Chem. Lett.* 11, 9248). rSCAN replaces the θ-step switches with smooth polynomials. +- **Members:** `MGGA_X_SCAN`, `MGGA_X_RSCAN`, `MGGA_X_R2SCAN`, `MGGA_X_R2SCAN01`, `MGGA_X_R4SCAN`, `MGGA_X_SCANL` (deorbitalized, Laplacian instead of τ), `MGGA_X_REVSCAN`, `MGGA_X_TASK`, `MGGA_X_RPPSCAN`. **Citation:** Sun, Ruzsinszky & Perdew, *PRL* 115, 036402 (2015); Furness, Kaplan, Ning, Perdew & Sun, *J. Phys. Chem. Lett.* 11, 8208 (2020), doi:10.1021/acs.jpclett.0c02405 (r2SCAN). **Stability note:** SCAN's `fx(α)` step switches are highly grid-sensitive near α≈1; rSCAN/r2SCAN explicitly regularize this (`α̃` with the η floor in density tails where τW/τ diverges). + +#### Kernel M-X3: Minnesota (M06/M11/MN) PBE×power-series + VS98 term +- **Variables:** ρ, σ, τ. **Reference:** M06-L (`MGGA_X_M06_L`). +- **Form:** `εx = Σσ [εx^PBE(ρσ,∇ρσ)·f(wσ) + εx^UEG(ρσ)·hx(xσ,zσ)]`, with `f(wσ) = Σ_{i=0}^{11} aᵢ·wσⁱ` a power series in `wσ = (t−1)/(t+1)`, `t = τLSDA/τ`; hx is the VS98 inhomogeneity rational. +- **Members / parameters:** the aᵢ vectors (and the VS98 a,b,c constants) differ per functional: `MGGA_X_M06_L`, `M06`, `M06_2X`, `M06_HF`, `M06_SX`, `M05`, `M05_2X`, `M08_HX`, `M08_SO`, `M11`, `M11_L`, `MN12_L`, `MN15_L`, `MN12_SX`, `MN15`. Parent: `MGGA_X_VS98`; also `MGGA_X_GVT4`. **Citation:** Zhao & Truhlar, *J. Chem. Phys.* 125, 194101 (2006); Van Voorhis & Scuseria, *J. Chem. Phys.* 109, 400 (1998). **Stability note:** Minnesota Fx shape is basis-set-dependent and slowly grid-convergent (Mardirossian & Head-Gordon, *JCTC* 9, 4453 (2013)); `MGGA_X_MK00` has a `τ−υ/4` denominator that can vanish or go negative. + +#### Kernel M-X4: Becke–Roussel exchange-hole inversion +- **Variables:** ρ, σ, τ (and ∇²ρ). **Form:** models the spherically-averaged exchange hole as a displaced exponential, solving a transcendental relation `x·e^(−x)/… = (curvature)` for the nonlinear hole parameter; energy density `= −(1−e^(−x)−½x·e^(−x))/b`. Requires a Newton root-find per grid point (or the analytic Proynov approximation). +- **Members:** `MGGA_X_BR89`, `MGGA_X_BR89_EXPLICIT` (Proynov–Gan–Kong 2008 analytic refit, γ usually 0.8 or 1.0), `MGGA_X_BR89_EXPLICIT_1`, `MGGA_X_B00`, `MGGA_X_MBR`, and the mBJ (Tran–Blaha modified Becke–Johnson) potential `MGGA_X_TB09`/`MBJ` (potential-only). **Citation:** Becke & Roussel, *Phys. Rev. A* 39, 3761 (1989); Proynov et al., *Chem. Phys. Lett.* 455, 103 (2008). **Stability note:** the explicit refit raises quadrature noise ~3 orders of magnitude vs the iterative BR89 (Lehtola et al., arXiv:2206.14062); needs the Laplacian and is affected by libxc's Fermi-hole-curvature handling (`--disable-fhc`). + +#### Kernel M-X5: MVS / MS (made-simple) interpolation +- **Form:** PBE-rational h-function interpolated via a *rational* fx(α) (no step functions). **Members:** `MGGA_X_MS0`, `MS1`, `MS2`, `MS2B`, `MVS`, `MGGA_X_MVSB`. **Citation:** Sun et al., *J. Chem. Phys.* 137, 051101 (2012); Sun, Perdew & Ruzsinszky, *PNAS* 112, 685 (2015). + +#### meta-GGA exchange one-offs +- `MGGA_X_PKZB`, `MGGA_X_2D_PRHG07`, `MGGA_X_GX`/`PBE_GX` (Loos; step-function, numerically ill-behaved per Lehtola 2022), `MGGA_X_TLDA`, `MGGA_X_EDMGGA`, `MGGA_X_MK00`/`MK00B`, `MGGA_X_RTPSS`, `MGGA_X_JK`. + +--- + +### META-GGA RUNG — CORRELATION + +#### Kernel M-C1: TPSS/revTPSS correlation (self-interaction-corrected PBE) +- **Form:** `εc^TPSS = εc^revPKZB·[1 + 2.8·z³·εc^revPKZB]`, with `εc^revPKZB = εc^PBE[1+C(ζ,ξ)z²] − [1+C(ζ,ξ)]z²·Σσ (ρσ/ρ)·max(εc^PBE,σ, εc^PBE)`, `z = τW/τ`. **Members:** `MGGA_C_TPSS` (enforces z = τW/τ ≤ 1), `MGGA_C_REVTPSS`, `MGGA_C_TPSSLOC`, `MGGA_C_BLOC`. **Citation:** Tao–Perdew–Staroverov–Scuseria, *PRL* 91, 146401 (2003). **Stability note:** libxc applies a `τW/τ ≤ 1` clamp; a historical factor-of-½-in-τ interface bug caused NaNs in some host codes (QE issue). + +#### Kernel M-C2: SCAN correlation +- **Form:** `εc = εc^1 + fc(α)[εc^0 − εc^1]`, interpolating a single-orbital `εc^0` (revised PW92-based) and a slowly-varying `εc^1` (PBE-like with β(rs)); fc(α) analogous to the exchange switch. **Members:** `MGGA_C_SCAN`, `MGGA_C_RSCAN`, `MGGA_C_R2SCAN`, `MGGA_C_R2SCAN01`, `MGGA_C_SCANL`, `MGGA_C_REVSCAN`, `MGGA_C_R2SCAN_VV10` (composition with VV10 NL). **Citation:** Sun et al., *PRL* 115, 036402 (2015). + +#### Kernel M-C3: M06-L / Minnesota correlation (BC95-derived) +- **Form:** opposite-spin + same-spin components each `= εc^UEG·g(x)` (a B97 power series in x) times a τ-dependent factor (`1 − τW/τ` for same-spin self-correlation removal), built on the BC95 (Becke 1995) τ-correlation. **Members:** `MGGA_C_M06_L`, `M06`, `M06_2X`, `M06_HF`, `M05`, `M05_2X`, `M08_HX/SO`, `M11`, `M11_L`, `MN12_L`, `MN15_L`, `MN15`, `MGGA_C_BC95` (parent), `MGGA_C_VSXC`. **Citation:** Becke, *J. Chem. Phys.* 104, 1040 (1996); Zhao & Truhlar. + +#### Kernel M-C4: Colle–Salvetti / B94 τ-correlation +- **Members:** `MGGA_C_CS` (Colle–Salvetti, the parent of LYP), `MGGA_C_B94`, `MGGA_C_KCIS`, `MGGA_C_TPSSLOC`. Distinct from LYP only in retaining τ. + +#### meta-GGA correlation one-offs +- `MGGA_C_PKZB`, `MGGA_C_M12C`, `MGGA_C_TM`, `MGGA_C_REVTM`, `MGGA_C_HLTAPW`, `MGGA_C_RREGTM`. + +#### meta-GGA kinetic +- `MGGA_K_PC07` (Perdew–Constantin Laplacian-level), `MGGA_K_CSK`, `MGGA_K_PGSL025` — Laplacian + τ kinetic functionals. + +--- + +### HYBRID & RANGE-SEPARATED LAYER (thin recipe layer — mixing data only) + +Exact exchange is host-supplied; libxc only stores the mixing coefficients and the semilocal kernel composition. The CAM/Ewald range partition convention is: + +`1/r₁₂ = [α + β·erf(ωr₁₂)]/r₁₂ (long-range) + [1 − α − β·erf(ωr₁₂)]/r₁₂ (short-range)`, + +so `EX = α·EX^HF(full) + β·EX^HF(LR,erf) + (1−α)EX^DFA(SR) + …`. As r₁₂→0 the HF fraction = α; as r₁₂→∞ the HF fraction = α+β. + +- **Global hybrids (single α):** `HYB_GGA_XC_PBEH`/PBE0 (α=0.25, PBE X+C), `HYB_GGA_XC_B3LYP` (α=0.20 HF; 0.08 LDA-X + 0.72 B88-X; 0.81 LYP-C + 0.19 VWN-C — note the VWN_RPA vs VWN5 ambiguity: B3LYP uses VWN3/RPA in Gaussian/ORCA/PySCF≥2.3, VWN5 in `B3LYP5`), `HYB_GGA_XC_B3LYP3` (VWN3), `B3LYP5` (VWN5), `BHANDH` (α=0.5), `BHANDHLYP`, `B3PW91`, `B3P86`, `X3LYP`, `O3LYP`, `PBE0_13` (α=1/3); `HYB_MGGA_XC_TPSSH` (α=0.10), `revTPSSH`, `SCAN0` (α=0.25); MGGA hybrids `B86B95`/`B88B95` (Becke τ-correlation, α per functional), and the M05/M06/M08/M11/MN15 hybrid members (HF % from 10% to 54%, and 100% for M06-HF). +- **Range-separated (α, β, ω):** + + | Functional | α | β | ω (bohr⁻¹) | convention | + |---|---|---|---|---| + | `HYB_GGA_XC_CAM_B3LYP` | 0.19 | 0.46 | 0.33 | CAM | + | LC-ωPBE (`LRC_WPBE`) | 0 | 1.0 | ~0.4 | LC | + | LC-ωPBEh (`LRC_WPBEH`) | 0.20 | 0.80 | ~0.2 | CAM | + | `HSE03` / `HSE06` | 0.25 (SR) | −0.25 | 0.11 (HSE06) | screened (HF SR, DFA LR) | + | `ωB97` | 0 | 1.0 | 0.40 | LC + B97 series | + | `ωB97X` | 0.157706 (SR) | →1.0 | 0.30 | CAM + B97 series | + | `ωB97X-V` / `ωB97M-V` | as ωB97X | | | + VV10 NL correlation | + | `CAMY-B3LYP` / `LCY-*` | per Yukawa | | | Yukawa screening (not erf) | + + Historical libxc bugs to be aware of: `LC_BLYP` used ω=0.3 vs intended 0.33; `CAM_QTP_01` had 81/19 vs intended 80/20 LYP/VWN5; `CAP0` had 75% vs intended 25% HF. +- **Yukawa range separation** (screening `exp(−ωr)/r` instead of erf): the CAMY / LCY / LC-Yukawa family — a different attenuation kernel, same mixing-data treatment. +- **Composition (recipe, not kernel):** many `_XC` entries are literally an X-kernel + C-kernel glued together: BLYP = B88 + LYP, BP86 = B88 + P86, OLYP = OPTX + LYP, PBE = PBE-X + PBE-C. These are recipes (two kernel references + coefficients), not new kernels. + +--- + +## Compression Summary + +- **600+ libxc entries → ~50–70 distinct analytic kernels.** Rough kernel counts by rung: **LDA ≈ 12** (Dirac-X + relativistic/dim variants; VWN, PW92, PZ, Chachiyo, Wigner, HL/GL, finite-T KSDT/GDSMFB, EPC, TF-kinetic); **GGA-X ≈ 12–15** (PBE-rational, RPBE-exp, B88, B86, PW91, OPTX, G96, AM05, LB-potential, HJS, WC, SOGGA, EV93, AK13); **GGA-C ≈ 6** (PBE-H, LYP, P86, Wilson–Levy, AM05, one-offs); **GGA-K ≈ 2** (TF-λ-vW, PBE-kinetic-shell); **mGGA-X ≈ 6** (TPSS, SCAN, Minnesota, BR89, MVS/MS, one-offs); **mGGA-C ≈ 5** (TPSS, SCAN, Minnesota/BC95, CS/B94); **hybrids/RSH = 0 new kernels** (pure mixing data). +- **The single biggest win:** the PBE rational shell `Fx = 1+κ−κ/(1+μs²/κ)` plus its argument-substitution descendants (TPSS, SCAN h1x) covers well over 100 functional entries across GGA-X, GGA-K and mGGA-X with parameter/argument changes only. + +--- + +## Recommendations (for the JAX clean-room reimplementation) + +1. **Stage 1 — implement the high-coverage kernels first.** In priority order by functional count: (a) Dirac–Slater LDA-X, (b) PW92 + VWN local correlation, (c) PBE rational exchange shell, (d) PBE H(t) correlation, (e) B88 + LYP (gets you BLYP/B3LYP), (f) the B97 u-power series. These six kernels plus their parameter tables already cover the majority of *production* functionals (PBE, PBEsol, revPBE, RPBE, BLYP, B3LYP, PBE0, the B97/HCTH family). +2. **Stage 2 — meta-GGA shells:** TPSS rational-with-iso-orbital, SCAN/r2SCAN interpolation (implement r2SCAN's regularized `α̃` first — it is the numerically robust default; use dp2 = 0.361), Minnesota PBE×power-series. Defer BR89 (it needs a per-point Newton root-find for the exchange-hole parameter — implement as a custom-JVP-wrapped solver). +3. **Represent parameters as a database, kernels as pure functions.** Schema per leaf: `{functional_id, kernel_family_id, param_tuple, base_LDA_ref, family/kind tags, citation/DOI}`. Hybrids/RSH get an extra mixing record `{alpha, beta, omega, screening_convention ∈ {erf, yukawa}, component_kernel_refs}`. +4. **Encode the numerical guards from libxc source as kernel-level config:** density floor (`dens_threshold`), s→∞ saturation (PBE rational saturates naturally; RPBE/B88 need overflow-safe exp/arcsinh), τ iso-orbital boundaries (α near 1 for SCAN — use r2SCAN regularization; z ≤ 1 clamp for TPSS correlation), and the Laplacian/Fermi-hole-curvature handling for the BR89 class (replicate `--disable-fhc` behavior). +5. **Validate against libxc reference data per functional** (energies + up to 4th derivatives), since libxc's Maple-generated code is the de-facto numerical standard. Mismatches usually trace to (i) the VWN3/VWN5 ambiguity, (ii) the τ factor-of-½ convention, or (iii) the (3π²) vs (6π²) reduced-gradient constants. + +**Thresholds that change the plan:** if you only need the ~20 most-used functionals, Stages 1–2 minus BR89/Minnesota suffice (~10 kernels). If you need full libxc parity including the long tail, budget for ~60–70 kernels and the transcendental BR89 solver; the marginal cost per additional functional then drops to a single database row. + +--- + +## Caveats + +- **Parameter precision:** three values were corroborated from the peer-reviewed Tran–Lehtola–Pittalis–Marques review (arXiv:2602.17333, 2026) and program docs rather than read verbatim from the `.mpl`/`.c` source blobs — specifically `GGA_X_XPBE` (κ=0.91954, μ=0.23889), `GGA_X_PBE_MOL` (μ≈0.27583), and the `GGA_X_LAMBDA_OC2_N` N-dependent κ. Cross-check these directly against `maple/gga_x_xpbe.mpl`, `gga_x_pbe.c`, and `gga_x_lambda*.mpl` before committing to a database. (GitLab raw blobs render via JavaScript and were not directly fetchable during this research.) The PBE/PBEsol/B88/LYP/PW92/SCAN/r2SCAN constants are confirmed (PBEsol μ=10/81, β=0.046 corroborated by the VASP wiki; r2SCAN dp2=0.361 per the published Correction to Furness 2020). +- **Kernel count is approximate** and depends on how finely you split "same form, different power p" cases (e.g. B86 p=1 vs p=4/5, or treating RPBE-exponential as distinct from PBE-rational). The 50–70 figure treats distinct analytic *expressions* as distinct; collapsing near-identical rationals would push it lower. +- **Potential-only functionals** (LB94, mBJ/TB09) do not derive from an energy expression and must be flagged specially — they return a vxc directly and have no consistent εxc; libxc lists them but they cannot be used in energy-only or post-SCF energy evaluations. +- **Version drift:** functional definitions and even default parametrizations change between libxc releases (the changelog documents corrected parameters for KCIS, OPTC, CAM_QTP_01, LC_BLYP ω, and others). Pin to a specific libxc version (devel as of mid-2026) when building the parameter database. +- **VWN3/VWN5 and PZ/PW92 LDA-base ambiguity** is the single most common source of cross-code disagreement; record the exact LDA-correlation base for every GGA/mGGA correlation functional and every hybrid explicitly. +- **Combined-XC vs split functionals:** in libxc, most functionals are split into separate `_X` and `_C` entries that the host combines; a minority (LYP-based, Tozer–Handy, the `_XC` Minnesota and range-separated forms) are delivered as monolithic XC. Your database should mark each leaf as exchange-only, correlation-only, or combined-XC to drive the composition layer correctly. diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..8ca27dc --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,20 @@ +[project] +name = "funxc" +version = "0.1.0" +description = "Differentiable exchange-correlation functionals in JAX: a clean-room libxc reimplementation" +readme = "README.md" +requires-python = ">=3.13" +dependencies = ["jax>=0.5", "numpy>=2.0"] + +[dependency-groups] +dev = ["pytest>=8"] + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["src/funxc"] + +[tool.pytest.ini_options] +testpaths = ["tests"] diff --git a/scratch_scan_cmp.py b/scratch_scan_cmp.py new file mode 100644 index 0000000..f4cbf31 --- /dev/null +++ b/scratch_scan_cmp.py @@ -0,0 +1,30 @@ +import sys, numpy as np +import jax; jax.config.update("jax_enable_x64", True) +sys.path.insert(0, "tests") +sys.path.insert(0, "libxc") +from libxc_reference import load_input, load_regression +import funxc + +inp = load_input("Li") +rho, sigma, lapl, tau = inp["rho"], inp["sigma"], inp["lapl"], inp["tau"] +ref1 = load_regression("mgga_x", "mgga_x_scan", "Li", "pol", 1) + +# installed library +from pylibxc import LibXCFunctional +lf = LibXCFunctional("MGGA_X_SCAN", "polarized") +out = lf.compute({ + "rho": np.ascontiguousarray(rho.reshape(-1)), + "sigma": np.ascontiguousarray(sigma.reshape(-1)), + "lapl": np.ascontiguousarray(lapl.reshape(-1)), + "tau": np.ascontiguousarray(tau.reshape(-1)), +}) +lib_vrho = out["vrho"].reshape(-1,2) + +f = funxc.functional("MGGA_X_SCAN", polarized=True) +fout = f.exc_vxc(rho, sigma, lapl, tau) +fx_vrho = np.asarray(fout["vrho"]) + +print("idx golden_vrhoB library_vrhoB funxc_vrhoB") +for i in [3,5]: + print(i, f"{ref1['vrho(b)'][i]: .10e} {lib_vrho[i,1]: .10e} {fx_vrho[i,1]: .10e}") + print(" rho:", rho[i], "sigma:", sigma[i], "tau:", tau[i]) diff --git a/scratch_scan_probe.py b/scratch_scan_probe.py new file mode 100644 index 0000000..8b910eb --- /dev/null +++ b/scratch_scan_probe.py @@ -0,0 +1,28 @@ +import sys, numpy as np +import jax; jax.config.update("jax_enable_x64", True) +sys.path.insert(0, "libxc") +import funxc +from pylibxc import LibXCFunctional + +# spread of points: vary rho, sigma, tau to hit small/moderate/large s and alpha +rho = np.array([[0.5,0.5],[0.3,0.2],[1.0,0.8],[0.1,0.1],[2.0,1.0],[0.05,0.04]]) +sigma = np.array([[0.1,0.05,0.1],[0.02,0.01,0.02],[0.5,0.3,0.4], + [0.001,0.0005,0.001],[3.0,2.0,2.5],[1e-4,5e-5,1e-4]]) +tau = np.array([[0.3,0.3],[0.2,0.15],[0.9,0.7],[0.05,0.05],[1.5,0.8],[0.02,0.015]]) +lapl = np.zeros_like(tau) + +lf = LibXCFunctional("MGGA_X_SCAN","polarized") +o = lf.compute({"rho":np.ascontiguousarray(rho.reshape(-1)), + "sigma":np.ascontiguousarray(sigma.reshape(-1)), + "lapl":np.ascontiguousarray(lapl.reshape(-1)), + "tau":np.ascontiguousarray(tau.reshape(-1))}) +lib_zk=o["zk"].reshape(-1); lib_vrho=o["vrho"].reshape(-1,2) +f=funxc.functional("MGGA_X_SCAN",polarized=True) +fx=f.exc_vxc(rho,sigma,lapl,tau) +fz=np.asarray(fx["zk"]).reshape(-1); fv=np.asarray(fx["vrho"]) +print("pt lib_zk fx_zk relerr_zk relerr_vrhoA relerr_vrhoB") +for i in range(len(rho)): + rz=(fz[i]-lib_zk[i])/lib_zk[i] + ra=(fv[i,0]-lib_vrho[i,0])/lib_vrho[i,0] + rb=(fv[i,1]-lib_vrho[i,1])/lib_vrho[i,1] + print(f"{i} {lib_zk[i]: .6e} {fz[i]: .6e} {rz: .2e} {ra: .2e} {rb: .2e}") diff --git a/scratch_scan_zk.py b/scratch_scan_zk.py new file mode 100644 index 0000000..57649d3 --- /dev/null +++ b/scratch_scan_zk.py @@ -0,0 +1,34 @@ +import sys, numpy as np +import jax; jax.config.update("jax_enable_x64", True) +sys.path.insert(0, "tests"); sys.path.insert(0, "libxc") +from libxc_reference import load_input +import funxc +from pylibxc import LibXCFunctional + +inp = load_input("Li") +rho, sigma, lapl, tau = inp["rho"], inp["sigma"], inp["lapl"], inp["tau"] +lf = LibXCFunctional("MGGA_X_SCAN", "polarized") +o = lf.compute({"rho":np.ascontiguousarray(rho.reshape(-1)), + "sigma":np.ascontiguousarray(sigma.reshape(-1)), + "lapl":np.ascontiguousarray(lapl.reshape(-1)), + "tau":np.ascontiguousarray(tau.reshape(-1))}) +lib_zk = o["zk"].reshape(-1) +f = funxc.functional("MGGA_X_SCAN", polarized=True) +fx = f.exc_vxc(rho, sigma, lapl, tau) +fx_zk = np.asarray(fx["zk"]).reshape(-1) +for i in [3,5]: + print(i, f"lib_zk={lib_zk[i]:.16e} fx_zk={fx_zk[i]:.16e} reldiff={(fx_zk[i]-lib_zk[i])/lib_zk[i]:.2e}") + +# Also finite-difference libxc zk*rho wrt rho_dn to get vrho(b) and compare +def e_of(rd, i): + r = rho.copy(); r[i,1]=rd + oo = lf.compute({"rho":np.ascontiguousarray(r.reshape(-1)), + "sigma":np.ascontiguousarray(sigma.reshape(-1)), + "lapl":np.ascontiguousarray(lapl.reshape(-1)), + "tau":np.ascontiguousarray(tau.reshape(-1))}) + n = r[i,0]+r[i,1] + return oo["zk"].reshape(-1)[i]*n +for i in [3,5]: + rd=rho[i,1]; h=rd*1e-6 + fd=(e_of(rd+h,i)-e_of(rd-h,i))/(2*h) + print(i,"lib_FD_vrhoB=",fd) diff --git a/src/funxc/__init__.py b/src/funxc/__init__.py new file mode 100644 index 0000000..62320ad --- /dev/null +++ b/src/funxc/__init__.py @@ -0,0 +1,25 @@ +"""funxc: differentiable exchange-correlation functionals in JAX. + +A clean-room reimplementation of libxc as composable, differentiable math +kernels parameterized by a TOML functional database. Numeric compatibility +with libxc 7.0.0 is the acceptance criterion. +""" + +# libxc parity is meaningless in float32; this must run before any jax +# arrays are created, so funxc must be imported before running computations. +import jax + +jax.config.update("jax_enable_x64", True) + +from .api import Functional, functional # noqa: E402 +from .libxc_compat import LibXCFunctional # noqa: E402 +from .registry import REGISTRY, FunctionalSpec, get_spec # noqa: E402 + +__all__ = [ + "Functional", + "functional", + "LibXCFunctional", + "REGISTRY", + "FunctionalSpec", + "get_spec", +] diff --git a/src/funxc/api.py b/src/funxc/api.py new file mode 100644 index 0000000..5c618de --- /dev/null +++ b/src/funxc/api.py @@ -0,0 +1,88 @@ +"""JAX-native public API. + + import funxc + f = funxc.functional("GGA_X_PBE") # polarized by default + zk = f.exc(rho, sigma) # rho (N,2), sigma (N,3) + out = f.exc_vxc(rho, sigma) # dict: zk, vrho, vsigma + +The per-point functions are also exposed (`f.eps`, `f.energy_density`) so +users can compose them with jax.grad / vmap / jit however they like — e.g. +differentiating a total energy directly with respect to density parameters. +""" + +import functools + +import jax +import jax.numpy as jnp + +from .evaluate import make_energy_density, make_eps +from .registry import get_spec + + +class Functional: + """A functional bound to a spin mode, with batched jit'd evaluators.""" + + # derivative names per family, in energy_density argument order + _DERIVS = {"lda": ("vrho",), "gga": ("vrho", "vsigma"), + "mgga": ("vrho", "vsigma", "vlapl", "vtau")} + + def __init__(self, name: str, polarized: bool = True): + self.spec = get_spec(name) + self.polarized = polarized + if self.spec.family not in self._DERIVS: + raise NotImplementedError(f"family {self.spec.family!r} not implemented yet") + + # Per-point pure functions (compose freely with jax transforms). + self.eps = make_eps(self.spec, polarized) + self.energy_density = make_energy_density(self.spec, polarized) + + self._derivs = self._DERIVS[self.spec.family] + nargs = len(self._derivs) + in_axes = (0,) * nargs + + self._exc = jax.jit(jax.vmap(self.eps, in_axes=in_axes)) + argnums = 0 if nargs == 1 else tuple(range(nargs)) + grad_e = jax.grad(self.energy_density, argnums=argnums) + self._vxc = jax.jit(jax.vmap(grad_e, in_axes=in_axes)) + + def exc(self, rho, sigma=None, lapl=None, tau=None): + """zk: XC energy per particle, shape (N,).""" + args = self._args(rho, sigma, lapl, tau) + return self._exc(*args) + + def exc_vxc(self, rho, sigma=None, lapl=None, tau=None): + """dict with zk and the first derivatives (vrho[, vsigma[, vlapl, vtau]]).""" + args = self._args(rho, sigma, lapl, tau) + out = {"zk": self._exc(*args)} + v = self._vxc(*args) + if len(self._derivs) == 1: + v = (v,) + out.update(zip(self._derivs, v)) + return out + + __call__ = exc + + def _args(self, rho, sigma, lapl, tau): + rho = jnp.asarray(rho) + family = self.spec.family + if family == "lda": + return (rho,) + if sigma is None: + raise ValueError(f"{self.spec.id} is a {family.upper()}: sigma is required") + if family == "gga": + return (rho, jnp.asarray(sigma)) + if tau is None: + raise ValueError(f"{self.spec.id} is a meta-GGA: tau is required") + if lapl is None: + lapl = jnp.zeros_like(rho) # vlapl is then identically zero + return (rho, jnp.asarray(sigma), jnp.asarray(lapl), jnp.asarray(tau)) + + def __repr__(self): + mode = "polarized" if self.polarized else "unpolarized" + return f"Functional({self.spec.id}, {mode}, kernel={self.spec.kernel})" + + +@functools.cache +def functional(name: str, polarized: bool = True) -> Functional: + """Look up a functional by libxc-style name (e.g. "GGA_X_PBE").""" + return Functional(name, polarized) diff --git a/src/funxc/constants.py b/src/funxc/constants.py new file mode 100644 index 0000000..5a8e108 --- /dev/null +++ b/src/funxc/constants.py @@ -0,0 +1,32 @@ +"""Shared mathematical constants, mirroring libxc/maple/util.mpl.""" + +import math + +# Reduced-gradient conversion: s = X2S * x, with x = |grad rho_sigma| / rho_sigma^(4/3) +X2S = 1.0 / (2.0 * (6.0 * math.pi**2) ** (1.0 / 3.0)) + +# LDA exchange prefactor: eps_x = -X_FACTOR_C * ((1+z)/2)^(4/3) * n^(1/3) summed over spins +X_FACTOR_C = 3.0 / 8.0 * (3.0 / math.pi) ** (1.0 / 3.0) * 4.0 ** (2.0 / 3.0) + +# Thomas-Fermi kinetic prefactor (per spin channel, libxc convention) +K_FACTOR_C = 3.0 / 10.0 * (6.0 * math.pi**2) ** (2.0 / 3.0) + +# Wigner-Seitz radius: rs = RS_FACTOR * n^(-1/3) +RS_FACTOR = (3.0 / (4.0 * math.pi)) ** (1.0 / 3.0) + +# Gradient-expansion and PBE constants +MU_GE = 10.0 / 81.0 +MU_PBE = 0.06672455060314922 * math.pi**2 / 3.0 +KAPPA_PBE = 0.8040 + +# libxc screening defaults (functionals.c): dens_threshold comes from each +# functional's info struct (1e-15 for all Stage 1 functionals); +# zeta_threshold is DBL_EPSILON; tau_threshold is a global 1e-20 +# ("many functionals shamelessly divide by tau", work_mgga_inc.c). +DENS_THRESHOLD = 1e-15 +ZETA_THRESHOLD = 2.220446049250313e-16 +TAU_THRESHOLD = 1e-20 + +# Floor used only to keep masked-out branches finite for autodiff; never +# affects unmasked values. +TINY = 1e-100 diff --git a/src/funxc/data/b97.toml b/src/funxc/data/b97.toml new file mode 100644 index 0000000..f8eb4e3 --- /dev/null +++ b/src/funxc/data/b97.toml @@ -0,0 +1,172 @@ +# B97 power-series family (kernel gga_xc_b97). Coefficient tables verified +# against libxc 7.0.0 src/gga_xc_b97.c; all members use dens_threshold 1e-14 +# (xc_func_info structs). Hybrid members (B97, B97-1/2/3, SB98, ...) share +# this kernel but need the exact-exchange recipe layer; deferred until then. + +[GGA_XC_B97_D] +name = "Becke 97-D (for use with -D dispersion)" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 170 +dens_threshold = 1e-14 +citation = "S. Grimme, J. Comput. Chem. 27, 1787 (2006)" +doi = "10.1002/jcc.20495" +[GGA_XC_B97_D.parameters] +c_x = [1.08662, -0.52127, 3.25429, 0.0, 0.0] +c_ss = [0.2234, -1.56208, 1.94293, 0.0, 0.0] +c_ab = [0.69041, 6.3027, -14.9712, 0.0, 0.0] + +[GGA_XC_B97_3C] +name = "Becke 97-3c (Grimme, mTZVP + D3 + SRB)" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 327 +dens_threshold = 1e-14 +citation = "J. G. Brandenburg, C. Bannwarth, A. Hansen, and S. Grimme, J. Chem. Phys. 148, 064104 (2018)" +doi = "10.1063/1.5012601" +[GGA_XC_B97_3C.parameters] +c_x = [1.076616, -0.469912, 3.322442, 0.0, 0.0] +c_ss = [0.543788, -1.444420, 1.637436, 0.0, 0.0] +c_ab = [0.635047, 5.532103, -15.301575, 0.0, 0.0] + +[GGA_XC_B97_GGA1] +name = "Becke 97 GGA-1" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 96 +dens_threshold = 1e-14 +citation = "A. J. Cohen and N. C. Handy, Chem. Phys. Lett. 316, 160 (2000)" +doi = "10.1016/S0009-2614(99)01273-7" +[GGA_XC_B97_GGA1.parameters] +c_x = [1.1068, -0.8765, 4.2639, 0.0, 0.0] +c_ss = [0.4883, -2.117, 2.3235, 0.0, 0.0] +c_ab = [0.7961, 5.706, -14.982, 0.0, 0.0] + +[GGA_XC_HCTH_93] +name = "HCTH/93" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 161 +dens_threshold = 1e-14 +citation = "F. A. Hamprecht, A. J. Cohen, D. J. Tozer, and N. C. Handy, J. Chem. Phys. 109, 6264 (1998)" +doi = "10.1063/1.477267" +[GGA_XC_HCTH_93.parameters] +c_x = [1.0932, -0.744056, 5.5992, -6.78549, 4.49357] +c_ss = [0.222601, -0.0338622, -0.012517, -0.802496, 1.55396] +c_ab = [0.729974, 3.35287, -11.543, 8.08564, -4.47857] + +[GGA_XC_HCTH_120] +name = "HCTH/120" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 162 +dens_threshold = 1e-14 +citation = "A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik, J. Chem. Phys. 112, 1670 (2000)" +doi = "10.1063/1.480732" +[GGA_XC_HCTH_120.parameters] +c_x = [1.09163, -0.747215, 5.07833, -4.10746, 1.17173] +c_ss = [0.489508, -0.260699, 0.432917, -1.99247, 2.48531] +c_ab = [0.51473, 6.92982, -24.7073, 23.1098, -11.3234] + +# libxc note (src/gga_xc_b97.c): the css1 = 0.0171436 coefficient has the +# wrong sign in the original paper; corrected per doi:10.1063/1.1589004 and +# https://gitlab.com/libxc/libxc/-/issues/205. +[GGA_XC_HCTH_147] +name = "HCTH/147" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 163 +dens_threshold = 1e-14 +citation = "A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik, J. Chem. Phys. 112, 1670 (2000)" +doi = "10.1063/1.480732" +[GGA_XC_HCTH_147.parameters] +c_x = [1.09025, -0.799194, 5.57212, -5.8676, 3.04544] +c_ss = [0.562576, 0.0171436, -1.30636, 1.05747, 0.885429] +c_ab = [0.542352, 7.01464, -28.3822, 35.0329, -20.4284] + +[GGA_XC_HCTH_407] +name = "HCTH/407" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 164 +dens_threshold = 1e-14 +citation = "A. D. Boese and N. C. Handy, J. Chem. Phys. 114, 5497 (2001)" +doi = "10.1063/1.1347371" +[GGA_XC_HCTH_407.parameters] +c_x = [1.08184, -0.518339, 3.42562, -2.62901, 2.28855] +c_ss = [1.18777, -2.40292, 5.61741, -9.17923, 6.24798] +c_ab = [0.589076, 4.42374, -19.2218, 42.5721, -42.0052] + +[GGA_XC_HCTH_P14] +name = "HCTH p=1/4" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 95 +dens_threshold = 1e-14 +citation = "G. Menconi, P. J. Wilson, and D. J. Tozer, J. Chem. Phys. 114, 3958 (2001)" +doi = "10.1063/1.1342776" +[GGA_XC_HCTH_P14.parameters] +c_x = [1.03161, -0.360781, 3.51994, -4.95944, 2.41165] +c_ss = [2.82414, 0.0318843, -1.78512, 2.39795, -0.876909] +c_ab = [0.0821827, 4.56466, -13.5529, 13.382, -3.17493] + +[GGA_XC_HCTH_P76] +name = "HCTH p=7/6" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 94 +dens_threshold = 1e-14 +citation = "G. Menconi, P. J. Wilson, and D. J. Tozer, J. Chem. Phys. 114, 3958 (2001)" +doi = "10.1063/1.1342776" +[GGA_XC_HCTH_P76.parameters] +c_x = [1.16525, -0.583033, 2.51769, 3.81278, -5.45906] +c_ss = [-3.92143, -1.10098, -0.091405, -0.859723, 2.07184] +c_ab = [0.192949, -5.73335, 50.8757, -135.475, 101.268] + +[GGA_XC_HCTH_407P] +name = "HCTH/407+" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 93 +dens_threshold = 1e-14 +citation = "A. D. Boese, A. Chandra, J. M. L. Martin, and D. Marx, J. Chem. Phys. 119, 5965 (2003)" +doi = "10.1063/1.1599338" +[GGA_XC_HCTH_407P.parameters] +c_x = [1.08018, -0.4117, 2.4368, 1.389, -1.3529] +c_ss = [0.80302, -1.0479, 4.9807, -12.89, 9.6446] +c_ab = [0.73604, 3.027, -10.075, 20.611, -29.418] + +[GGA_XC_HLE16] +name = "High local exchange 2016" +maple_name = "gga_xc_b97" +family = "gga" +kind = "xc" +kernel = "gga_xc_b97" +libxc_id = 545 +dens_threshold = 1e-14 +citation = "P. Verma and D. G. Truhlar, J. Phys. Chem. Lett. 8, 380 (2017)" +doi = "10.1021/acs.jpclett.6b02757" +[GGA_XC_HLE16.parameters] +c_x = [1.3523, -0.64792375, 4.282025, -3.2862625, 2.8606875] +c_ss = [0.593885, -1.20146, 2.808705, -4.589615, 3.12399] +c_ab = [0.294538, 2.21187, -9.6109, 21.28605, -21.0026] diff --git a/src/funxc/data/gga.toml b/src/funxc/data/gga.toml new file mode 100644 index 0000000..510beca --- /dev/null +++ b/src/funxc/data/gga.toml @@ -0,0 +1,255 @@ +# GGA functionals. Parameter values verified against libxc 7.0.0 +# src/gga_x_pbe.c, gga_x_rpbe.c, gga_x_b88.c, gga_c_pbe.c, gga_c_lyp.c. +# +# mu_PBE = 0.06672455060314922 * pi^2 / 3 = 0.2195149727645171 +# mu_GE = 10/81 = 0.12345679012345678 +# gamma_PBE = (1 - ln 2)/pi^2 = 0.031090690869654895034 + +[GGA_X_PBE] +name = "Perdew, Burke & Ernzerhof exchange" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 101 +citation = "J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)" +doi = "10.1103/PhysRevLett.77.3865" +[GGA_X_PBE.parameters] +kappa = 0.8040 +mu = 0.2195149727645171 + +[GGA_X_PBE_R] +name = "revPBE: Perdew, Burke & Ernzerhof exchange, revised by Zhang & Yang" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 102 +citation = "Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998)" +doi = "10.1103/PhysRevLett.80.890" +[GGA_X_PBE_R.parameters] +kappa = 1.245 +mu = 0.2195149727645171 + +[GGA_X_PBE_SOL] +name = "PBEsol: Perdew, Burke & Ernzerhof exchange for solids" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 116 +citation = "J. P. Perdew et al., Phys. Rev. Lett. 100, 136406 (2008)" +doi = "10.1103/PhysRevLett.100.136406" +[GGA_X_PBE_SOL.parameters] +kappa = 0.804 +mu = 0.12345679012345678 + +[GGA_X_XPBE] +name = "xPBE exchange of Xu & Goddard" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 123 +citation = "X. Xu and W. A. Goddard III, J. Chem. Phys. 121, 4068 (2004)" +doi = "10.1063/1.1771632" +[GGA_X_XPBE.parameters] +kappa = 0.91954 +mu = 0.23214 + +[GGA_X_APBE] +name = "APBE exchange: mu fixed from the semiclassical neutral atom" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 184 +citation = "L. A. Constantin, E. Fabiano, S. Laricchia, and F. Della Sala, Phys. Rev. Lett. 106, 186406 (2011)" +doi = "10.1103/PhysRevLett.106.186406" +[GGA_X_APBE.parameters] +kappa = 0.8040 +mu = 0.260 + +[GGA_X_PBE_TCA] +name = "PBE revised by Tognetti, Cortona & Adamo" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 59 +citation = "V. Tognetti, P. Cortona, and C. Adamo, Chem. Phys. Lett. 460, 536 (2008)" +doi = "10.1016/j.cplett.2008.06.032" +[GGA_X_PBE_TCA.parameters] +kappa = 1.227 +mu = 0.2195149727645171 + +[GGA_X_PBE_MOL] +name = "PBEmol exchange of Del Campo, Gazquez, Trickey & Vela" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 49 +citation = "J. M. del Campo, J. L. Gazquez, S. B. Trickey, and A. Vela, J. Chem. Phys. 136, 104108 (2012)" +doi = "10.1063/1.3691197" +[GGA_X_PBE_MOL.parameters] +kappa = 0.8040 +mu = 0.27583 + +[GGA_X_PBEFE] +name = "PBEfe exchange for formation energies" +maple_name = "gga_x_pbe" +family = "gga" +kind = "x" +kernel = "gga_x_pbe" +libxc_id = 265 +citation = "R. Sarmiento-Perez, S. Botti, and M. A. L. Marques, J. Chem. Theory Comput. 11, 3844 (2015)" +doi = "10.1021/acs.jctc.5b00529" +[GGA_X_PBEFE.parameters] +kappa = 0.437 +mu = 0.346 + +[GGA_X_RPBE] +name = "RPBE exchange of Hammer, Hansen & Norskov" +maple_name = "gga_x_rpbe" +family = "gga" +kind = "x" +kernel = "gga_x_rpbe" +libxc_id = 117 +citation = "B. Hammer, L. B. Hansen, and J. K. Norskov, Phys. Rev. B 59, 7413 (1999)" +doi = "10.1103/PhysRevB.59.7413" +[GGA_X_RPBE.parameters] +kappa = 0.8040 +mu = 0.2195149727645171 + +[GGA_X_B88] +name = "Becke 88 exchange" +maple_name = "gga_x_b88" +family = "gga" +kind = "x" +kernel = "gga_x_b88" +libxc_id = 106 +citation = "A. D. Becke, Phys. Rev. A 38, 3098 (1988)" +doi = "10.1103/PhysRevA.38.3098" +[GGA_X_B88.parameters] +beta = 0.0042 +gamma = 6.0 + +[GGA_X_B88_6311G] +name = "Becke 88 exchange reoptimized with the 6-311G** basis set" +maple_name = "gga_x_b88" +family = "gga" +kind = "x" +kernel = "gga_x_b88" +libxc_id = 179 +citation = "J. M. Ugalde, C. Sarasola, and M. Aguado, J. Phys. Chem. 98, 4064 (1994)" +doi = "10.1021/j100066a026" +[GGA_X_B88_6311G.parameters] +beta = 0.0051 +gamma = 6.0 + +[GGA_C_PBE] +dens_threshold = 1e-12 +name = "Perdew, Burke & Ernzerhof correlation" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 130 +citation = "J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)" +doi = "10.1103/PhysRevLett.77.3865" +[GGA_C_PBE.parameters] +beta = 0.06672455060314922 +gamma = 0.031090690869654895034 +B = 1.0 + +[GGA_C_PBE_SOL] +dens_threshold = 1e-12 +name = "PBEsol correlation" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 133 +citation = "J. P. Perdew et al., Phys. Rev. Lett. 100, 136406 (2008)" +doi = "10.1103/PhysRevLett.100.136406" +[GGA_C_PBE_SOL.parameters] +beta = 0.046 +gamma = 0.031090690869654895034 +B = 1.0 + +[GGA_C_XPBE] +dens_threshold = 1e-12 +name = "xPBE correlation of Xu & Goddard" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 136 +citation = "X. Xu and W. A. Goddard III, J. Chem. Phys. 121, 4068 (2004)" +doi = "10.1063/1.1771632" +[GGA_C_XPBE.parameters] +beta = 0.089809 +gamma = 0.02043355766025040154 +B = 1.0 + +[GGA_C_APBE] +dens_threshold = 1e-12 +name = "APBE correlation: mu fixed from the semiclassical neutral atom" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 186 +citation = "L. A. Constantin, E. Fabiano, S. Laricchia, and F. Della Sala, Phys. Rev. Lett. 106, 186406 (2011)" +doi = "10.1103/PhysRevLett.106.186406" +[GGA_C_APBE.parameters] +beta = 0.079030523241023461723 +gamma = 0.031090690869654895034 +B = 1.0 + +[GGA_C_PBEFE] +dens_threshold = 1e-12 +name = "PBEfe correlation for formation energies" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 258 +citation = "R. Sarmiento-Perez, S. Botti, and M. A. L. Marques, J. Chem. Theory Comput. 11, 3844 (2015)" +doi = "10.1021/acs.jctc.5b00529" +[GGA_C_PBEFE.parameters] +beta = 0.043 +gamma = 0.031090690869654895034 +B = 1.0 + +[GGA_C_PBE_MOL] +dens_threshold = 1e-12 +name = "PBEmol correlation of Del Campo, Gazquez, Trickey & Vela" +maple_name = "gga_c_pbe" +family = "gga" +kind = "c" +kernel = "gga_c_pbe" +libxc_id = 272 +citation = "J. M. del Campo, J. L. Gazquez, S. B. Trickey, and A. Vela, J. Chem. Phys. 136, 104108 (2012)" +doi = "10.1063/1.3691197" +[GGA_C_PBE_MOL.parameters] +beta = 0.08384 +gamma = 0.031090690869654895034 +B = 1.0 + +[GGA_C_LYP] +dens_threshold = 1e-14 +name = "Lee, Yang & Parr correlation" +maple_name = "gga_c_lyp" +family = "gga" +kind = "c" +kernel = "gga_c_lyp" +libxc_id = 131 +citation = "C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988); B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Chem. Phys. Lett. 157, 200 (1989)" +doi = "10.1103/PhysRevB.37.785" +[GGA_C_LYP.parameters] +a = 0.04918 +b = 0.132 +c = 0.2533 +d = 0.349 diff --git a/src/funxc/data/lda.toml b/src/funxc/data/lda.toml new file mode 100644 index 0000000..d43ca9c --- /dev/null +++ b/src/funxc/data/lda.toml @@ -0,0 +1,109 @@ +# LDA functionals. Parameter values verified against libxc 7.0.0 src/lda_*.c. + +[LDA_X] +name = "Slater exchange" +maple_name = "lda_x" +family = "lda" +kind = "x" +kernel = "lda_x" +libxc_id = 1 +citation = "P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930); F. Bloch, Z. Phys. 57, 545 (1929)" +doi = "10.1017/S0305004100016108" +[LDA_X.parameters] +alpha = 1.0 + +[LDA_C_VWN] +name = "Vosko, Wilk & Nusair correlation (VWN5)" +maple_name = "lda_c_vwn" +family = "lda" +kind = "c" +kernel = "lda_c_vwn" +libxc_id = 7 +citation = "S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)" +doi = "10.1139/p80-159" + +[LDA_C_VWN_RPA] +name = "Vosko, Wilk & Nusair correlation (RPA parametrization)" +maple_name = "lda_c_vwn_rpa" +family = "lda" +kind = "c" +kernel = "lda_c_vwn_rpa" +libxc_id = 8 +citation = "S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)" +doi = "10.1139/p80-159" + +[LDA_C_PW] +name = "Perdew & Wang 1992 correlation" +maple_name = "lda_c_pw" +family = "lda" +kind = "c" +kernel = "lda_c_pw" +libxc_id = 12 +citation = "J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992)" +doi = "10.1103/PhysRevB.45.13244" +[LDA_C_PW.parameters] +pp = [1.0, 1.0, 1.0] +a = [0.031091, 0.015545, 0.016887] +alpha1 = [0.21370, 0.20548, 0.11125] +beta1 = [7.5957, 14.1189, 10.357] +beta2 = [3.5876, 6.1977, 3.6231] +beta3 = [1.6382, 3.3662, 0.88026] +beta4 = [0.49294, 0.62517, 0.49671] +fz20 = 1.709921 + +[LDA_C_PW_MOD] +name = "Perdew & Wang 1992 correlation, modified (higher-precision constants)" +maple_name = "lda_c_pw" +family = "lda" +kind = "c" +kernel = "lda_c_pw" +libxc_id = 13 +citation = "J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992)" +doi = "10.1103/PhysRevB.45.13244" +[LDA_C_PW_MOD.parameters] +pp = [1.0, 1.0, 1.0] +a = [0.0310907, 0.01554535, 0.0168869] +alpha1 = [0.21370, 0.20548, 0.11125] +beta1 = [7.5957, 14.1189, 10.357] +beta2 = [3.5876, 6.1977, 3.6231] +beta3 = [1.6382, 3.3662, 0.88026] +beta4 = [0.49294, 0.62517, 0.49671] +fz20 = 1.709920934161365617563962776245 + +[LDA_C_OB_PW] +name = "Ortiz & Ballone refit of PW92 correlation" +maple_name = "lda_c_pw" +family = "lda" +kind = "c" +kernel = "lda_c_pw" +libxc_id = 14 +citation = "G. Ortiz and P. Ballone, Phys. Rev. B 50, 1391 (1994)" +doi = "10.1103/PhysRevB.50.1391" +[LDA_C_OB_PW.parameters] +pp = [1.0, 1.0, 1.0] +a = [0.031091, 0.015545, 0.016887] +alpha1 = [0.026481, 0.022465, 0.11125] +beta1 = [7.5957, 14.1189, 10.357] +beta2 = [3.5876, 6.1977, 3.6231] +beta3 = [-0.46647, -0.56043, 0.88026] +beta4 = [0.13354, 0.11313, 0.49671] +fz20 = 1.709921 + +[LDA_C_PW_RPA] +name = "Perdew & Wang 1992 fit of the RPA" +maple_name = "lda_c_pw" +family = "lda" +kind = "c" +kernel = "lda_c_pw" +libxc_id = 25 +citation = "J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992)" +doi = "10.1103/PhysRevB.45.13244" +[LDA_C_PW_RPA.parameters] +pp = [0.75, 0.75, 1.0] +a = [0.031091, 0.015545, 0.016887] +alpha1 = [0.082477, 0.035374, 0.028829] +beta1 = [5.1486, 6.4869, 10.357] +beta2 = [1.6483, 1.3083, 3.6231] +beta3 = [0.23647, 0.15180, 0.47990] +beta4 = [0.20614, 0.082349, 0.12279] +fz20 = 1.709921 diff --git a/src/funxc/data/mgga.toml b/src/funxc/data/mgga.toml new file mode 100644 index 0000000..425e4ee --- /dev/null +++ b/src/funxc/data/mgga.toml @@ -0,0 +1,42 @@ +# meta-GGA functionals. Parameter values and flags verified against +# libxc 7.0.0 src/mgga_*.c (dens_threshold and XC_FLAGS_* in the +# xc_func_info structs; none of the LTA family sets ENFORCE_FHC). + +[MGGA_X_LTA] +name = "Local tau approximation" +maple_name = "mgga_x_lta" +family = "mgga" +kind = "x" +kernel = "mgga_x_lta" +libxc_id = 201 +needs_tau = true +citation = "M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 111, 911 (1999)" +doi = "10.1063/1.479374" +[MGGA_X_LTA.parameters] +ltafrac = 1.0 + +[MGGA_X_TLDA] +name = "LDA-type exchange with tau-dependent potential" +maple_name = "mgga_x_lta" +family = "mgga" +kind = "x" +kernel = "mgga_x_lta" +libxc_id = 685 +needs_tau = true +citation = "F. G. Eich and M. Hellgren, J. Chem. Phys. 141, 224107 (2014)" +doi = "10.1063/1.4903273" +[MGGA_X_TLDA.parameters] +ltafrac = 0.25 + +[MGGA_X_HLTA] +name = "Half-and-half meta-LDAized LDA exchange" +maple_name = "mgga_x_lta" +family = "mgga" +kind = "x" +kernel = "mgga_x_lta" +libxc_id = 698 +needs_tau = true +citation = "S. Lehtola and M. A. L. Marques, J. Chem. Theory Comput. 17, 943 (2021)" +doi = "10.1021/acs.jctc.0c01147" +[MGGA_X_HLTA.parameters] +ltafrac = 0.5 diff --git a/src/funxc/data/r2scan.toml b/src/funxc/data/r2scan.toml new file mode 100644 index 0000000..065a275 --- /dev/null +++ b/src/funxc/data/r2scan.toml @@ -0,0 +1,80 @@ +# r2SCAN meta-GGA family. Parameter values and flags verified against +# libxc 7.0.0 src/mgga_x_r2scan.c and src/mgga_c_r2scan.c (dens_threshold and +# XC_FLAGS_* in the xc_func_info structs). All four carry XC_FLAGS_NEEDS_TAU +# only among the mGGA input flags (no NEEDS_LAPLACIAN, no ENFORCE_FHC). +# +# Exchange params order in par_r2scan[] = {c1, c2, d, k1, eta, dp2}: +# r2scan = {0.667, 0.8, 1.24, 0.065, 0.001, 0.361} +# r2scan01 = {0.667, 0.8, 1.24, 0.065, 0.01, 0.361} (only eta differs) +# Correlation exposes only eta (dp2=0.361 hardcoded in the maple): +# r2scan = 0.001, r2scan01 = 0.01 +# +# LDA-correlation base for the correlation kernel: PW92 in the PW_MOD +# parametrization (lda_c_pw with lda_c_pw_params AND lda_c_pw_modified_params), +# i.e. gga_c._PW_MOD; the SCAN e0 branch uses its own eclda0 Pade form. + +[MGGA_X_R2SCAN] +name = "Re-regularized SCAN exchange by Furness et al" +maple_name = "mgga_x_r2scan" +family = "mgga" +kind = "x" +kernel = "mgga_x_r2scan" +libxc_id = 497 +dens_threshold = 1e-11 +needs_tau = true +citation = "J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020); erratum ibid. 11, 9248 (2020)" +doi = "10.1021/acs.jpclett.0c02405" +[MGGA_X_R2SCAN.parameters] +c1 = 0.667 +c2 = 0.8 +d = 1.24 +k1 = 0.065 +eta = 0.001 +dp2 = 0.361 + +[MGGA_X_R2SCAN01] +name = "Re-regularized SCAN exchange by Furness et al with larger value for eta" +maple_name = "mgga_x_r2scan" +family = "mgga" +kind = "x" +kernel = "mgga_x_r2scan" +libxc_id = 645 +dens_threshold = 1e-11 +needs_tau = true +citation = "J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020); erratum ibid. 11, 9248 (2020)" +doi = "10.1021/acs.jpclett.0c02405" +[MGGA_X_R2SCAN01.parameters] +c1 = 0.667 +c2 = 0.8 +d = 1.24 +k1 = 0.065 +eta = 0.01 +dp2 = 0.361 + +[MGGA_C_R2SCAN] +name = "Re-regularized SCAN correlation by Furness et al" +maple_name = "mgga_c_r2scan" +family = "mgga" +kind = "c" +kernel = "mgga_c_r2scan" +libxc_id = 498 +dens_threshold = 1e-15 +needs_tau = true +citation = "J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020); erratum ibid. 11, 9248 (2020)" +doi = "10.1021/acs.jpclett.0c02405" +[MGGA_C_R2SCAN.parameters] +eta = 0.001 + +[MGGA_C_R2SCAN01] +name = "Re-regularized SCAN correlation with larger value for eta" +maple_name = "mgga_c_r2scan" +family = "mgga" +kind = "c" +kernel = "mgga_c_r2scan" +libxc_id = 642 +dens_threshold = 1e-15 +needs_tau = true +citation = "J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020); erratum ibid. 11, 9248 (2020)" +doi = "10.1021/acs.jpclett.0c02405" +[MGGA_C_R2SCAN01.parameters] +eta = 0.01 diff --git a/src/funxc/data/scan.toml b/src/funxc/data/scan.toml new file mode 100644 index 0000000..1167648 --- /dev/null +++ b/src/funxc/data/scan.toml @@ -0,0 +1,82 @@ +# SCAN meta-GGA family. Parameter values and flags verified against +# libxc 7.0.0 src/mgga_x_scan.c and src/mgga_c_scan.c (xc_func_info structs): +# all carry XC_FLAGS_3D | XC_FLAGS_NEEDS_TAU (no ENFORCE_FHC, no lapl), +# dens_threshold 1e-15. Exchange params are the {c1, c2, d, k1} tables +# par_scan / par_revscan; the correlation f(alpha) constants (c1=0.64, +# c2=1.5, d=0.7) are hardcoded in the maple, not exposed as parameters. +# +# Correlation LDA base: the ec1 (PBE-like) limit uses the *modified* PW92 +# LDA correlation (lda_c_pw_modified_params, i.e. _PW_MOD) with the regtpss +# Hu-Langreth rs-dependent beta override; the ec0 (LDA-like) limit uses the +# SCAN eclda0 rational form. No VWN anywhere. + +[MGGA_X_SCAN] +name = "SCAN exchange of Sun, Ruzsinszky, and Perdew" +maple_name = "mgga_x_scan" +family = "mgga" +kind = "x" +kernel = "mgga_x_scan" +libxc_id = 263 +needs_tau = true +needs_lapl = false +enforce_fhc = false +dens_threshold = 1e-15 +citation = "J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015)" +doi = "10.1103/PhysRevLett.115.036402" +[MGGA_X_SCAN.parameters] +c1 = 0.667 +c2 = 0.8 +d = 1.24 +k1 = 0.065 + +[MGGA_X_REVSCAN] +name = "revised SCAN exchange" +maple_name = "mgga_x_scan" +family = "mgga" +kind = "x" +kernel = "mgga_x_scan" +libxc_id = 581 +needs_tau = true +needs_lapl = false +enforce_fhc = false +dens_threshold = 1e-15 +citation = "P. D. Mezei, G. I. Csonka, and M. Kállay, J. Chem. Theory Comput. 14, 2469 (2018)" +doi = "10.1021/acs.jctc.8b00072" +[MGGA_X_REVSCAN.parameters] +c1 = 0.607 +c2 = 0.7 +d = 1.37 +k1 = 0.065 + +[MGGA_C_SCAN] +name = "SCAN correlation of Sun, Ruzsinszky, and Perdew" +maple_name = "mgga_c_scan" +family = "mgga" +kind = "c" +kernel = "mgga_c_scan" +libxc_id = 267 +needs_tau = true +needs_lapl = false +enforce_fhc = false +dens_threshold = 1e-15 +citation = "J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015)" +doi = "10.1103/PhysRevLett.115.036402" +[MGGA_C_SCAN.parameters] + +[MGGA_C_SCAN_RVV10] +name = "SCAN + rVV10 correlation" +maple_name = "mgga_c_scan" +family = "mgga" +kind = "c" +kernel = "mgga_c_scan" +libxc_id = 292 +needs_tau = true +needs_lapl = false +enforce_fhc = false +dens_threshold = 1e-15 +# Semilocally identical to MGGA_C_SCAN; the rVV10 piece is a nonlocal +# correction (nlc_b=15.7, nlc_C=0.0093) applied outside the kernel and not +# exercised by the pointwise regression harness. +citation = "H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Phys. Rev. X 6, 041005 (2016)" +doi = "10.1103/PhysRevX.6.041005" +[MGGA_C_SCAN_RVV10.parameters] diff --git a/src/funxc/evaluate.py b/src/funxc/evaluate.py new file mode 100644 index 0000000..1ee575f --- /dev/null +++ b/src/funxc/evaluate.py @@ -0,0 +1,140 @@ +"""Evaluation core: (rho, sigma[, lapl, tau]) -> Maple internal variables -> zk, v*. + +Mirrors libxc's work_lda_inc.c / work_gga_inc.c / work_mgga_inc.c +preprocessing exactly: + +* points with total raw density < dens_threshold produce zero outputs; +* per-spin rho is floored at dens_threshold, sigma_uu/dd at + dens_threshold^(8/3), sigma_ud is clipped to +-(sigma_uu+sigma_dd)/2; +* for functionals that need tau, per-spin tau is floored at TAU_THRESHOLD, + and (only under the ENFORCE_FHC flag) sigma_ss is capped at 8*rho_s*tau_s + so the Fermi-hole curvature 1 - xs^2/(8 ts) stays positive; the lapl + channel is passed through unguarded ("lapl can have any values"); +* derivatives are taken with respect to the *floored* inputs (libxc never + chains through the floors), hence the straight-through floor_st/clip_st; +* 1 +- zeta are computed as 2*rho_sigma/n (cancellation-free, as in libxc's + generated code) and packed into a guards.Spin context together with the + exact per-channel screening masks. + +Derivatives are autodiff of the energy density e = n * eps, exactly the +quantities libxc calls vrho, vsigma, vlapl, vtau. Batches are vmap'd in api.py. + +Conventions (matching libxc): + polarized: rho = (rho_up, rho_dn), sigma = (sigma_uu, sigma_ud, sigma_dd), + lapl = (lapl_up, lapl_dn), tau = (tau_up, tau_dn) + unpolarized: rho, sigma = |grad rho|^2, lapl, tau are totals; each spin + channel sees rho/2, sigma/4, lapl/2, tau/2 (so + x_sigma = 2^(1/3)*x_t and t_sigma = 2^(2/3)*tau/rho^(5/3)). +tau is libxc's kinetic-energy density including the factor 1/2. +""" + +import jax.numpy as jnp + +from .constants import RS_FACTOR, TAU_THRESHOLD, TINY +from .guards import Spin, clip_st, floor_st, masked, safe_max + + +def make_eps(spec, polarized: bool): + """Per-point eps(rho[, sigma[, lapl, tau]]) — XC energy per particle (zk).""" + kfn = spec.kernel_fn + p = dict(spec.parameters) + thr = spec.dens_threshold + sigma_floor = (thr ** (4.0 / 3.0)) ** 2 + is_lda = spec.family == "lda" + is_mgga = spec.family == "mgga" + + if polarized: + + def eps(rho, sigma=None, lapl=None, tau=None): + dead = rho[0] + rho[1] < thr # total-density screen, raw inputs + ru = floor_st(rho[0], thr) + rd = floor_st(rho[1], thr) + n = ru + rd + rs = RS_FACTOR / n ** (1.0 / 3.0) + # z, 1+z, 1-z formed exactly as libxc's generated code does + # (t42 = (rho0-rho1)*(1/n); t43 = 1+t42; t50 = 1-t42). The + # subtraction 1-z loses relative accuracy near z = 1, but libxc + # is the parity oracle, so we reproduce its arithmetic including + # this cancellation rather than using the stabler 2*rho_dn/n. + inv_n = 1.0 / n + z = (ru - rd) * inv_n + sp = Spin( + z=z, + opz=1.0 + z, + omz=1.0 - z, + screen_up=ru <= thr, + screen_dn=rd <= thr, + ) + if is_lda: + return masked(dead, kfn(p, rs, sp)) + suu = floor_st(sigma[0], sigma_floor) + sdd = floor_st(sigma[2], sigma_floor) + if is_mgga and spec.needs_tau: + tu = floor_st(tau[0], TAU_THRESHOLD) + td = floor_st(tau[1], TAU_THRESHOLD) + if spec.enforce_fhc: + suu = clip_st(suu, None, 8.0 * ru * tu) + sdd = clip_st(sdd, None, 8.0 * rd * td) + s_ave = 0.5 * (suu + sdd) + sud = clip_st(sigma[1], -s_ave, s_ave) + st = safe_max(suu + 2.0 * sud + sdd, TINY**2) + xt = jnp.sqrt(st) / n ** (4.0 / 3.0) + xs0 = jnp.sqrt(suu) / ru ** (4.0 / 3.0) + xs1 = jnp.sqrt(sdd) / rd ** (4.0 / 3.0) + if not is_mgga: + return masked(dead, kfn(p, rs, sp, xt, xs0, xs1)) + ru53 = ru ** (5.0 / 3.0) + rd53 = rd ** (5.0 / 3.0) + ts0 = tu / ru53 if spec.needs_tau else 0.0 + ts1 = td / rd53 if spec.needs_tau else 0.0 + us0 = lapl[0] / ru53 if spec.needs_lapl else 0.0 + us1 = lapl[1] / rd53 if spec.needs_lapl else 0.0 + return masked(dead, kfn(p, rs, sp, xt, xs0, xs1, us0, us1, ts0, ts1)) + + else: + + def eps(rho, sigma=None, lapl=None, tau=None): + dead = rho < thr + n = floor_st(rho, thr) + rs = RS_FACTOR / n ** (1.0 / 3.0) + screen = n / 2.0 <= thr # generated code screens on rho/2 + sp = Spin(z=0.0, opz=1.0, omz=1.0, screen_up=screen, screen_dn=screen) + if is_lda: + return masked(dead, kfn(p, rs, sp)) + st = floor_st(sigma, sigma_floor) + if is_mgga and spec.needs_tau: + tt = floor_st(tau, TAU_THRESHOLD) + if spec.enforce_fhc: + st = clip_st(st, None, 8.0 * n * tt) + xt = jnp.sqrt(st) / n ** (4.0 / 3.0) + xs = 2.0 ** (1.0 / 3.0) * xt + if not is_mgga: + return masked(dead, kfn(p, rs, sp, xt, xs, xs)) + # each channel sees tau/2, lapl/2, rho/2: 2^(2/3)*total/n^(5/3) + n53 = n ** (5.0 / 3.0) + ts = 2.0 ** (2.0 / 3.0) * tt / n53 if spec.needs_tau else 0.0 + us = 2.0 ** (2.0 / 3.0) * lapl / n53 if spec.needs_lapl else 0.0 + return masked(dead, kfn(p, rs, sp, xt, xs, xs, us, us, ts, ts)) + + return eps + + +def make_energy_density(spec, polarized: bool): + """Per-point e(rho[, sigma[, lapl, tau]]) = n * eps; its partial + derivatives are libxc's vrho, vsigma, vlapl, vtau. Uses the floored n, + as libxc does.""" + eps = make_eps(spec, polarized) + thr = spec.dens_threshold + + if polarized: + + def e(rho, sigma=None, lapl=None, tau=None): + n = floor_st(rho[0], thr) + floor_st(rho[1], thr) + return n * eps(rho, sigma, lapl, tau) + + else: + + def e(rho, sigma=None, lapl=None, tau=None): + return floor_st(rho, thr) * eps(rho, sigma, lapl, tau) + + return e diff --git a/src/funxc/guards.py b/src/funxc/guards.py new file mode 100644 index 0000000..509c8a0 --- /dev/null +++ b/src/funxc/guards.py @@ -0,0 +1,121 @@ +"""NaN-safe primitives and libxc-parity helpers shared by all kernels. + +Numerical conventions found in libxc (and mirrored here): + +* Input preprocessing lives in libxc's work_lda_inc.c / work_gga_inc.c, NOT in + the kernels: per-spin rho is floored at dens_threshold, sigma_uu/dd at + sigma_threshold^2 = dens_threshold^(8/3), sigma_ud is clipped to + +-(sigma_uu+sigma_dd)/2, and points with *total* raw density below + dens_threshold are skipped (all outputs zero). Derivatives are reported + with respect to the floored inputs — a straight-through floor in JAX terms. + +* The generated code never forms 1 +- zeta by subtraction (catastrophic + cancellation near zeta = +-1); it uses 2*rho_sigma/n directly. Kernels here + therefore receive a `Spin` context carrying stably computed opz = 1+zeta, + omz = 1-zeta and exact per-channel screening masks, instead of deriving + them from zeta. + +* Guards follow the double-where idiom: inputs are masked to safe values + *before* the dangerous op and results masked after, so no derivative order + sees inf/NaN through the dead branch. +""" + +from typing import NamedTuple + +import jax +import jax.numpy as jnp + +from .constants import RS_FACTOR, X_FACTOR_C, ZETA_THRESHOLD + + +class Spin(NamedTuple): + """Spin-polarization context, computed cancellation-free in evaluate.py. + + z: zeta = (rho_up - rho_dn)/n + opz: 1 + zeta = 2*rho_up/n + omz: 1 - zeta = 2*rho_dn/n + screen_up/screen_dn: rho_sigma <= dens_threshold (exact, from raw inputs) + """ + + z: jnp.ndarray + opz: jnp.ndarray + omz: jnp.ndarray + screen_up: jnp.ndarray + screen_dn: jnp.ndarray + + def flipped(self): + return Spin(-self.z, self.omz, self.opz, self.screen_dn, self.screen_up) + + +def safe_max(x, floor): + return jnp.maximum(x, floor) + + +def masked(cond, value): + """Zero out `value` where `cond` holds. `value` must already be finite.""" + return jnp.where(cond, 0.0, value) + + +def floor_st(x, floor): + """max(x, floor) in value, identity in gradient (straight-through). + + Matches libxc, which evaluates kernels at floored inputs but reports + derivatives with respect to them without chaining through the max. + """ + return x + jax.lax.stop_gradient(jnp.maximum(x, floor) - x) + + +def clip_st(x, lo, hi): + """clip in value, identity in gradient (straight-through).""" + return x + jax.lax.stop_gradient(jnp.clip(x, lo, hi) - x) + + +def safe_exp(x): + """exp with the argument clipped to avoid overflow to inf.""" + return jnp.exp(jnp.clip(x, -700.0, 700.0)) + + +def opz_pow_n(opz, n): + """(1+z)^n with libxc's zeta-threshold floor (util.mpl opz_pow_n). + + `opz` is the stably computed 1+z (pass Spin.omz for (1-z)^n). + """ + safe = safe_max(opz, ZETA_THRESHOLD) + return jnp.where(opz <= ZETA_THRESHOLD, ZETA_THRESHOLD**n, safe**n) + + +def opz_thr(opz, omz): + """1 + z_thr(z) (util.mpl z_thr in the 1+z representation): + clamps 1+z into [ZETA_THRESHOLD, 2 - ZETA_THRESHOLD].""" + out = jnp.where(opz <= ZETA_THRESHOLD, ZETA_THRESHOLD, opz) + return jnp.where(omz <= ZETA_THRESHOLD, 2.0 - ZETA_THRESHOLD, out) + + +def f_zeta(sp: Spin): + """Spin interpolation f(zeta) of Perdew-Wang Eq. (9).""" + return (opz_pow_n(sp.opz, 4.0 / 3.0) + opz_pow_n(sp.omz, 4.0 / 3.0) - 2.0) / ( + 2.0 ** (4.0 / 3.0) - 2.0 + ) + + +def mphi(sp: Spin): + """phi(zeta) = ((1+z)^(2/3) + (1-z)^(2/3))/2, used in PBE-type correlation.""" + return (opz_pow_n(sp.opz, 2.0 / 3.0) + opz_pow_n(sp.omz, 2.0 / 3.0)) / 2.0 + + +def lda_x_spin(rs, opz): + """Per-spin LDA exchange energy density (util.mpl lda_x_spin), 3D. + + `opz` is the (possibly opz_thr-clamped) 1+z of the channel. + """ + return ( + -X_FACTOR_C + * opz_pow_n(opz, 4.0 / 3.0) + * 2.0 ** (-4.0 / 3.0) + * (RS_FACTOR / rs) + ) + + +def tt(rs, sp: Spin, xt): + """PBE correlation reduced gradient t (util.mpl tt).""" + return xt / (4.0 * 2.0 ** (1.0 / 3.0) * mphi(sp) * jnp.sqrt(rs)) diff --git a/src/funxc/kernels/__init__.py b/src/funxc/kernels/__init__.py new file mode 100644 index 0000000..45a92e2 --- /dev/null +++ b/src/funxc/kernels/__init__.py @@ -0,0 +1,34 @@ +"""Math-kernel registry. + +A kernel is a pure JAX function of libxc's Maple internal variables: + + LDA: eps = kernel(p, rs, sp) + GGA: eps = kernel(p, rs, sp, xt, xs0, xs1) + mGGA: eps = kernel(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1) + +where ``p`` is the parameter dict from the functional database, ``rs`` the +Wigner-Seitz radius, ``sp`` the guards.Spin context, ``xt``/``xs*`` the total +and per-spin reduced gradients x = |grad rho|/rho^(4/3), ``u*`` and ``t*`` +the per-spin lapl/rho^(5/3) and tau/rho^(5/3). Kernels return the XC energy +per particle (libxc's ``zk``); all derivatives come from autodiff. + +Every parameter-only functional variant maps to one of these kernels through +the TOML database — new variants should not need new code here. +""" + +KERNELS: dict[str, callable] = {} + + +def kernel(name): + """Register a math kernel under `name` (the database `kernel` key).""" + + def deco(fn): + if name in KERNELS: + raise ValueError(f"duplicate kernel name: {name}") + KERNELS[name] = fn + return fn + + return deco + + +from . import b97, gga_c, gga_x, lda, mgga_x, r2scan, scan, tpss # noqa: E402,F401 (populate KERNELS) diff --git a/src/funxc/kernels/b97.py b/src/funxc/kernels/b97.py new file mode 100644 index 0000000..099fdb3 --- /dev/null +++ b/src/funxc/kernels/b97.py @@ -0,0 +1,93 @@ +"""B97 power-series kernel (maple/b97.mpl + gga_exc/gga_xc_b97.mpl). + +The B97 form weights Stoll-partitioned LDA energies with an inhomogeneity +power series g(gamma, cc, x) = sum_i cc[i] * u^i, u = gamma*x^2/(1+gamma*x^2): + + exc = e_x^par(up)*g_x(xs0) + e_x^par(dn)*g_x(xs1) + + e_c^par(up)*g_ss(xs0) + e_c^par(dn)*g_ss(xs1) + + e_c^perp * g_ab(sqrt((xs0^2+xs1^2)/2)) + +where the Stoll partition of an LDA is + e^par(z) = (1+z)/2 * lda(rs * 2^(1/3) * (1+z)^(-1/3), zeta=1) + e^perp = lda(rs, z) - e^par(z) - e^par(-z) +(util.mpl lda_stoll_par / lda_stoll_perp). Exchange uses Slater LDA and has +no perpendicular term; correlation uses PW92 with the *original* constants +(gga_xc_b97.mpl defines lda_c_pw_params, not the modified set PBE uses). + +This one kernel covers B97/B97-1/2/3, B97-D/3c/GGA1, the HCTH family and +HLE16 — they differ only in the c_x/c_ss/c_ab coefficient tables. +""" + +from ..constants import ZETA_THRESHOLD +from ..guards import Spin, lda_x_spin, masked, opz_pow_n +from . import kernel +from .lda import pw92_f + +# Series gammas are part of the functional form (hard-coded in b97.mpl). +GAMMA_X = 0.004 +GAMMA_SS = 0.2 +GAMMA_AB = 0.006 + +CBRT2 = 2.0 ** (1.0 / 3.0) + +# Original PW92 constants (lda_c_pw.mpl, lda_c_pw_params branch). +_PW_ORIG = { + "pp": (1.0, 1.0, 1.0), + "a": (0.031091, 0.015545, 0.016887), + "alpha1": (0.21370, 0.20548, 0.11125), + "beta1": (7.5957, 14.1189, 10.357), + "beta2": (3.5876, 6.1977, 3.6231), + "beta3": (1.6382, 3.3662, 0.88026), + "beta4": (0.49294, 0.62517, 0.49671), + "fz20": 1.709921, +} + +# Fully polarized spin context for evaluating an LDA at zeta = 1; the +# (1-z)^n factors inside hit the ZETA_THRESHOLD clamp exactly as libxc's +# generated code does. +_SP_FULL = Spin(z=1.0, opz=2.0, omz=0.0, screen_up=False, screen_dn=False) + + +def b97_g(gamma, cc, x2): + """The B97 series g as a function of x^2 (u saturates to 1 as x -> inf).""" + u = gamma * x2 / (1.0 + gamma * x2) + u2 = u * u + return cc[0] + cc[1] * u + cc[2] * u2 + cc[3] * u2 * u + cc[4] * u2 * u2 + + +def _stoll_par(lda_at_full, rs, opz, screen_channel): + """util.mpl lda_stoll_par: (1+z)/2 * lda(rs*2^(1/3)*(1+z)^(-1/3), 1), + screened where the channel density or 1+z is below threshold.""" + rs_s = rs * CBRT2 * opz_pow_n(opz, -1.0 / 3.0) + screen = screen_channel | (opz <= ZETA_THRESHOLD) + return masked(screen, 0.5 * opz * lda_at_full(rs_s)) + + +@kernel("gga_xc_b97") +def gga_xc_b97(p, rs, sp: Spin, xt, xs0, xs1): + cx, css, cab = p["c_x"], p["c_ss"], p["c_ab"] + xs0_2 = xs0 * xs0 + xs1_2 = xs1 * xs1 + + # Exchange: per-channel Stoll LDA (the perpendicular part is zero). + def slater_full(r): + return lda_x_spin(r, 2.0) + + ex_up = _stoll_par(slater_full, rs, sp.opz, sp.screen_up) + ex_dn = _stoll_par(slater_full, rs, sp.omz, sp.screen_dn) + ex = ex_up * b97_g(GAMMA_X, cx, xs0_2) + ex_dn * b97_g(GAMMA_X, cx, xs1_2) + + # Correlation: PW92 split into parallel and perpendicular components. + def pw_full(r): + return pw92_f(_PW_ORIG, r, _SP_FULL) + + ec_par_up = _stoll_par(pw_full, rs, sp.opz, sp.screen_up) + ec_par_dn = _stoll_par(pw_full, rs, sp.omz, sp.screen_dn) + ec_perp = pw92_f(_PW_ORIG, rs, sp) - ec_par_up - ec_par_dn + + ec = ( + ec_par_up * b97_g(GAMMA_SS, css, xs0_2) + + ec_par_dn * b97_g(GAMMA_SS, css, xs1_2) + + ec_perp * b97_g(GAMMA_AB, cab, 0.5 * (xs0_2 + xs1_2)) + ) + return ex + ec diff --git a/src/funxc/kernels/gga_c.py b/src/funxc/kernels/gga_c.py new file mode 100644 index 0000000..a37f7d4 --- /dev/null +++ b/src/funxc/kernels/gga_c.py @@ -0,0 +1,82 @@ +"""GGA correlation kernels: PBE H(t) gradient correction and LYP. + +Forms follow libxc's maple/gga_exc/gga_c_pbe.mpl and gga_c_lyp.mpl verbatim. +""" + +import math + +import jax.numpy as jnp + +from ..constants import RS_FACTOR, TINY +from ..guards import Spin, mphi, opz_pow_n, safe_max, tt +from . import kernel +from .lda import pw92_f + +# gga_c_pbe.mpl includes lda_c_pw.mpl with *both* lda_c_pw_params and +# lda_c_pw_modified_params defined, i.e. the PW_MOD parametrization. +_PW_MOD = { + "pp": (1.0, 1.0, 1.0), + "a": (0.0310907, 0.01554535, 0.0168869), + "alpha1": (0.21370, 0.20548, 0.11125), + "beta1": (7.5957, 14.1189, 10.357), + "beta2": (3.5876, 6.1977, 3.6231), + "beta3": (1.6382, 3.3662, 0.88026), + "beta4": (0.49294, 0.62517, 0.49671), + "fz20": 1.709920934161365617563962776245, +} + + +@kernel("gga_c_pbe") +def gga_c_pbe(p, rs, sp: Spin, xt, xs0, xs1): + """PBE correlation: eps_c^PW92 + H(rs, z, t) (Eqs. 7-8 of PBE).""" + beta, gamma, bb = p["beta"], p["gamma"], p["B"] + phi = mphi(sp) + ec = pw92_f(_PW_MOD, rs, sp) + t = tt(rs, sp, xt) + + # A = beta/gamma / (exp(-ec/(gamma*phi^3)) - 1); ec <= 0 so expm1 >= 0. + # The denominator vanishes only in the ec -> 0- limit (rs -> inf), where + # H -> 0; the floor keeps that branch finite for autodiff. + denom = safe_max(jnp.expm1(-ec / (gamma * phi**3)), TINY) + a = beta / gamma / denom + + f1 = t**2 + bb * a * t**4 + f2 = beta * f1 / (gamma * (1.0 + a * f1)) + h = gamma * phi**3 * jnp.log1p(f2) + return ec + h + + +LYP_CF = 3.0 / 10.0 * (3.0 * math.pi**2) ** (2.0 / 3.0) +LYP_AUX6 = 1.0 / 2.0 ** (8.0 / 3.0) +LYP_AUX4 = LYP_AUX6 / 4.0 +LYP_AUX5 = LYP_AUX4 / 18.0 + + +@kernel("gga_c_lyp") +def gga_c_lyp(p, rs, sp: Spin, xt, xs0, xs1): + """Lee-Yang-Parr correlation, Miehlich form (maple/gga_exc/gga_c_lyp.mpl).""" + a, b, c, d = p["a"], p["b"], p["c"], p["d"] + rr = rs / RS_FACTOR # = rho^(-1/3) + + omega = b * jnp.exp(-c * rr) / (1.0 + d * rr) + delta = (c + d / (1.0 + d * rr)) * rr + + # 1 - z^2 = (1+z)(1-z), computed cancellation-free + omz2 = sp.opz * sp.omz + opz83 = opz_pow_n(sp.opz, 8.0 / 3.0) + omz83 = opz_pow_n(sp.omz, 8.0 / 3.0) + + t1 = -omz2 / (1.0 + d * rr) + t2 = -(xt**2) * (omz2 * (47.0 - 7.0 * delta) / 72.0 - 2.0 / 3.0) + t3 = -LYP_CF / 2.0 * omz2 * (opz83 + omz83) + t4 = LYP_AUX4 * omz2 * (5.0 / 2.0 - delta / 18.0) * (xs0**2 * opz83 + xs1**2 * omz83) + t5 = LYP_AUX5 * omz2 * (delta - 11.0) * ( + xs0**2 * opz_pow_n(sp.opz, 11.0 / 3.0) + xs1**2 * opz_pow_n(sp.omz, 11.0 / 3.0) + ) + t6 = -LYP_AUX6 * ( + 2.0 / 3.0 * (xs0**2 * opz83 + xs1**2 * omz83) + - opz_pow_n(sp.opz, 2.0) * xs1**2 * omz83 / 4.0 + - opz_pow_n(sp.omz, 2.0) * xs0**2 * opz83 / 4.0 + ) + + return a * (t1 + omega * (t2 + t3 + t4 + t5 + t6)) diff --git a/src/funxc/kernels/gga_x.py b/src/funxc/kernels/gga_x.py new file mode 100644 index 0000000..c18668e --- /dev/null +++ b/src/funxc/kernels/gga_x.py @@ -0,0 +1,68 @@ +"""GGA exchange kernels: PBE rational shell, RPBE exponential, B88 arcsinh. + +Forms follow libxc's maple/gga_exc/*.mpl verbatim. Each enhancement factor +F(x) is combined with LDA exchange through the `gga_exchange` spin-scaling +combinator of maple/util.mpl. Kernel signature: +eps = f(p, rs, sp, xt, xs0, xs1) with sp a guards.Spin context. +""" + +import jax.numpy as jnp + +from ..constants import X2S, X_FACTOR_C +from ..guards import Spin, lda_x_spin, masked, opz_thr +from . import kernel + + +def gga_exchange(enhancement, p, rs, sp: Spin, xs0, xs1): + """util.mpl gga_exchange: per-spin eps_x^LDA * F(x_sigma). + + Each channel is screened on its density and its zeta is clamped by + z_thr (unlike plain LDA exchange, which clamps only through opz_pow_n). + """ + up = masked( + sp.screen_up, + lda_x_spin(rs, opz_thr(sp.opz, sp.omz)) * enhancement(p, xs0), + ) + dn = masked( + sp.screen_dn, + lda_x_spin(rs, opz_thr(sp.omz, sp.opz)) * enhancement(p, xs1), + ) + return up + dn + + +def pbe_enhancement(p, x): + """F(s) = 1 + kappa*(1 - kappa/(kappa + mu*s^2)) — saturates as s -> inf.""" + s = X2S * x + kappa, mu = p["kappa"], p["mu"] + return 1.0 + kappa * (1.0 - kappa / (kappa + mu * s**2)) + + +@kernel("gga_x_pbe") +def gga_x_pbe(p, rs, sp, xt, xs0, xs1): + """PBE-type rational exchange (maple/gga_exc/gga_x_pbe.mpl).""" + return gga_exchange(pbe_enhancement, p, rs, sp, xs0, xs1) + + +def rpbe_enhancement(p, x): + """F(s) = 1 + kappa*(1 - exp(-mu*s^2/kappa)); exponent <= 0, overflow-free.""" + s = X2S * x + kappa, mu = p["kappa"], p["mu"] + return 1.0 + kappa * (1.0 - jnp.exp(-mu * s**2 / kappa)) + + +@kernel("gga_x_rpbe") +def gga_x_rpbe(p, rs, sp, xt, xs0, xs1): + """RPBE exponential exchange (maple/gga_exc/gga_x_rpbe.mpl).""" + return gga_exchange(rpbe_enhancement, p, rs, sp, xs0, xs1) + + +def b88_enhancement(p, x): + """F(x) = 1 + (beta/Cx) x^2 / (1 + gamma*beta*x*arcsinh(x)).""" + beta, gamma = p["beta"], p["gamma"] + return 1.0 + beta / X_FACTOR_C * x**2 / (1.0 + gamma * beta * x * jnp.arcsinh(x)) + + +@kernel("gga_x_b88") +def gga_x_b88(p, rs, sp, xt, xs0, xs1): + """Becke 88 exchange (maple/gga_exc/gga_x_b88.mpl).""" + return gga_exchange(b88_enhancement, p, rs, sp, xs0, xs1) diff --git a/src/funxc/kernels/lda.py b/src/funxc/kernels/lda.py new file mode 100644 index 0000000..7d0f955 --- /dev/null +++ b/src/funxc/kernels/lda.py @@ -0,0 +1,108 @@ +"""LDA kernels: Dirac-Slater exchange, PW92 and VWN local correlation. + +Forms follow libxc's maple/lda_exc/*.mpl verbatim. Kernel signature: +eps = f(p, rs, sp) with sp a guards.Spin context. +""" + +import math + +import jax.numpy as jnp + +from ..guards import Spin, f_zeta, lda_x_spin, masked +from . import kernel + + +@kernel("lda_x") +def lda_x(p, rs, sp: Spin): + """Dirac-Slater exchange (maple/lda_exc/lda_x.mpl). + + No zeta clamp here (the maple uses raw z); screening is per channel. + """ + alpha = p["alpha"] + up = masked(sp.screen_up, lda_x_spin(rs, sp.opz)) + dn = masked(sp.screen_dn, lda_x_spin(rs, sp.omz)) + return alpha * (up + dn) + + +def pw92_g(a, alpha1, beta1, beta2, beta3, beta4, pp, rs): + """PW92 Eq. (10). Parametrizes -alpha_c in the spin-stiffness channel.""" + g_aux = ( + beta1 * jnp.sqrt(rs) + + beta2 * rs + + beta3 * rs**1.5 + + beta4 * rs ** (pp + 1.0) + ) + return -2.0 * a * (1.0 + alpha1 * rs) * jnp.log1p(1.0 / (2.0 * a * g_aux)) + + +def pw92_f(p, rs, sp: Spin): + """PW92 Eq. (8) spin interpolation over the three channels (ec0, ec1, -alpha_c).""" + g = [ + pw92_g( + p["a"][k], p["alpha1"][k], p["beta1"][k], p["beta2"][k], + p["beta3"][k], p["beta4"][k], p["pp"][k], rs, + ) + for k in range(3) + ] + fz = f_zeta(sp) + z4 = sp.z**4 + return g[0] + z4 * fz * (g[1] - g[0] + g[2] / p["fz20"]) - fz * g[2] / p["fz20"] + + +@kernel("lda_c_pw") +def lda_c_pw(p, rs, sp: Spin): + """Perdew-Wang 1992 local correlation (maple/lda_exc/lda_c_pw.mpl).""" + return pw92_f(p, rs, sp) + + +# --- VWN --------------------------------------------------------------- +# Constants from maple/vwn.mpl (original Rydberg values halved to Hartree). +# These define the kernel's structure, not database parameters: libxc +# exposes no external parameters for the VWN family either. + +VWN_A = (0.0310907, 0.01554535, -1.0 / (6.0 * math.pi**2)) +VWN_B = (3.72744, 7.06042, 1.13107) +VWN_C = (12.9352, 18.0578, 13.0045) +VWN_X0 = (-0.10498, -0.32500, -0.0047584) + +RPA_A = (0.0310907, 0.01554535, -1.0 / (6.0 * math.pi**2)) +RPA_B = (13.0720, 20.1231, 1.06835) +RPA_C = (42.7198, 101.578, 11.4813) +RPA_X0 = (-0.409286, -0.743294, -0.228344) + +FPP_VWN = 4.0 / (9.0 * (2.0 ** (1.0 / 3.0) - 1.0)) + + +def vwn_f_aux(a, b, c, x0, rs): + """The G(sqrt(rs)) Pade/log/arctan form of maple/vwn.mpl f_aux.""" + q = math.sqrt(4.0 * c - b * b) + f1 = 2.0 * b / q + f2 = b * x0 / (x0 * x0 + b * x0 + c) + f3 = 2.0 * (2.0 * x0 + b) / q + sqrt_rs = jnp.sqrt(rs) + fx = rs + b * sqrt_rs + c + return a * ( + jnp.log(rs / fx) + + (f1 - f2 * f3) * jnp.arctan(q / (2.0 * sqrt_rs + b)) + - f2 * jnp.log((sqrt_rs - x0) ** 2 / fx) + ) + + +@kernel("lda_c_vwn") +def lda_c_vwn(p, rs, sp: Spin): + """VWN5 correlation (maple/lda_exc/lda_c_vwn.mpl).""" + ec_p = vwn_f_aux(VWN_A[0], VWN_B[0], VWN_C[0], VWN_X0[0], rs) + ec_f = vwn_f_aux(VWN_A[1], VWN_B[1], VWN_C[1], VWN_X0[1], rs) + alpha = vwn_f_aux(VWN_A[2], VWN_B[2], VWN_C[2], VWN_X0[2], rs) + fz = f_zeta(sp) + z4 = sp.z**4 + return ec_p + alpha * fz * (1.0 - z4) / FPP_VWN + (ec_f - ec_p) * fz * z4 + + +@kernel("lda_c_vwn_rpa") +def lda_c_vwn_rpa(p, rs, sp: Spin): + """VWN RPA correlation (maple/lda_exc/lda_c_vwn_rpa.mpl).""" + ec_p = vwn_f_aux(RPA_A[0], RPA_B[0], RPA_C[0], RPA_X0[0], rs) + ec_f = vwn_f_aux(RPA_A[1], RPA_B[1], RPA_C[1], RPA_X0[1], rs) + fz = f_zeta(sp) + return ec_p * (1.0 - fz) + ec_f * fz diff --git a/src/funxc/kernels/mgga_x.py b/src/funxc/kernels/mgga_x.py new file mode 100644 index 0000000..b59eebf --- /dev/null +++ b/src/funxc/kernels/mgga_x.py @@ -0,0 +1,42 @@ +"""meta-GGA exchange kernels. + +Kernel signature (mirroring the Maple mgga_exc convention): + + eps = f(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1) + +with u_s = lapl_s/rho_s^(5/3) and t_s = tau_s/rho_s^(5/3) (tau includes the +factor 1/2). Enhancement factors F(x, u, t) are combined with LDA exchange +through the `mgga_exchange` spin-scaling combinator of maple/util.mpl. +""" + +from ..constants import K_FACTOR_C +from ..guards import Spin, lda_x_spin, masked, opz_thr +from . import kernel + + +def mgga_exchange(enhancement, p, rs, sp: Spin, xs0, xs1, u0, u1, t0, t1): + """util.mpl mgga_exchange: per-spin eps_x^LDA * F(x_s, u_s, t_s).""" + up = masked( + sp.screen_up, + lda_x_spin(rs, opz_thr(sp.opz, sp.omz)) * enhancement(p, xs0, u0, t0), + ) + dn = masked( + sp.screen_dn, + lda_x_spin(rs, opz_thr(sp.omz, sp.opz)) * enhancement(p, xs1, u1, t1), + ) + return up + dn + + +def lta_enhancement(p, x, u, t): + """F = (t/C_TF)^(4*ltafrac/5) — pure tau dependence (Ernzerhof-Scuseria). + + t is floored at TAU_THRESHOLD/rho^(5/3) upstream, so the fractional + power and its derivatives stay finite. + """ + return (t / K_FACTOR_C) ** (0.8 * p["ltafrac"]) + + +@kernel("mgga_x_lta") +def mgga_x_lta(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1): + """Local tau approximation (maple/mgga_exc/mgga_x_lta.mpl).""" + return mgga_exchange(lta_enhancement, p, rs, sp, xs0, xs1, u0, u1, t0, t1) diff --git a/src/funxc/kernels/r2scan.py b/src/funxc/kernels/r2scan.py new file mode 100644 index 0000000..9186ded --- /dev/null +++ b/src/funxc/kernels/r2scan.py @@ -0,0 +1,233 @@ +"""r2SCAN meta-GGA exchange and correlation kernels. + +Ground truth: libxc/maple/mgga_exc/mgga_x_r2scan.mpl and mgga_c_r2scan.mpl +(which pull in mgga_x_scan / mgga_x_rscan and mgga_c_scan / gga_c_scan_e0). + +r2SCAN (Furness, Kaplan, Ning, Perdew, Sun 2020 + erratum) re-regularizes +SCAN: a smooth regularized alpha-tilde, + + alpha_tilde = (t - x^2/8) / (K_FACTOR_C + eta*x^2/8) (exchange) + alpha_tilde = (tau_tot - x_t^2/8) / (K*tau_tot(z,1,1) + eta*x_t^2/8) (corr) + +and a *piecewise* interpolation f(alpha) replacing SCAN's singular exp switch: +exp(-c1*a/(1-a)) for a<=0, a degree-7 polynomial for 02.5. The three branches are evaluated on clamped +arguments (double-where idiom) so values and first/second derivatives are all +finite at the branch boundaries, at a=0, and as a->+-inf. + +The gradient-expansion coefficient in exchange gets a Gaussian damping +exp(-p^2/dp2^4) (dp2 = 0.361, the published erratum value, matching what +libxc 7.0.0 ships); correlation carries the analogous exp(-s^4/dp2^4) damping +in its S34 gradient correction term. +""" + +import math + +import jax +import jax.numpy as jnp + +from ..constants import K_FACTOR_C, MU_GE, X2S +from ..guards import f_zeta, mphi, opz_pow_n, safe_max, tt +from . import kernel +from .gga_c import _PW_MOD +from .lda import pw92_f +from .mgga_x import mgga_exchange + +# Floor for masked-out branches; never affects unmasked values. +TINY = 1e-100 + +TWO13 = 2.0 ** (1.0 / 3.0) + + +def _poly(coeffs, a): + """Horner evaluation of sum_i coeffs[i]*a^i.""" + out = coeffs[-1] + for c in reversed(coeffs[:-1]): + out = out * a + c + return out + + +def _r2scan_f_alpha(a, coeffs, c1, c2, d): + """my_piecewise5 interpolation f(alpha), NaN-safe via clamped inputs. + + Each branch uses an argument clamped into its own safe domain, so every + sub-expression (and its derivatives) is finite everywhere before the + where-select picks the active branch: + a <= 0 -> exp(-c1*a/(1-a)) (m_min(a,0): 1-a >= 1) + 0 < a<=2.5 -> polynomial(m_min(a,2.5)) + a > 2.5 -> -d*exp(c2/(1-a)) (m_max(a,2.5): 1-a <= -1.5) + """ + a_neg = jnp.minimum(a, 0.0) + fn_neg = jnp.exp(-c1 * a_neg / (1.0 - a_neg)) + a_sml = jnp.minimum(a, 2.5) + fn_sml = _poly(coeffs, a_sml) + a_lrg = jnp.maximum(a, 2.5) + fn_lrg = -d * jnp.exp(c2 / (1.0 - a_lrg)) + return jnp.where(a <= 0.0, fn_neg, jnp.where(a <= 2.5, fn_sml, fn_lrg)) + + +# ============================================================================ +# Exchange (maple/mgga_exc/mgga_x_r2scan.mpl) +# ============================================================================ + +SCAN_H0X = 1.174 +SCAN_A1 = 4.9479 + +# rSCAN switching polynomial: c[i] = coeff of alpha^i = maple rscan_fx[8-i]. +_RSCAN_FX = ( + 1.0, -0.667, -0.4445555, -0.663086601049, + 1.451297044490, -0.887998041597, 0.234528941479, -0.023185843322, +) +# C2 = -sum_{i=1..8} i*ff[9-i] * (1 - h0x) = -sum_{j=0..7}(j+1)*c[j] * (1-h0x). +_C2X = -sum((j + 1) * _RSCAN_FX[j] for j in range(8)) * (1.0 - SCAN_H0X) + + +def _scan_h1x(y, k1): + return 1.0 + k1 * (1.0 - k1 / (k1 + y)) + + +def _scan_gx(x): + """1 - exp(-a1/sqrt(X2S*x)); -> 1 (and finite gradient) as x -> 0.""" + inner = X2S * x + safe = jnp.where(inner > 0.0, inner, 1.0) + gx = 1.0 - jnp.exp(-SCAN_A1 / jnp.sqrt(safe)) + return jnp.where(inner > 0.0, gx, 1.0) + + +def _r2scan_x_enhancement(p, x, u, t): + """r2SCAN exchange enhancement F(x, t) (u = lapl is unused).""" + c1, c2, d, k1 = p["c1"], p["c2"], p["d"], p["k1"] + eta, dp2 = p["eta"], p["dp2"] + + pp = X2S**2 * x**2 # scan_p(x) + cn = 20.0 / 27.0 + eta * 5.0 / 3.0 + r2x = (cn * _C2X * jnp.exp(-(pp**2) / dp2**4) + MU_GE) * pp + h1 = _scan_h1x(r2x, k1) + + alpha = (t - x**2 / 8.0) / (K_FACTOR_C + eta * x**2 / 8.0) + fa = _r2scan_f_alpha(alpha, _RSCAN_FX, c1, c2, d) + + return (h1 + fa * (SCAN_H0X - h1)) * _scan_gx(x) + + +@kernel("mgga_x_r2scan") +def mgga_x_r2scan(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1): + """Re-regularized SCAN exchange (Furness et al 2020).""" + return mgga_exchange(_r2scan_x_enhancement, p, rs, sp, xs0, xs1, u0, u1, t0, t1) + + +# ============================================================================ +# Correlation (maple/mgga_exc/mgga_c_r2scan.mpl) +# LDA-correlation base: PW92 in the PW_MOD parametrization (lda_c_pw with both +# lda_c_pw_params and lda_c_pw_modified_params defined), reused from gga_c. +# ============================================================================ + +MGAMMA = (1.0 - math.log(2.0)) / math.pi**2 +DP2_C = 0.361 +# f(alpha) parameters set inside mgga_c_scan.mpl (distinct from exchange). +C1_C, C2_C, D_C = 0.64, 1.5, 0.7 + +# beta(rs) = beta_a*(1+beta_b*rs)/(1+beta_c*rs) (== regtpss / Hu-Langreth). +BETA_A, BETA_B, BETA_C = 0.066724550603149220, 0.1, 0.1778 + +# SCAN e0 LDA-0 pieces (gga_c_scan_e0.mpl / mgga_c_scan.mpl). +SCAN_B1C, SCAN_B2C, SCAN_B3C = 0.0285764, 0.0889, 0.125541 +SCAN_CHI_INF = 0.12802585262625815 +SCAN_G_CNST = 2.363 + +# rSCAN correlation switching polynomial: c[i] = coeff of alpha^i = fc[8-i]. +_RSCAN_FC = ( + 1.0, -0.64, -0.4352, -1.535685604549, + 3.061560252175, -1.915710236206, 0.516884468372, -0.051848879792, +) +# r2scan_dfc2 = sum_{i=1..7} i*ff[8-i] = sum_{j=1..7} j*c[j]. +_DFC2 = sum(j * _RSCAN_FC[j] for j in range(1, 8)) + + +def _scan_eclda0(rs): + return -SCAN_B1C / (1.0 + SCAN_B2C * jnp.sqrt(rs) + SCAN_B3C * rs) + + +def _scan_Gc(sp): + return (1.0 - SCAN_G_CNST * (TWO13 - 1.0) * f_zeta(sp)) * (1.0 - sp.z**12) + + +def _scan_g_infty(s): + return (1.0 + 4.0 * SCAN_CHI_INF * s**2) ** (-0.25) + + +def _scan_H0(rs, s): + ec0 = _scan_eclda0(rs) + return SCAN_B1C * jnp.log1p(jnp.expm1(-ec0 / SCAN_B1C) * (1.0 - _scan_g_infty(s))) + + +def _scan_e0(rs, sp, s): + return (_scan_eclda0(rs) + _scan_H0(rs, s)) * _scan_Gc(sp) + + +def _mbeta(rs): + return BETA_A * (1.0 + BETA_B * rs) / (1.0 + BETA_C * rs) + + +def _r2scan_d(sp): + return (opz_pow_n(sp.opz, 5.0 / 3.0) + opz_pow_n(sp.omz, 5.0 / 3.0)) / 2.0 + + +def _t_total(sp, ts0, ts1): + return ts0 * (sp.opz / 2.0) ** (5.0 / 3.0) + ts1 * (sp.omz / 2.0) ** (5.0 / 3.0) + + +@kernel("mgga_c_r2scan") +def mgga_c_r2scan(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1): + """Re-regularized SCAN correlation (Furness et al 2020).""" + eta = p["eta"] + phi = mphi(sp) + phi3 = phi**3 + s = X2S * TWO13 * xt # reduced gradient + t = tt(rs, sp, xt) # PBE reduced gradient t + + fpw = pw92_f(_PW_MOD, rs, sp) # LSDA1 = PW92 (S25) + + # w1 = exp(-ec/(gamma*phi^3)) - 1 >= 0; -> 0 as rs -> inf (S30). + w1 = jnp.expm1(-fpw / (MGAMMA * phi3)) + w1s = safe_max(w1, TINY) + + # LSDA0 and LSDA1 and their d/drs (S25, S28, S34). + def _elsda0(r): + return _scan_eclda0(r) * _scan_Gc(sp) + + def _elsda1(r): + return pw92_f(_PW_MOD, r, sp) + + elsda0 = _elsda0(rs) + elsda1 = fpw + delsda0 = jax.grad(_elsda0)(rs) + delsda1 = jax.grad(_elsda1)(rs) + + dz = _r2scan_d(sp) + # r2scan_dy: gradient correction to y (S34), Gaussian-damped in s. + dy = ( + _DFC2 + / (27.0 * MGAMMA * dz * phi3 * w1s) + * (20.0 * rs * (delsda0 - delsda1) - 45.0 * eta * (elsda0 - elsda1)) + * s**2 + * jnp.exp(-(s**4) / DP2_C**4) + ) + # r2scan_y (S32) and r2scan_g (S31). + y = _mbeta(rs) * t**2 / (MGAMMA * w1s) + base = safe_max(1.0 + 4.0 * (y - dy), TINY) + g = base ** (-0.25) + + # fH (S29); guard the log argument to > -1 for edge NaN-safety only. + inner = jnp.maximum(w1 * (1.0 - g), TINY - 1.0) + fH = MGAMMA * phi3 * jnp.log1p(inner) + + ec1 = fpw + fH # r2scan_ec1 (S24) + ec0 = _scan_e0(rs, sp, s) # r2scan_ec0 (S35-S41, == SCAN e0) + + alpha = (_t_total(sp, t0, t1) - xt**2 / 8.0) / ( + K_FACTOR_C * _t_total(sp, 1.0, 1.0) + eta * xt**2 / 8.0 + ) + fa = _r2scan_f_alpha(alpha, _RSCAN_FC, C1_C, C2_C, D_C) + + return ec1 + fa * (ec0 - ec1) diff --git a/src/funxc/kernels/scan.py b/src/funxc/kernels/scan.py new file mode 100644 index 0000000..6425b79 --- /dev/null +++ b/src/funxc/kernels/scan.py @@ -0,0 +1,174 @@ +"""SCAN meta-GGA exchange and correlation kernels. + +Forms follow libxc's maple/mgga_exc/mgga_x_scan.mpl and mgga_c_scan.mpl +verbatim (with the includes gga_c_scan_e0.mpl -> gga_c_regtpss.mpl -> +gga_c_pbe.mpl for the correlation's PBE-like ec1 limit). + +The numerically dangerous part of SCAN is the iso-orbital interpolation +f(alpha) = exp(-c1*alpha/(1-alpha)) (alpha <= 1, "left") + / -d*exp(c2/(1-alpha)) (alpha > 1, "right") +which is singular at alpha = 1. libxc truncates each branch to zero once the +exponential drops below machine epsilon; we reproduce that truncation with +the double-where idiom (arguments clamped *before* the reciprocal so values +AND all derivative orders stay finite at alpha -> 1, alpha = 0, alpha -> inf). +""" + +import math + +import jax.numpy as jnp + +from ..constants import K_FACTOR_C, MU_GE, TINY, X2S +from ..guards import Spin, f_zeta, mphi, safe_max, tt +from . import kernel +from .gga_c import _PW_MOD +from .lda import pw92_f +from .mgga_x import mgga_exchange + +# libxc's DBL_EPSILON and the |log eps| used to place the f(alpha) cutoffs. +DBL_EPSILON = 2.220446049250313e-16 +_MLOG_EPS = -math.log(DBL_EPSILON) # ~= 36.0437 + + +def scan_f_alpha(a, c1, c2, d): + """SCAN iso-orbital interpolation f(alpha) (maple mgga_x_scan.mpl). + + Left branch (a <= 1): exp(-c1*a/(1-a)), truncated to 0 for a above the + cutoff L/(L+c1) < 1 where the exponential reaches eps. + Right branch (a > 1): -d*exp(c2/(1-a)), truncated to 0 for a below the + cutoff (L'+c2)/L' > 1 (L' = -log(eps/|d|)). + + The min/max clamps keep 1-a bounded away from 0 in each branch, so the + exp arguments (and their derivatives) never see the singularity — the + double-where idiom, using the clamps as the input mask. + """ + left_cut = _MLOG_EPS / (_MLOG_EPS + c1) + ld = -math.log(DBL_EPSILON / abs(d)) + right_cut = (ld + c2) / ld + + a_l = jnp.minimum(a, left_cut) # <= left_cut < 1 => 1 - a_l > 0 + left0 = jnp.exp(-c1 * a_l / (1.0 - a_l)) + left = jnp.where(a > left_cut, 0.0, left0) + + a_r = jnp.maximum(a, right_cut) # >= right_cut > 1 => 1 - a_r < 0 + right0 = -d * jnp.exp(c2 / (1.0 - a_r)) + right = jnp.where(a < right_cut, 0.0, right0) + + return jnp.where(a <= 1.0, left, right) + + +# --- exchange ---------------------------------------------------------- + +# scan_y coefficients (maple mgga_x_scan.mpl); b4 depends on k1, computed +# per-call. b1, b2, b3 are pure constants. +SCAN_B2 = math.sqrt(5913.0 / 405000.0) +SCAN_B1 = (511.0 / 13500.0) / (2.0 * SCAN_B2) +SCAN_B3 = 0.5 +SCAN_A1 = 4.9479 +SCAN_H0X = 1.174 +# below this reduced gradient, gx = 1 - exp(-a1/sqrt(X2S*x)) is 1.0 to machine +# precision (the exp underflows exactly as it does in libxc); masking here also +# keeps the 1/sqrt(x) derivative finite (double-where). +_GX_X_MIN = 1e-10 + + +def scan_x_enhancement(p, x, u, t): + """SCAN exchange enhancement factor F(x, u, t) (u unused).""" + k1 = p["k1"] + c1, c2, d = p["c1"], p["c2"], p["d"] + + pp = X2S**2 * x**2 # = s^2 + alpha = (t - x**2 / 8.0) / K_FACTOR_C + b4 = MU_GE**2 / k1 - 1606.0 / 18225.0 - SCAN_B1**2 + + oma = 1.0 - alpha + y = ( + MU_GE * pp + + b4 * pp**2 * jnp.exp(-b4 * pp / MU_GE) + + (SCAN_B1 * pp + SCAN_B2 * oma * jnp.exp(-SCAN_B3 * oma**2)) ** 2 + ) + h1x = 1.0 + k1 * (1.0 - k1 / (k1 + y)) + fa = scan_f_alpha(alpha, c1, c2, d) + + # gx = 1 - exp(-a1/sqrt(X2S*x)); guard x against 0/subnormal. + xg = jnp.where(x > _GX_X_MIN, x, 1.0) + gx = jnp.where( + x > _GX_X_MIN, 1.0 - jnp.exp(-SCAN_A1 / jnp.sqrt(X2S * xg)), 1.0 + ) + + return (h1x * (1.0 - fa) + SCAN_H0X * fa) * gx + + +@kernel("mgga_x_scan") +def mgga_x_scan(p, rs, sp: Spin, xt, xs0, xs1, u0, u1, t0, t1): + """SCAN exchange (maple/mgga_exc/mgga_x_scan.mpl). + + Also serves MGGA_X_REVSCAN (parameter-only: c1/c2/d differ).""" + return mgga_exchange(scan_x_enhancement, p, rs, sp, xs0, xs1, u0, u1, t0, t1) + + +# --- correlation ------------------------------------------------------- + +# ec0 (LDA-like) limit constants, maple mgga_c_scan.mpl +SCAN_B1C = 0.0285764 +SCAN_B2C = 0.0889 +SCAN_B3C = 0.125541 +SCAN_CHI_INF = 0.12802585262625815 +SCAN_G_CNST = 2.363 # "in the paper it is 2.3631" + +# ec1 (PBE-like) limit: gga_c_pbe with the regtpss (Hu-Langreth, rs-dependent) +# beta override, over the *modified* PW92 LDA correlation base (_PW_MOD). +C_GAMMA = (1.0 - math.log(2.0)) / math.pi**2 +BETA_A = 0.066724550603149220 +BETA_B = 0.1 +BETA_C = 0.1778 + +# f(alpha) constants for correlation are hardcoded in mgga_c_scan.mpl +# (distinct from the exchange param table). +C_C1 = 0.64 +C_C2 = 1.5 +C_D = 0.7 + +_CBRT2 = 2.0 ** (1.0 / 3.0) + + +def _f_pbe_scan(rs, sp: Spin, xt): + """PBE-form correlation with regtpss rs-dependent beta (the ec1 limit).""" + ec = pw92_f(_PW_MOD, rs, sp) # modified-PW92 LDA correlation base + phi = mphi(sp) + t = tt(rs, sp, xt) + mbeta = BETA_A * (1.0 + BETA_B * rs) / (1.0 + BETA_C * rs) + denom = safe_max(jnp.expm1(-ec / (C_GAMMA * phi**3)), TINY) + a = mbeta / (C_GAMMA * denom) + f1 = t**2 + a * t**4 # BB = 1 + f2 = mbeta * f1 / (C_GAMMA * (1.0 + a * f1)) + h = C_GAMMA * phi**3 * jnp.log1p(f2) + return ec + h + + +@kernel("mgga_c_scan") +def mgga_c_scan(p, rs, sp: Spin, xt, xs0, xs1, u0, u1, t0, t1): + """SCAN correlation (maple/mgga_exc/mgga_c_scan.mpl). + + Interpolates the PBE-like ec1 and the SCAN ec0 limit with the same + iso-orbital f(alpha) switch. Also serves MGGA_C_SCAN_RVV10 (identical at + the semilocal level; VV10 is a nonlocal add-on outside the kernel).""" + ec1 = _f_pbe_scan(rs, sp, xt) + + # ec0 (LDA-like) limit: scan_e0(rs, z, s = X2S*2^(1/3)*xt) + s = X2S * _CBRT2 * xt + eclda0 = -SCAN_B1C / (1.0 + SCAN_B2C * jnp.sqrt(rs) + SCAN_B3C * rs) + g_inf = (1.0 + 4.0 * SCAN_CHI_INF * s**2) ** (-0.25) + h0 = SCAN_B1C * jnp.log1p((jnp.exp(-eclda0 / SCAN_B1C) - 1.0) * (1.0 - g_inf)) + fz = f_zeta(sp) + gc = (1.0 - SCAN_G_CNST * (_CBRT2 - 1.0) * fz) * (1.0 - sp.z**12) + e0 = (eclda0 + h0) * gc + + # alpha = (t_total - xt^2/8) / (K_FACTOR_C * t_total(z, 1, 1)) + opz2 = sp.opz / 2.0 + omz2 = sp.omz / 2.0 + t_tot = t0 * opz2 ** (5.0 / 3.0) + t1 * omz2 ** (5.0 / 3.0) + t_unif = opz2 ** (5.0 / 3.0) + omz2 ** (5.0 / 3.0) + alpha = (t_tot - xt**2 / 8.0) / (K_FACTOR_C * t_unif) + fa = scan_f_alpha(alpha, C_C1, C_C2, C_D) + + return ec1 + fa * (e0 - ec1) diff --git a/src/funxc/kernels/tpss.py b/src/funxc/kernels/tpss.py new file mode 100644 index 0000000..41675af --- /dev/null +++ b/src/funxc/kernels/tpss.py @@ -0,0 +1,171 @@ +"""TPSS meta-GGA exchange and correlation kernels. + +Ground truth: libxc/maple/mgga_exc/mgga_x_tpss.mpl + tpss_x.mpl (exchange), +mgga_exc/mgga_c_tpss.mpl + tpss_c.mpl (correlation), and util.mpl helpers. + +Exchange enhancement F(x, t) (Tao-Perdew-Staroverov-Scuseria 2003, Eq. 10) +plugs into the mgga_exchange spin-scaling combinator. Correlation is built on +PBE correlation (reused from gga_c.py) with the TPSS self-interaction +correction of Eq. 24-25: an "perp" term scaling PBE by (1 + C0*z^2) and a +"par" term subtracting the spin-resolved fully-polarized PBE pieces, gated by +the iso-orbital ratio z = tau_W/tau <= 1 (tpss_aux) and the C0(zeta) prefactor +whose z -> +-1 singularity is handled by libxc's 1-|z| <= 1e-12 piecewise. +""" + +import math + +import jax.numpy as jnp + +from ..constants import K_FACTOR_C, MU_GE, X2S, ZETA_THRESHOLD +from ..guards import Spin, opz_thr +from . import kernel +from .gga_c import gga_c_pbe +from .mgga_x import mgga_exchange + +X2S2 = X2S * X2S # X2S^2, converts x^2 -> s^2 (tpss_p) + + +# --- TPSS exchange ----------------------------------------------------------- + + +def _tpss_x_enhancement(p, x, u, t): + """F(x, t) from tpss_x.mpl (lapl u is unused). t = tau_s/rho_s^(5/3).""" + b, c, e = p["b"], p["c"], p["e"] + kappa, mu = p["kappa"], p["mu"] + bloc_a, bloc_b = p["BLOC_a"], p["BLOC_b"] + + x2 = x * x + p_s = X2S2 * x2 # tpss_p = s^2 + z = x2 / (8.0 * t) # tpss_z = x^2/(8 t) (t floored upstream, > 0) + alpha = (t - x2 / 8.0) / K_FACTOR_C + + qb = (9.0 / 20.0) * (alpha - 1.0) / jnp.sqrt( + 1.0 + b * alpha * (alpha - 1.0) + ) + (2.0 / 3.0) * p_s + + # z^tpss_ff(z) with tpss_ff = BLOC_a + BLOC_b*z. For the parameter-only + # variants implemented here BLOC_b == 0, so the exponent is the static + # float BLOC_a and z=0 (x->0) differentiates cleanly; the general + # BLOC_b != 0 branch guards z > 0 for the y-gradient of z^ff. + if bloc_b == 0.0: + z_pow = z ** bloc_a + else: + zc = jnp.where(z > 0.0, z, 1.0) + z_pow = jnp.where(z > 0.0, zc ** (bloc_a + bloc_b * z), 0.0) + + fxnum = ( + (MU_GE + c * z_pow / (1.0 + z * z) ** 2) * p_s + + (146.0 / 2025.0) * qb * qb + - (73.0 / 405.0) * qb * jnp.sqrt( + 0.5 * ((9.0 / 25.0) * z * z + p_s * p_s) + ) + + MU_GE**2 / kappa * p_s * p_s + + 2.0 * jnp.sqrt(e) * MU_GE * (9.0 / 25.0) * z * z + + e * mu * p_s**3 + ) + fxden = (1.0 + jnp.sqrt(e) * p_s) ** 2 + fx = fxnum / fxden + + # tpss_a1 = kappa/(kappa + fx); F = 1 + kappa*(1 - a1) = 1 + kappa*fx/(kappa+fx) + return 1.0 + kappa * fx / (kappa + fx) + + +@kernel("mgga_x_tpss") +def mgga_x_tpss(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1): + """TPSS exchange (mgga_x_tpss.mpl); covers MGGA_X_TPSS, MGGA_X_MODTPSS.""" + return mgga_exchange(_tpss_x_enhancement, p, rs, sp, xs0, xs1, u0, u1, t0, t1) + + +# --- TPSS correlation -------------------------------------------------------- + +# tpss_c.mpl sets params_a_gamma := (1 - log 2)/pi^2 and params_a_BB := 1 for +# the PBE base; beta comes from the functional params. The LDA-correlation +# base is therefore PW92 with the *modified* parametrization (PW_MOD), exactly +# as PBE correlation uses (gga_c.py::_PW_MOD) -- reused via gga_c_pbe below. +_TPSS_GAMMA = (1.0 - math.log(2.0)) / math.pi**2 +# xi2 denominator (2*(3 pi^2)^(1/3))^2 = 4 (3 pi^2)^(2/3) +_XI2_DEN = (2.0 * (3.0 * math.pi**2) ** (1.0 / 3.0)) ** 2 + + +def _pbe(pbe_p, rs, z, xt, xs0, xs1): + """PBE correlation eps_c evaluated at an arbitrary zeta scalar/array. + + gga_c_pbe reads only sp.z / sp.opz / sp.omz (opz_pow_n clamps zeta=+-1), + never the screen masks, so dummy False screens are safe here. + """ + sp = Spin(z=z, opz=1.0 + z, omz=1.0 - z, screen_up=False, screen_dn=False) + return gga_c_pbe(pbe_p, rs, sp, xt, xs0, xs1) + + +def _t_total(a0, a1, opz, omz): + """util.mpl t_total(z, a0, a1) = a0*((1+z)/2)^(5/3) + a1*((1-z)/2)^(5/3).""" + return a0 * (opz / 2.0) ** (5.0 / 3.0) + a1 * (omz / 2.0) ** (5.0 / 3.0) + + +@kernel("mgga_c_tpss") +def mgga_c_tpss(p, rs, sp, xt, xs0, xs1, u0, u1, t0, t1): + """TPSS correlation (mgga_c_tpss.mpl); MGGA_C_TPSS, MGGA_C_TPSS_GAUSSIAN.""" + pbe_p = {"beta": p["beta"], "gamma": _TPSS_GAMMA, "B": 1.0} + cc = p["C0_c"] + d = p["d"] + + z = sp.z + opz, omz = sp.opz, sp.omz + + # z_thr(z) and its 1+z_thr, 1-z_thr forms (opz_thr clamps into + # [zeta_threshold, 2 - zeta_threshold]). + opz_z = opz_thr(opz, omz) + omz_z = 2.0 - opz_z + z_t = opz_z - 1.0 + + # --- tpss_aux = min(xt^2/(8 t_total(z, t0, t1)), 1) (tau_W/tau clamp) --- + # t_total > 0 (per-spin tau floored upstream), so the ratio is finite; the + # min caps the iso-orbital ratio at 1 and its gradient there is 0. + t_tot_tau = _t_total(t0, t1, opz, omz) + aux = jnp.minimum(xt * xt / (8.0 * t_tot_tau), 1.0) + aux2 = aux * aux + aux3 = aux2 * aux + + # --- C0(zeta, xt, xs0, xs1) (Eqs. 33-34) with the 1-|z| <= 1e-12 branch --- + # The else branch has a singular first derivative as |z| -> 1; libxc caps + # it with a hard piecewise. Double-where on z (z_safe) keeps the discarded + # branch's value *and* gradient finite in the polarized region. + polar = (1.0 - jnp.abs(z)) <= 1e-12 + z_s = jnp.where(polar, 0.0, z) + opz_s = 1.0 + z_s + omz_s = 1.0 - z_s + opz_zs = opz_thr(opz_s, omz_s) + omz_zs = 2.0 - opz_zs + t_tot_x = _t_total(xs0 * xs0, xs1 * xs1, opz_s, omz_s) + xi2 = (1.0 - z_s * z_s) * (t_tot_x - xt * xt) / _XI2_DEN + c0_den = 1.0 + xi2 * (opz_zs ** (-4.0 / 3.0) + omz_zs ** (-4.0 / 3.0)) / 2.0 + c00 = cc[0] + cc[1] * z_s**2 + cc[2] * z_s**4 + cc[3] * z_s**6 + C0 = jnp.where(polar, cc[0] + cc[1] + cc[2] + cc[3], c00 / c0_den**4) + + # --- PBE base at raw zeta (perp) and at z_thr (par) --- + f_raw = _pbe(pbe_p, rs, z, xt, xs0, xs1) + f_zt = _pbe(pbe_p, rs, z_t, xt, xs0, xs1) + + # tpss_perp (first line of Eq. 25) + perp = (1.0 + C0 * aux2) * f_raw + + # tpss_par (second line of Eq. 25); per-channel Stoll-like decomposition + # with fully-polarized PBE evaluated at scaled rs, gated by density/zeta. + screen_up_z = jnp.logical_or(sp.screen_up, opz <= ZETA_THRESHOLD) + screen_dn_z = jnp.logical_or(sp.screen_dn, omz <= ZETA_THRESHOLD) + + rs_up = rs * (2.0 / opz_z) ** (1.0 / 3.0) + rs_dn = rs * (2.0 / omz_z) ** (1.0 / 3.0) + par_s0 = jnp.maximum( + _pbe(pbe_p, rs_up, 1.0, xs0, xs0, 0.0), f_zt + ) * (opz_z / 2.0) + par_s1 = jnp.maximum( + _pbe(pbe_p, rs_dn, -1.0, xs1, 0.0, xs1), f_zt + ) * (omz_z / 2.0) + term_up = jnp.where(screen_up_z, f_zt * (opz / 2.0), par_s0) + term_dn = jnp.where(screen_dn_z, f_zt * (omz / 2.0), par_s1) + par = -(1.0 + C0) * aux2 * (term_up + term_dn) + + # tpss_f0 (Eq. 25) and tpss_f (Eq. 24) + f0 = par + perp + return f0 * (1.0 + d * f0 * aux3) diff --git a/src/funxc/libxc_compat.py b/src/funxc/libxc_compat.py new file mode 100644 index 0000000..673c631 --- /dev/null +++ b/src/funxc/libxc_compat.py @@ -0,0 +1,74 @@ +"""Thin pylibxc-compatible shim. + +Mirrors pylibxc.LibXCFunctional's `compute` interface (dict in, dict of +numpy-shaped arrays out) so funxc can be dropped into hosts that already +speak pylibxc: + + from funxc.libxc_compat import LibXCFunctional + f = LibXCFunctional("gga_x_pbe", "polarized") + out = f.compute({"rho": rho, "sigma": sigma}) + out["zk"], out["vrho"], out["vsigma"] + +Array shapes follow pylibxc: zk (N, 1); polarized vrho (N, 2), vsigma (N, 3); +unpolarized vrho (N, 1), vsigma (N, 1). Flat interleaved input arrays are +accepted, as with pylibxc. +""" + +import numpy as np + +from .api import functional + + +class LibXCFunctional: + def __init__(self, func_name: str, spin): + if spin in (1, "unpolarized", "XC_UNPOLARIZED"): + self._polarized = False + elif spin in (2, "polarized", "XC_POLARIZED"): + self._polarized = True + else: + raise ValueError(f"invalid spin: {spin!r}") + self._f = functional(func_name, polarized=self._polarized) + + def get_name(self): + return self._f.spec.name + + def get_family(self): + return self._f.spec.family + + def get_number(self): + return self._f.spec.libxc_id + + def compute(self, inp, do_exc: bool = True, do_vxc: bool = True): + if not isinstance(inp, dict): + inp = {"rho": inp} + family = self._f.spec.family + rho = np.asarray(inp["rho"], dtype=np.float64) + if self._polarized: + rho = rho.reshape(-1, 2) + else: + rho = rho.reshape(-1) + npoints = rho.shape[0] + + def spin_field(key, ncol): + arr = np.asarray(inp[key], dtype=np.float64) + return arr.reshape(npoints, ncol) if self._polarized else arr.reshape(npoints) + + sigma = lapl = tau = None + if family != "lda": + sigma = spin_field("sigma", 3) + if family == "mgga": + tau = spin_field("tau", 2) + if "lapl" in inp: + lapl = spin_field("lapl", 2) + + out = {} + if do_vxc: + res = self._f.exc_vxc(rho, sigma, lapl, tau) + if do_exc: + out["zk"] = np.asarray(res["zk"]).reshape(npoints, 1) + for key in ("vrho", "vsigma", "vlapl", "vtau"): + if key in res: + out[key] = np.asarray(res[key]).reshape(npoints, -1) + elif do_exc: + out["zk"] = np.asarray(self._f.exc(rho, sigma, lapl, tau)).reshape(npoints, 1) + return out diff --git a/src/funxc/registry.py b/src/funxc/registry.py new file mode 100644 index 0000000..cb7248b --- /dev/null +++ b/src/funxc/registry.py @@ -0,0 +1,89 @@ +"""Functional database: TOML files in funxc/data -> FunctionalSpec registry. + +Adding a functional that reuses an existing kernel is one TOML entry; no code. +""" + +import tomllib +from dataclasses import dataclass, field +from importlib import resources +from types import MappingProxyType + +from .kernels import KERNELS + + +@dataclass(frozen=True) +class FunctionalSpec: + id: str # canonical libxc-style identifier, e.g. "GGA_X_PBE" + name: str # human-readable description + maple_name: str # libxc maple/ file the form was taken from + family: str # lda | gga | mgga + kind: str # x | c | xc | k + kernel: str # key into funxc.kernels.KERNELS + libxc_id: int # libxc functional number (regression-file header) + parameters: MappingProxyType = field(default_factory=lambda: MappingProxyType({})) + citation: str = "" + doi: str = "" + dens_threshold: float = 1e-15 + # meta-GGA flags, mirroring libxc's xc_func_info flags (work_mgga_inc.c): + # NEEDS_TAU gates the tau floor, NEEDS_LAPLACIAN marks real lapl use, + # ENFORCE_FHC clamps sigma <= 8*rho*tau (positive Fermi-hole curvature). + needs_tau: bool = False + needs_lapl: bool = False + enforce_fhc: bool = False + + @property + def kernel_fn(self): + return KERNELS[self.kernel] + + +def _freeze(params: dict) -> MappingProxyType: + return MappingProxyType( + {k: tuple(v) if isinstance(v, list) else v for k, v in params.items()} + ) + + +def _load_registry() -> dict[str, FunctionalSpec]: + registry: dict[str, FunctionalSpec] = {} + data_dir = resources.files("funxc.data") + for entry in sorted(data_dir.iterdir(), key=lambda e: e.name): + if not entry.name.endswith(".toml"): + continue + table = tomllib.loads(entry.read_text()) + for fid, spec in table.items(): + if fid in registry: + raise ValueError(f"duplicate functional id: {fid}") + if spec["kernel"] not in KERNELS: + raise ValueError(f"{fid}: unknown kernel {spec['kernel']!r}") + registry[fid] = FunctionalSpec( + id=fid, + name=spec["name"], + maple_name=spec["maple_name"], + family=spec["family"], + kind=spec["kind"], + kernel=spec["kernel"], + libxc_id=spec["libxc_id"], + parameters=_freeze(spec.get("parameters", {})), + citation=spec.get("citation", ""), + doi=spec.get("doi", ""), + dens_threshold=spec.get("dens_threshold", 1e-15), + needs_tau=spec.get("needs_tau", False), + needs_lapl=spec.get("needs_lapl", False), + enforce_fhc=spec.get("enforce_fhc", False), + ) + return registry + + +REGISTRY: dict[str, FunctionalSpec] = _load_registry() + + +def get_spec(name: str) -> FunctionalSpec: + """Look up a functional by libxc-style name, case-insensitively.""" + key = name.upper() + if key.startswith("XC_"): + key = key[3:] + try: + return REGISTRY[key] + except KeyError: + raise KeyError( + f"unknown functional {name!r}; available: {', '.join(sorted(REGISTRY))}" + ) from None diff --git a/tests/libxc_reference.py b/tests/libxc_reference.py new file mode 100644 index 0000000..7065cae --- /dev/null +++ b/tests/libxc_reference.py @@ -0,0 +1,67 @@ +"""Loaders for libxc's vendored testsuite data (golden reference). + +Input files (libxc/testsuite/input/): first line npoints, then +9 columns: rhoa rhob sigmaaa sigmaab sigmabb lapla laplb taua taub. + +Regression files (libxc/testsuite/regression//....bz2): +first line "func_id npoints order", then a column-header line, then values. +""" + +import bz2 +from pathlib import Path + +import numpy as np + +LIBXC_ROOT = Path(__file__).resolve().parent.parent / "libxc" +INPUT_DIR = LIBXC_ROOT / "testsuite" / "input" +REGRESSION_DIR = LIBXC_ROOT / "testsuite" / "regression" + +SYSTEMS = ("BrOH", "BrOH+", "H", "Li") + + +def load_input(system: str) -> dict[str, np.ndarray]: + """Return the 9-column grid data for one test system.""" + raw = np.loadtxt(INPUT_DIR / system, skiprows=1) + return { + "rho": raw[:, 0:2], # (N, 2): rhoa, rhob + "sigma": raw[:, 2:5], # (N, 3): sigmaaa, sigmaab, sigmabb + "lapl": raw[:, 5:7], + "tau": raw[:, 7:9], + } + + +def unpolarize(inp: dict[str, np.ndarray]) -> dict[str, np.ndarray]: + """Collapse spin channels the way xc-regression.c does for nspin=1.""" + return { + "rho": inp["rho"].sum(axis=1), + "sigma": inp["sigma"][:, 0] + 2.0 * inp["sigma"][:, 1] + inp["sigma"][:, 2], + "lapl": inp["lapl"].sum(axis=1), + "tau": inp["tau"].sum(axis=1), + } + + +def load_regression( + family_dir: str, functional: str, system: str, spin: str, order: int +) -> dict[str, np.ndarray]: + """Return reference outputs keyed by column name (zk, vrho(a), ...).""" + path = ( + REGRESSION_DIR + / family_dir + / f"{functional}.{system}.{spin}.{order}.bz2" + ) + with bz2.open(path, "rt") as fh: + header = fh.readline().split() + func_id, npoints, forder = int(header[0]), int(header[1]), int(header[2]) + assert forder == order + names = fh.readline().split() + data = np.loadtxt(fh) + data = np.asarray(data).reshape(npoints, len(names)) + out = {name: data[:, i] for i, name in enumerate(names)} + out["_func_id"] = func_id + return out + + +def has_regression(family_dir: str, functional: str) -> bool: + return any( + (REGRESSION_DIR / family_dir).glob(f"{functional}.*.bz2") + ) diff --git a/tests/test_api.py b/tests/test_api.py new file mode 100644 index 0000000..357c9f0 --- /dev/null +++ b/tests/test_api.py @@ -0,0 +1,84 @@ +"""API surface tests: registry integrity and the pylibxc-compatible shim.""" + +import numpy as np +import pytest + +import funxc +from funxc.libxc_compat import LibXCFunctional + + +def test_registry_specs_wellformed(): + for fid, spec in funxc.REGISTRY.items(): + assert spec.id == fid + assert spec.family in ("lda", "gga", "mgga") + assert spec.kind in ("x", "c", "xc", "k") + assert spec.libxc_id > 0 + assert spec.kernel_fn is not None + assert spec.doi + + +def test_get_spec_is_case_insensitive(): + assert funxc.get_spec("gga_x_pbe") is funxc.get_spec("GGA_X_PBE") + assert funxc.get_spec("XC_GGA_X_PBE") is funxc.get_spec("GGA_X_PBE") + with pytest.raises(KeyError): + funxc.get_spec("GGA_X_DOES_NOT_EXIST") + + +def test_shim_polarized_shapes(): + f = LibXCFunctional("gga_x_pbe", "polarized") + n = 4 + rng = np.random.default_rng(0) + inp = { + "rho": rng.uniform(0.1, 1.0, (n, 2)), + "sigma": rng.uniform(0.0, 0.5, (n, 3)), + } + out = f.compute(inp) + assert out["zk"].shape == (n, 1) + assert out["vrho"].shape == (n, 2) + assert out["vsigma"].shape == (n, 3) + + +def test_shim_unpolarized_shapes_and_flat_input(): + f = LibXCFunctional("lda_x", 1) + rho = np.linspace(0.1, 2.0, 5) + out = f.compute({"rho": rho}) + assert out["zk"].shape == (5, 1) + assert out["vrho"].shape == (5, 1) + out_exc_only = f.compute({"rho": rho}, do_vxc=False) + assert set(out_exc_only) == {"zk"} + + +def test_shim_mgga_shapes(): + f = LibXCFunctional("mgga_x_lta", "polarized") + n = 4 + rng = np.random.default_rng(1) + inp = { + "rho": rng.uniform(0.1, 1.0, (n, 2)), + "sigma": rng.uniform(0.0, 0.5, (n, 3)), + "lapl": rng.uniform(-1.0, 1.0, (n, 2)), + "tau": rng.uniform(0.1, 1.0, (n, 2)), + } + out = f.compute(inp) + assert out["zk"].shape == (n, 1) + assert out["vrho"].shape == (n, 2) + assert out["vsigma"].shape == (n, 3) + assert out["vlapl"].shape == (n, 2) + assert out["vtau"].shape == (n, 2) + # LTA has no lapl dependence: vlapl must be exactly zero + assert np.all(out["vlapl"] == 0.0) + + +def test_shim_metadata(): + f = LibXCFunctional("gga_c_lyp", "polarized") + assert f.get_number() == 131 + assert f.get_family() == "gga" + + +def test_composition_blyp(): + """X and C functionals compose additively (BLYP = B88 + LYP).""" + rho = np.array([[0.4, 0.3]]) + sigma = np.array([[0.02, 0.01, 0.015]]) + b88 = funxc.functional("GGA_X_B88").exc(rho, sigma) + lyp = funxc.functional("GGA_C_LYP").exc(rho, sigma) + total = np.asarray(b88) + np.asarray(lyp) + assert np.all(np.isfinite(total)) and total[0] < 0 diff --git a/tests/test_guards.py b/tests/test_guards.py new file mode 100644 index 0000000..4effc58 --- /dev/null +++ b/tests/test_guards.py @@ -0,0 +1,100 @@ +"""NaN-safety and differentiability tests for the guard layer. + +The regression tests establish numeric parity; these establish that funxc +stays finite (values *and* gradients) at the singular corners libxc's grids +avoid: rho -> 0 tails, sigma = 0, full spin polarization, and combinations. +""" + +import jax +import jax.numpy as jnp +import numpy as np +import pytest + +import funxc + +GGA_IDS = [fid for fid in funxc.REGISTRY if funxc.REGISTRY[fid].family == "gga"] +LDA_IDS = [fid for fid in funxc.REGISTRY if funxc.REGISTRY[fid].family == "lda"] + +# Adversarial polarized points: +# (rho_up, rho_dn, sigma_uu, sigma_ud, sigma_dd, tau_up, tau_dn) +EDGE_POINTS = [ + (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # vacuum + (1e-300, 0.0, 0.0, 0.0, 0.0, 1e-300, 0.0), # denormal tail + (0.3, 0.0, 7e-3, 0.0, 0.0, 2.9e-3, 0.0), # fully polarized, tau ~ tau_W + (0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0), # sigma = tau = 0 exactly + (0.3, 1e-16, 1e-2, 0.0, 1e-30, 5e-2, 1e-30), # one channel below thr + (1e-14, 1e-14, 1e-25, -1e-25, 1e-25, 1e-25, 1e-25), # negative sigma_ud + (1e3, 1e3, 1e8, 1e8, 1e8, 1e8, 1e8), # high density, huge gradient + (0.5, 0.5, 1e12, 0.0, 1e12, 1e-5, 1e-5), # s -> inf, tau << tau_W (FHC) +] + + +def _edge_args(spec, f): + rho = np.array([[p[0], p[1]] for p in EDGE_POINTS]) + sigma = np.array([[p[2], p[3], p[4]] for p in EDGE_POINTS]) + tau = np.array([[p[5], p[6]] for p in EDGE_POINTS]) + lapl = np.zeros_like(tau) + if spec.family == "lda": + return (rho,) + if spec.family == "gga": + return (rho, sigma) + return (rho, sigma, lapl, tau) + + +@pytest.mark.parametrize("fid", sorted(funxc.REGISTRY)) +def test_values_and_gradients_finite_at_edges(fid): + spec = funxc.REGISTRY[fid] + f = funxc.functional(fid, polarized=True) + out = f.exc_vxc(*_edge_args(spec, f)) + for key, val in out.items(): + assert np.all(np.isfinite(val)), f"{fid}: non-finite {key} at edge points" + + +@pytest.mark.parametrize( + "fid", + ["LDA_X", "LDA_C_PW", "GGA_X_PBE", "GGA_C_PBE", "GGA_C_LYP", + "GGA_XC_B97_D", "MGGA_X_LTA"], +) +def test_second_derivatives_finite(fid): + """fxc-level derivatives must also stay finite (autodiff twice).""" + spec = funxc.REGISTRY[fid] + f = funxc.functional(fid, polarized=True) + if spec.family == "lda": + hess = jax.hessian(f.energy_density, argnums=0) + args = lambda p: (jnp.array(p[:2]),) # noqa: E731 + elif spec.family == "gga": + hess = jax.hessian(f.energy_density, argnums=(0, 1)) + args = lambda p: (jnp.array(p[:2]), jnp.array(p[2:5])) # noqa: E731 + else: + hess = jax.hessian(f.energy_density, argnums=(0, 1, 2, 3)) + args = lambda p: ( # noqa: E731 + jnp.array(p[:2]), jnp.array(p[2:5]), jnp.zeros(2), jnp.array(p[5:7]) + ) + for p in EDGE_POINTS: + h = hess(*args(p)) + leaves = jax.tree.leaves(h) + assert all(jnp.all(jnp.isfinite(x)) for x in leaves), ( + f"{fid}: non-finite second derivative at {p}" + ) + + +def test_vacuum_outputs_are_zero(): + f = funxc.functional("GGA_X_PBE", polarized=True) + out = f.exc_vxc(np.zeros((1, 2)), np.zeros((1, 3))) + for key, val in out.items(): + assert np.all(np.asarray(val) == 0.0), f"nonzero {key} in vacuum" + + +def test_jit_and_vmap_compose(): + """The per-point eps must survive user-side jit/vmap/grad composition.""" + f = funxc.functional("GGA_C_PBE", polarized=True) + g = jax.jit(jax.vmap(jax.grad(f.energy_density, argnums=(0, 1)))) + rho = jnp.array([[0.3, 0.2], [1.0, 1.0]]) + sigma = jnp.array([[0.01, 0.0, 0.02], [0.1, 0.05, 0.1]]) + vrho, vsigma = g(rho, sigma) + assert vrho.shape == (2, 2) and vsigma.shape == (2, 3) + assert jnp.all(jnp.isfinite(vrho)) and jnp.all(jnp.isfinite(vsigma)) + + +def test_float64_active(): + assert jnp.asarray(1.0).dtype == jnp.float64 diff --git a/tests/test_regression.py b/tests/test_regression.py new file mode 100644 index 0000000..af6bfcf --- /dev/null +++ b/tests/test_regression.py @@ -0,0 +1,134 @@ +"""Golden-data tests: funxc vs libxc 7.0.0 vendored regression files. + +Every registered functional is compared against libxc/testsuite/regression +for all four systems, both spin modes, orders 0 (zk) and 1 (vrho, vsigma). + +Tolerance policy: |got - ref| <= ATOL + RTOL*|ref|. +Typical agreement is machine precision (zk ~1e-15 relative). The absolute +floor covers two known noise sources that libxc itself carries: + * roundoff residue of analytically-cancelling terms at zeta = +-1 (the + reference files store values like -1.1e-16 where the exact result is 0); + * expression-grouping noise ~(A*t^2)^2 * eps in the saturated large-t + limit of the PBE H(t) correlation term, visible only at density-tail + points (rho ~ 1e-11) where |vsigma| ~ 1e-6 and the noise is ~4e-10. +""" + +import numpy as np +import pytest + +import funxc + +from libxc_reference import ( + SYSTEMS, + has_regression, + load_input, + load_regression, + unpolarize, +) + +ATOL = 2e-9 +RTOL = 1e-7 + +# Points where the vendored regression *file* disagrees with the libxc 7.0.0 +# *library* itself (module libxc/7.0.0): the testsuite files ship with the +# source tree and some predate behavior changes. Each entry was verified by +# comparing the installed .so against the file — funxc matches the library +# (e.g. LTA Li vrho(b) point 4: library -2.38273e-8, funxc -2.38273e-8, +# stale file -1.05e-14, a leftover of an older density screen). +STALE_GOLDEN = { + ("MGGA_X_LTA", "Li", "pol", "vrho(b)", 4), +} + + +def family_dir(fid: str) -> str: + return "_".join(fid.lower().split("_")[:2]) + + +FUNCTIONALS = [ + fid for fid in sorted(funxc.REGISTRY) if has_regression(family_dir(fid), fid.lower()) +] + + +def check(got, ref, label, stale=()): + got = np.asarray(got).reshape(-1) + err = np.abs(got - ref) + tol = ATOL + RTOL * np.abs(ref) + bad = err > tol + if len(stale): + bad[np.asarray(stale, dtype=int)] = False + assert not bad.any(), ( + f"{label}: {bad.sum()}/{bad.size} points exceed tolerance; " + f"worst at i={err.argmax()}: got={got[err.argmax()]:.12e} " + f"ref={ref[err.argmax()]:.12e}" + ) + + +@pytest.mark.parametrize("spin", ["pol", "unpol"]) +@pytest.mark.parametrize("system", SYSTEMS) +@pytest.mark.parametrize("fid", FUNCTIONALS) +def test_against_libxc(fid, system, spin): + fdir = family_dir(fid) + f = funxc.functional(fid, polarized=(spin == "pol")) + + inp = load_input(system) + if spin == "unpol": + inp = unpolarize(inp) + rho, sigma = inp["rho"], inp["sigma"] + + if f.spec.family == "lda": + out = f.exc_vxc(rho) + elif f.spec.family == "gga": + out = f.exc_vxc(rho, sigma) + else: + out = f.exc_vxc(rho, sigma, inp["lapl"], inp["tau"]) + + try: + ref0 = load_regression(fdir, fid.lower(), system, spin, 0) + except FileNotFoundError: + # libxc keeps some combinations in regression//disabled + # (e.g. gga_c_lyp on the fully polarized H atom); follow suit. + pytest.skip(f"libxc testsuite has no {fid} {system} {spin} data") + assert ref0["_func_id"] == f.spec.libxc_id + check(out["zk"], ref0["zk"], f"{fid} {system} {spin} zk") + + try: + ref1 = load_regression(fdir, fid.lower(), system, spin, 1) + except FileNotFoundError: + # Some combinations disable only the derivative data (e.g. + # gga_xc_b97_3c on polarized H); zk was still checked above. + pytest.skip(f"libxc testsuite has no {fid} {system} {spin} order-1 data") + if spin == "pol": + columns = { + "vrho": ("(a)", "(b)"), + "vsigma": ("(aa)", "(ab)", "(bb)"), + "vlapl": ("(a)", "(b)"), + "vtau": ("(a)", "(b)"), + } + for key, comps in columns.items(): + if key not in out or f"{key}{comps[0]}" not in ref1: + continue + got = np.asarray(out[key]) + for i, comp in enumerate(comps): + stale = [ + idx + for (sf, ss, sm, sc, idx) in STALE_GOLDEN + if (sf, ss, sm, sc) == (fid, system, spin, f"{key}{comp}") + ] + check( + got[:, i], ref1[f"{key}{comp}"], + f"{fid} {system} {key}{comp}", stale=stale, + ) + else: + for key in ("vrho", "vsigma", "vlapl", "vtau"): + if key in out and key in ref1: + stale = [ + idx + for (sf, ss, sm, sc, idx) in STALE_GOLDEN + if (sf, ss, sm, sc) == (fid, system, spin, key) + ] + check(out[key], ref1[key], f"{fid} {system} {key}", stale=stale) + + +def test_all_registered_functionals_have_golden_data(): + missing = [fid for fid in sorted(funxc.REGISTRY) if fid not in FUNCTIONALS] + assert not missing, f"no regression data for: {missing}" diff --git a/uv.lock b/uv.lock new file mode 100644 index 0000000..e2a3d7b --- /dev/null +++ b/uv.lock @@ -0,0 +1,250 @@ +version = 1 +revision = 3 +requires-python = ">=3.13" + +[[package]] +name = "colorama" +version = "0.4.6" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d8/53/6f443c9a4a8358a93a6792e2acffb9d9d5cb0a5cfd8802644b7b1c9a02e4/colorama-0.4.6.tar.gz", hash = "sha256:08695f5cb7ed6e0531a20572697297273c47b8cae5a63ffc6d6ed5c201be6e44", size = 27697, upload-time = "2022-10-25T02:36:22.414Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/d1/d6/3965ed04c63042e047cb6a3e6ed1a63a35087b6a609aa3a15ed8ac56c221/colorama-0.4.6-py2.py3-none-any.whl", hash = "sha256:4f1d9991f5acc0ca119f9d443620b77f9d6b33703e51011c16baf57afb285fc6", size = 25335, upload-time = "2022-10-25T02:36:20.889Z" }, +] + +[[package]] +name = "funxc" +version = "0.1.0" +source = { editable = "." } +dependencies = [ + { name = "jax" }, + { name = "numpy" }, +] + +[package.dev-dependencies] +dev = [ + { name = "pytest" }, +] + +[package.metadata] +requires-dist = [ + { name = "jax", specifier = ">=0.5" }, + { name = "numpy", specifier = ">=2.0" }, +] + +[package.metadata.requires-dev] +dev = [{ name = "pytest", specifier = ">=8" }] + +[[package]] +name = "iniconfig" +version = "2.3.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/72/34/14ca021ce8e5dfedc35312d08ba8bf51fdd999c576889fc2c24cb97f4f10/iniconfig-2.3.0.tar.gz", hash = "sha256:c76315c77db068650d49c5b56314774a7804df16fee4402c1f19d6d15d8c4730", size = 20503, upload-time = "2025-10-18T21:55:43.219Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl", hash = "sha256:f631c04d2c48c52b84d0d0549c99ff3859c98df65b3101406327ecc7d53fbf12", size = 7484, upload-time = "2025-10-18T21:55:41.639Z" }, +] + +[[package]] +name = "jax" +version = "0.10.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "jaxlib" }, + { name = "ml-dtypes" }, + { name = "numpy" }, + { name = "opt-einsum" }, + { name = "scipy" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/d2/73/eb91d98fcadfa2cbcfdd4e417ab116e47eb20882acc5ee678e47c35d6b57/jax-0.10.2.tar.gz", hash = "sha256:bf77428a8c2e6904c4f46d5ab12aa5cfc6cad2179f07f7e4c0fc75ac86ef0639", size = 2775110, upload-time = "2026-06-17T23:44:57.818Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/89/82/5ab5211079a151b6f661529369c0c8e98ec64cabf5c0cf22a0a05af124d8/jax-0.10.2-py3-none-any.whl", hash = "sha256:724d73c4678d8b06f6a6ab4db1b8a2fea8cd4f1e2c2564f99601634ec7b8d1c6", size = 3219515, upload-time = "2026-06-17T23:42:41.259Z" }, +] + +[[package]] +name = "jaxlib" +version = "0.10.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "ml-dtypes" }, + { name = "numpy" }, + { name = "scipy" }, +] +wheels = [ + { url = "https://files.pythonhosted.org/packages/c2/68/eaa4cebe253359196a8e80a33b242959e27d8d2a6ae3d09339f21da2acb8/jaxlib-0.10.2-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:4df530afa354a22dc1747a5d560640450cbb895d49889338a3f58c76a4c76c8e", size = 61434805, upload-time = "2026-06-17T23:44:10.511Z" }, + { url = "https://files.pythonhosted.org/packages/25/c1/4b884ea5962b6beb3c0f93742db54246bbf8b3274e48b0aca47908e454be/jaxlib-0.10.2-cp313-cp313-manylinux_2_27_aarch64.whl", hash = "sha256:45b28b0238697ab74bbcf20411aafb6db42acc31836cc2fd711e5cf056bf9556", size = 81084260, upload-time = "2026-06-17T23:44:14.065Z" }, + { url = "https://files.pythonhosted.org/packages/36/ac/4ee28c65861605945223145fbcd3c9362ec2255ddf7d917574e205548c82/jaxlib-0.10.2-cp313-cp313-manylinux_2_27_x86_64.whl", hash = "sha256:9e4818b4a8756fd3918766ca2aa5342125809f4f08a6fe46026d4386e7c23644", size = 85467706, upload-time = "2026-06-17T23:44:17.471Z" }, + { url = "https://files.pythonhosted.org/packages/79/54/9918b0f77a25a1299818c0610305ca2bea38ed90584f4489b60357e2dd39/jaxlib-0.10.2-cp313-cp313-win_amd64.whl", hash = "sha256:c75d6f1df1c9cff08e110b4a21c79560fdc502f4288972d6b117d25dafd44352", size = 65897894, upload-time = "2026-06-17T23:44:21.153Z" }, + { url = "https://files.pythonhosted.org/packages/56/5b/70df11da52a8b1a826184cccc05a3fec8aed76058a980021873fba3069cb/jaxlib-0.10.2-cp313-cp313t-macosx_11_0_arm64.whl", hash = "sha256:4c202d8ff7c1f3b5049dbd8f1e30e52759cd4e0a5835f0b3c7ae076a05818e28", size = 61564746, upload-time = "2026-06-17T23:44:24.421Z" }, + { url = "https://files.pythonhosted.org/packages/23/5c/184a648ea5db6c8b1a08fc5784c157b4c557255e009fb56091393df3c6de/jaxlib-0.10.2-cp313-cp313t-manylinux_2_27_aarch64.whl", hash = "sha256:b7b029bb95d981566750475b9719a9d6b66ed5dd2748851667899b6cfe075299", size = 81204888, upload-time = "2026-06-17T23:44:27.57Z" }, + { url = "https://files.pythonhosted.org/packages/6c/5c/539596a55265711d74147913278bcdc38412980be7d74d9c9d860297c486/jaxlib-0.10.2-cp313-cp313t-manylinux_2_27_x86_64.whl", hash = "sha256:e8b126097d609b0c6e6786e89f6dd6978adc02ebd5f63a1c61293fbac7821305", size = 85583810, upload-time = "2026-06-17T23:44:30.962Z" }, + { url = "https://files.pythonhosted.org/packages/f2/0d/27471ec9f1d04674f6e62de809412371e097aed3eca7d9483e677c54c214/jaxlib-0.10.2-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:72eba28b12fee02616fa42aa4b881b4ab62d7757c7843c462401d3fb34a27be4", size = 61446097, upload-time = "2026-06-17T23:44:34.196Z" }, + { url = "https://files.pythonhosted.org/packages/af/c8/941a7f7f37510f51290a5bd1a413aeef977fb8ba8adc0cfe8391233a764c/jaxlib-0.10.2-cp314-cp314-manylinux_2_27_aarch64.whl", hash = "sha256:f18f56fee90699cfba9b6627045a7a299702cb0e2af82ce180d9a6a7c8048093", size = 81096546, upload-time = "2026-06-17T23:44:37.486Z" }, + { url = "https://files.pythonhosted.org/packages/69/77/ac054882c220872512df28d16aeb648fe0e651efbb5be4fd7c4817fd88b0/jaxlib-0.10.2-cp314-cp314-manylinux_2_27_x86_64.whl", hash = "sha256:ca34f363197fb0ac4082582ca755007910369e33f8a8ba3d35ed94b71070107d", size = 85472993, upload-time = "2026-06-17T23:44:40.904Z" }, + { url = "https://files.pythonhosted.org/packages/0d/7d/c592d1fa69c210be0d2743fffc598dfc2f54efa9671c5f6f5d1e151c6f4a/jaxlib-0.10.2-cp314-cp314-win_amd64.whl", hash = "sha256:99818b0a18adc0b899abf4873795e8d65169441d87ab2e5cbb228e73d0f25808", size = 68376553, upload-time = "2026-06-17T23:44:44.968Z" }, + { url = "https://files.pythonhosted.org/packages/54/9b/91b00ec74985d29708b50420b4103c1f651c8f1c253d4fcb49d1bbb532cd/jaxlib-0.10.2-cp314-cp314t-macosx_11_0_arm64.whl", hash = "sha256:fc62997fce8831819551a2a5469a818169b09582b5b648c102d11ac7205bb812", size = 61564620, upload-time = "2026-06-17T23:44:48.217Z" }, + { url = "https://files.pythonhosted.org/packages/a1/c7/49d2b19c3b3105c30e1d3af2062e82e1977fb239d3dc3cbb583ef676dda7/jaxlib-0.10.2-cp314-cp314t-manylinux_2_27_aarch64.whl", hash = "sha256:a24d6e3cba263978293eae8b41330d5ccf24d6cdd1a6bcd4e82aff34e767620d", size = 81206365, upload-time = "2026-06-17T23:44:51.419Z" }, + { url = "https://files.pythonhosted.org/packages/bf/99/006cedf443f4a01f2088651facce79b2105bfb4905bfe9162eb0920a6dfb/jaxlib-0.10.2-cp314-cp314t-manylinux_2_27_x86_64.whl", hash = "sha256:5a2ac7aed7c4e661f67600bbcdec9e589151c1efec91f4cdb8d484af1a45c895", size = 85584458, upload-time = "2026-06-17T23:44:55.377Z" }, +] + +[[package]] +name = "ml-dtypes" +version = "0.5.4" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/0e/4a/c27b42ed9b1c7d13d9ba8b6905dece787d6259152f2309338aed29b2447b/ml_dtypes-0.5.4.tar.gz", hash = "sha256:8ab06a50fb9bf9666dd0fe5dfb4676fa2b0ac0f31ecff72a6c3af8e22c063453", size = 692314, upload-time = "2025-11-17T22:32:31.031Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/d9/a1/4008f14bbc616cfb1ac5b39ea485f9c63031c4634ab3f4cf72e7541f816a/ml_dtypes-0.5.4-cp313-cp313-macosx_10_13_universal2.whl", hash = "sha256:8c760d85a2f82e2bed75867079188c9d18dae2ee77c25a54d60e9cc79be1bc48", size = 676888, upload-time = "2025-11-17T22:31:56.907Z" }, + { url = "https://files.pythonhosted.org/packages/d3/b7/dff378afc2b0d5a7d6cd9d3209b60474d9819d1189d347521e1688a60a53/ml_dtypes-0.5.4-cp313-cp313-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:ce756d3a10d0c4067172804c9cc276ba9cc0ff47af9078ad439b075d1abdc29b", size = 5036993, upload-time = "2025-11-17T22:31:58.497Z" }, + { url = "https://files.pythonhosted.org/packages/eb/33/40cd74219417e78b97c47802037cf2d87b91973e18bb968a7da48a96ea44/ml_dtypes-0.5.4-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:533ce891ba774eabf607172254f2e7260ba5f57bdd64030c9a4fcfbd99815d0d", size = 5010956, upload-time = "2025-11-17T22:31:59.931Z" }, + { url = "https://files.pythonhosted.org/packages/e1/8b/200088c6859d8221454825959df35b5244fa9bdf263fd0249ac5fb75e281/ml_dtypes-0.5.4-cp313-cp313-win_amd64.whl", hash = "sha256:f21c9219ef48ca5ee78402d5cc831bd58ea27ce89beda894428bc67a52da5328", size = 212224, upload-time = "2025-11-17T22:32:01.349Z" }, + { url = "https://files.pythonhosted.org/packages/8f/75/dfc3775cb36367816e678f69a7843f6f03bd4e2bcd79941e01ea960a068e/ml_dtypes-0.5.4-cp313-cp313-win_arm64.whl", hash = "sha256:35f29491a3e478407f7047b8a4834e4640a77d2737e0b294d049746507af5175", size = 160798, upload-time = "2025-11-17T22:32:02.864Z" }, + { url = "https://files.pythonhosted.org/packages/4f/74/e9ddb35fd1dd43b1106c20ced3f53c2e8e7fc7598c15638e9f80677f81d4/ml_dtypes-0.5.4-cp313-cp313t-macosx_10_13_universal2.whl", hash = "sha256:304ad47faa395415b9ccbcc06a0350800bc50eda70f0e45326796e27c62f18b6", size = 702083, upload-time = "2025-11-17T22:32:04.08Z" }, + { url = "https://files.pythonhosted.org/packages/74/f5/667060b0aed1aa63166b22897fdf16dca9eb704e6b4bbf86848d5a181aa7/ml_dtypes-0.5.4-cp313-cp313t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:6a0df4223b514d799b8a1629c65ddc351b3efa833ccf7f8ea0cf654a61d1e35d", size = 5354111, upload-time = "2025-11-17T22:32:05.546Z" }, + { url = "https://files.pythonhosted.org/packages/40/49/0f8c498a28c0efa5f5c95a9e374c83ec1385ca41d0e85e7cf40e5d519a21/ml_dtypes-0.5.4-cp313-cp313t-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:531eff30e4d368cb6255bc2328d070e35836aa4f282a0fb5f3a0cd7260257298", size = 5366453, upload-time = "2025-11-17T22:32:07.115Z" }, + { url = "https://files.pythonhosted.org/packages/8c/27/12607423d0a9c6bbbcc780ad19f1f6baa2b68b18ce4bddcdc122c4c68dc9/ml_dtypes-0.5.4-cp313-cp313t-win_amd64.whl", hash = "sha256:cb73dccfc991691c444acc8c0012bee8f2470da826a92e3a20bb333b1a7894e6", size = 225612, upload-time = "2025-11-17T22:32:08.615Z" }, + { url = "https://files.pythonhosted.org/packages/e5/80/5a5929e92c72936d5b19872c5fb8fc09327c1da67b3b68c6a13139e77e20/ml_dtypes-0.5.4-cp313-cp313t-win_arm64.whl", hash = "sha256:3bbbe120b915090d9dd1375e4684dd17a20a2491ef25d640a908281da85e73f1", size = 164145, upload-time = "2025-11-17T22:32:09.782Z" }, + { url = "https://files.pythonhosted.org/packages/72/4e/1339dc6e2557a344f5ba5590872e80346f76f6cb2ac3dd16e4666e88818c/ml_dtypes-0.5.4-cp314-cp314-macosx_10_13_universal2.whl", hash = "sha256:2b857d3af6ac0d39db1de7c706e69c7f9791627209c3d6dedbfca8c7e5faec22", size = 673781, upload-time = "2025-11-17T22:32:11.364Z" }, + { url = "https://files.pythonhosted.org/packages/04/f9/067b84365c7e83bda15bba2b06c6ca250ce27b20630b1128c435fb7a09aa/ml_dtypes-0.5.4-cp314-cp314-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:805cef3a38f4eafae3a5bf9ebdcdb741d0bcfd9e1bd90eb54abd24f928cd2465", size = 5036145, upload-time = "2025-11-17T22:32:12.783Z" }, + { url = "https://files.pythonhosted.org/packages/c6/bb/82c7dcf38070b46172a517e2334e665c5bf374a262f99a283ea454bece7c/ml_dtypes-0.5.4-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:14a4fd3228af936461db66faccef6e4f41c1d82fcc30e9f8d58a08916b1d811f", size = 5010230, upload-time = "2025-11-17T22:32:14.38Z" }, + { url = "https://files.pythonhosted.org/packages/e9/93/2bfed22d2498c468f6bcd0d9f56b033eaa19f33320389314c19ef6766413/ml_dtypes-0.5.4-cp314-cp314-win_amd64.whl", hash = "sha256:8c6a2dcebd6f3903e05d51960a8058d6e131fe69f952a5397e5dbabc841b6d56", size = 221032, upload-time = "2025-11-17T22:32:15.763Z" }, + { url = "https://files.pythonhosted.org/packages/76/a3/9c912fe6ea747bb10fe2f8f54d027eb265db05dfb0c6335e3e063e74e6e8/ml_dtypes-0.5.4-cp314-cp314-win_arm64.whl", hash = "sha256:5a0f68ca8fd8d16583dfa7793973feb86f2fbb56ce3966daf9c9f748f52a2049", size = 163353, upload-time = "2025-11-17T22:32:16.932Z" }, + { url = "https://files.pythonhosted.org/packages/cd/02/48aa7d84cc30ab4ee37624a2fd98c56c02326785750cd212bc0826c2f15b/ml_dtypes-0.5.4-cp314-cp314t-macosx_10_13_universal2.whl", hash = "sha256:bfc534409c5d4b0bf945af29e5d0ab075eae9eecbb549ff8a29280db822f34f9", size = 702085, upload-time = "2025-11-17T22:32:18.175Z" }, + { url = "https://files.pythonhosted.org/packages/5a/e7/85cb99fe80a7a5513253ec7faa88a65306be071163485e9a626fce1b6e84/ml_dtypes-0.5.4-cp314-cp314t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:2314892cdc3fcf05e373d76d72aaa15fda9fb98625effa73c1d646f331fcecb7", size = 5355358, upload-time = "2025-11-17T22:32:19.7Z" }, + { url = "https://files.pythonhosted.org/packages/79/2b/a826ba18d2179a56e144aef69e57fb2ab7c464ef0b2111940ee8a3a223a2/ml_dtypes-0.5.4-cp314-cp314t-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:0d2ffd05a2575b1519dc928c0b93c06339eb67173ff53acb00724502cda231cf", size = 5366332, upload-time = "2025-11-17T22:32:21.193Z" }, + { url = "https://files.pythonhosted.org/packages/84/44/f4d18446eacb20ea11e82f133ea8f86e2bf2891785b67d9da8d0ab0ef525/ml_dtypes-0.5.4-cp314-cp314t-win_amd64.whl", hash = "sha256:4381fe2f2452a2d7589689693d3162e876b3ddb0a832cde7a414f8e1adf7eab1", size = 236612, upload-time = "2025-11-17T22:32:22.579Z" }, + { url = "https://files.pythonhosted.org/packages/ad/3f/3d42e9a78fe5edf792a83c074b13b9b770092a4fbf3462872f4303135f09/ml_dtypes-0.5.4-cp314-cp314t-win_arm64.whl", hash = "sha256:11942cbf2cf92157db91e5022633c0d9474d4dfd813a909383bd23ce828a4b7d", size = 168825, upload-time = "2025-11-17T22:32:23.766Z" }, +] + +[[package]] +name = "numpy" +version = "2.5.1" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/22/fd/89965aa4ac08c74998539fcbf24fa3540f3e15237fbeb6bcf9c908f4aade/numpy-2.5.1.tar.gz", hash = "sha256:a48a113e6afea91f5608793bafa7ef2ad481fefbda87ec5069f483de61cb9fa3", size = 20755553, upload-time = "2026-07-04T17:08:00.933Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/eb/07/ec2a3f0c91761581d4b7104a740791800025983f9a4dc4e73f91a99aeac4/numpy-2.5.1-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:0bfebd8695f9863592fe744be833a258120b14a9f39da255e8aa8fade2c0ddd1", size = 16796419, upload-time = "2026-07-04T17:06:40.37Z" }, + { url = "https://files.pythonhosted.org/packages/ab/ab/ddb499fc4f8780354395face5b65c7fd107bcd6e1d667a5f07d046956f6f/numpy-2.5.1-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:30b44a6b53a7ae63c54c089a8726e5563ed302716c5b7ccc85afade40b0e7ff6", size = 11765832, upload-time = "2026-07-04T17:06:42.768Z" }, + { url = "https://files.pythonhosted.org/packages/88/b3/3c28c558a09fc72100c646dac6d2fce8e834c471b0edca01a29996706117/numpy-2.5.1-cp313-cp313-macosx_14_0_arm64.whl", hash = "sha256:6165343f81b56ef8f514f396989e529b61d9dc709b99421b07e9f3e698e2287d", size = 5325143, upload-time = "2026-07-04T17:06:45.466Z" }, + { url = "https://files.pythonhosted.org/packages/5e/0e/ce19b985bb15c596f4f05954e76cccc77c845083b3b8f938a6c68e523128/numpy-2.5.1-cp313-cp313-macosx_14_0_x86_64.whl", hash = "sha256:4939237038ada79308dda3204ac6462df056b5672b2e25db1149cf873668b3e1", size = 6659749, upload-time = "2026-07-04T17:06:47.288Z" }, + { url = "https://files.pythonhosted.org/packages/2e/20/1ee6614d64332a1bba6411f38e68cb79eec1b2459e20a623777c5c5492a2/numpy-2.5.1-cp313-cp313-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:1c6759f538fb912fc46de0a6b1758ccf7b57bc7c7ebebc23974fdac3de8db0cd", size = 15164716, upload-time = "2026-07-04T17:06:49.494Z" }, + { url = "https://files.pythonhosted.org/packages/ed/a7/2bcd3fdbb87804755c35b729bf8709d62025c5f4cfd7d5b2415997097515/numpy-2.5.1-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:9726558e8db4a5bf7929a70ae50f63abda4daf0efe810e3bfbab95976f75fc1a", size = 16661440, upload-time = "2026-07-04T17:06:52.061Z" }, + { url = "https://files.pythonhosted.org/packages/fc/d7/a41e3310c886fe457d36e670bbf24fae411aca8a7b6ad92a32afd924077c/numpy-2.5.1-cp313-cp313-musllinux_1_2_aarch64.whl", hash = "sha256:3935f3b419b244a02732676fa5317a9193cc596a4c0646db07e5b421229ac9f7", size = 16526305, upload-time = "2026-07-04T17:06:54.605Z" }, + { url = "https://files.pythonhosted.org/packages/53/75/4333a9a707c1edd3a4e1a0c58eca52c0f31e55089fa80db02b5565b24df7/numpy-2.5.1-cp313-cp313-musllinux_1_2_x86_64.whl", hash = "sha256:dc932a65ded7ce9013d120845a2514dcccb1a67bfc8deb8d37633762951904a6", size = 18423008, upload-time = "2026-07-04T17:06:57.54Z" }, + { url = "https://files.pythonhosted.org/packages/ee/90/e314a32b1c11a2ffe818ddad3a57b50b4b6e1b6c487192eb50cdef0415d0/numpy-2.5.1-cp313-cp313-win32.whl", hash = "sha256:4b4ff1608417eb7a59da7b967bbb798cacfe071d2caf526a24281cd562072ed9", size = 6063885, upload-time = "2026-07-04T17:07:00.14Z" }, + { url = "https://files.pythonhosted.org/packages/10/70/800b3fca480af32df9e8ea9f3d4a0c8feb4b32d7f195d174eabbda4829ad/numpy-2.5.1-cp313-cp313-win_amd64.whl", hash = "sha256:6c3fe51bc6a16453d452997053454f309e8e0ed7b42d6b361ce4ac8c32913d74", size = 12425674, upload-time = "2026-07-04T17:07:02.387Z" }, + { url = "https://files.pythonhosted.org/packages/8b/0b/196350c122f50f6ca56846f2d71efd5e0d24b7b2e07355e019b2e2c7a11e/numpy-2.5.1-cp313-cp313-win_arm64.whl", hash = "sha256:f7feb014281029e628ba2d5a007407443b06e418b6fe451d1e2adcbc8eba0107", size = 10350256, upload-time = "2026-07-04T17:07:04.878Z" }, + { url = "https://files.pythonhosted.org/packages/db/f4/731b6085a83faf6ca843394cbd5e217280c214399f7e8b21b9f552af0ae2/numpy-2.5.1-cp314-cp314-macosx_10_15_x86_64.whl", hash = "sha256:7c786fe9a5bbe360022e584c5a34cf6b54265c71bd7ec8ac3d8fec38968071f8", size = 16795063, upload-time = "2026-07-04T17:07:07.374Z" }, + { url = "https://files.pythonhosted.org/packages/bf/64/0e215f2048dd11a55bb989ed41b3585ef57452404e638d703a211a3e4157/numpy-2.5.1-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:32985c896d897419ef8da6917872d80b78ad0ea26d85b23245c7366ffde76d75", size = 11776652, upload-time = "2026-07-04T17:07:09.907Z" }, + { url = "https://files.pythonhosted.org/packages/b5/59/2b844c7a6e9deff69b404a66221e1542937734f65d5e6e39411876053862/numpy-2.5.1-cp314-cp314-macosx_14_0_arm64.whl", hash = "sha256:efd736408cc97c79b9e6917338dfc8f06013b2274f992e96b1d9a81a71e2a2c2", size = 5335944, upload-time = "2026-07-04T17:07:12.227Z" }, + { url = "https://files.pythonhosted.org/packages/86/51/9bf7cb2cabcebc9e017e4ec7e6322b378317a542c08b4cb68479c1efc716/numpy-2.5.1-cp314-cp314-macosx_14_0_x86_64.whl", hash = "sha256:ab84dc6b074fa881cae55bea94cc4f68e285181ba7f32497bf7dee6b1496165b", size = 6656266, upload-time = "2026-07-04T17:07:14.368Z" }, + { url = "https://files.pythonhosted.org/packages/83/3e/fb7615b211b82a32f44d5180a6d421b61f84d4fadd578b48ba4ac34e189f/numpy-2.5.1-cp314-cp314-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:caf3e317d33d60c37986b452613f4ab51246d0691350c03d0cb4a898627f4a95", size = 15179720, upload-time = "2026-07-04T17:07:16.272Z" }, + { url = "https://files.pythonhosted.org/packages/41/5f/0f992cb24560673496c5d68de61913b57166ce530ffda07c1f280e0cc464/numpy-2.5.1-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:54ad769f17bc2d833b620851989f62054fb9ab93c969d9e1dc3c8e3d56beea21", size = 16664835, upload-time = "2026-07-04T17:07:19.021Z" }, + { url = "https://files.pythonhosted.org/packages/a2/2f/97d6475ee91afe2587797d09446f9d3e475ad4cb681662d824809327b75a/numpy-2.5.1-cp314-cp314-musllinux_1_2_aarch64.whl", hash = "sha256:c12afb53450fa976d4c681c50a7423729a4c51c0465ed9f32b8a9cabbc472373", size = 16539135, upload-time = "2026-07-04T17:07:22.015Z" }, + { url = "https://files.pythonhosted.org/packages/c4/5b/4db81e4ba0be7e2776b1de68c82aa862c7f8ec27e1b4927d4ae075e20678/numpy-2.5.1-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:e8c11c405efc5ff6816d5983c96cdfa215bab3428961243af3ff59b228490438", size = 18426684, upload-time = "2026-07-04T17:07:24.941Z" }, + { url = "https://files.pythonhosted.org/packages/1f/64/c0ba2d90724d450279a7df8f32057241070250a26a7e2b5337d77347f481/numpy-2.5.1-cp314-cp314-win32.whl", hash = "sha256:f2479a47f8d5932d1718168a681ad6e536a9df484c83cfcf9de365e164537ace", size = 6116103, upload-time = "2026-07-04T17:07:27.622Z" }, + { url = "https://files.pythonhosted.org/packages/c1/1a/837f9ed7405adcd7a40538792eb169eddd8fa5630c16a1ef49dae71a30f4/numpy-2.5.1-cp314-cp314-win_amd64.whl", hash = "sha256:24d0eb82c0541d3415a33425db64ae439dffccd7b4dbcb30e7c35120205c506a", size = 12562177, upload-time = "2026-07-04T17:07:29.887Z" }, + { url = "https://files.pythonhosted.org/packages/22/ed/49707938b6dd0a78a9178dd93227dc89e4c11af47f5c798d70366e8d0483/numpy-2.5.1-cp314-cp314-win_arm64.whl", hash = "sha256:5a4c988b38d261deeeaad9954e3deb091ad905c94e8bb6708654ef1d97f286b0", size = 10627739, upload-time = "2026-07-04T17:07:32.568Z" }, + { url = "https://files.pythonhosted.org/packages/a6/c7/bb4b882cfe7f299cbc8b66e42e7dd78cf9d14e40f9469fc5e3db7e15b3bd/numpy-2.5.1-cp314-cp314t-macosx_11_0_arm64.whl", hash = "sha256:a33276be12fa045805f477f22482088b66bb758ffbe89a9d21457de863a32e22", size = 11894709, upload-time = "2026-07-04T17:07:34.941Z" }, + { url = "https://files.pythonhosted.org/packages/40/3f/5af7f4a7f6224aef48017aa82bb6174c7a659d724be0c75017b7e64a55b4/numpy-2.5.1-cp314-cp314t-macosx_14_0_arm64.whl", hash = "sha256:f089d7b00756190aacf1f5d34bdf38c3c430ac82b4f868f8cede73380460fce7", size = 5453810, upload-time = "2026-07-04T17:07:37.495Z" }, + { url = "https://files.pythonhosted.org/packages/20/c9/3474309bc94d634d3f9c3eddf03250ecb8c22cd948ef16fef69a77cc5d7b/numpy-2.5.1-cp314-cp314t-macosx_14_0_x86_64.whl", hash = "sha256:09e9bfd8d2cf479c7d174804fb3811c53a8e9f20a37444008606b57d6b7a826d", size = 6761189, upload-time = "2026-07-04T17:07:39.563Z" }, + { url = "https://files.pythonhosted.org/packages/90/8a/558ae39fdd55d7e7f7fef9a84a6e964ac6b23edbd2a07e52bb084500507d/numpy-2.5.1-cp314-cp314t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:e68d8dd1e7eba712948f2053a29ec86917bc70ba1358df869d9f06649ef9cf09", size = 15225039, upload-time = "2026-07-04T17:07:41.682Z" }, + { url = "https://files.pythonhosted.org/packages/63/27/ca7392b2d030277bdf0273e7d23255b3ee57d57a7c170a6f4fb3981e1e5d/numpy-2.5.1-cp314-cp314t-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:99d5095fa265a0c4152e7bb12759e14381ef5496152f1ce58f44bdf55c44beb4", size = 16701306, upload-time = "2026-07-04T17:07:44.611Z" }, + { url = "https://files.pythonhosted.org/packages/02/42/03d53ae7996c44d4374a8262e9dc41671fd56cbb98f7d47ef85cf5da4c6b/numpy-2.5.1-cp314-cp314t-musllinux_1_2_aarch64.whl", hash = "sha256:ab87a91b3cc3382b8956095bd8f95e00cf679bb81554339be1a2ba404a1473c1", size = 16589955, upload-time = "2026-07-04T17:07:47.694Z" }, + { url = "https://files.pythonhosted.org/packages/7b/15/6c1784ae469640e65db111e9a34b3d0f14d91e8a38b9ce34810ced370dbb/numpy-2.5.1-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:224ca51130ef7da85bea2191625181cb4f337f9cb64b471f10c1a12aa8b60077", size = 18464252, upload-time = "2026-07-04T17:07:50.684Z" }, + { url = "https://files.pythonhosted.org/packages/94/a8/f98e50356cf167df656c526c2dfeec2d7dde182f2a3da4b458a5938e2776/numpy-2.5.1-cp314-cp314t-win32.whl", hash = "sha256:6eab239876581b2b3c5a242281b6007bbdbcd1c7085d7709bb57c5929b11e6bf", size = 6263298, upload-time = "2026-07-04T17:07:53.445Z" }, + { url = "https://files.pythonhosted.org/packages/72/ac/96ae880cdecad0b3275d9359fcec72667b49a4863c9f12942e43679dda02/numpy-2.5.1-cp314-cp314t-win_amd64.whl", hash = "sha256:83ce9c80d5b521b0d77ddcbe5447c218d247929b6cc056ca5351342accfff0af", size = 12748623, upload-time = "2026-07-04T17:07:55.384Z" }, + { url = "https://files.pythonhosted.org/packages/a1/5a/4d2b1601df3602dba7a14f3348ba9bfe94a18adb428e693df6154c293831/numpy-2.5.1-cp314-cp314t-win_arm64.whl", hash = "sha256:5a6db61f9aaa57e369905c67d852045d3c4f7126405b29d09b19dec118e9c9cb", size = 10697674, upload-time = "2026-07-04T17:07:58.506Z" }, +] + +[[package]] +name = "opt-einsum" +version = "3.4.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/8c/b9/2ac072041e899a52f20cf9510850ff58295003aa75525e58343591b0cbfb/opt_einsum-3.4.0.tar.gz", hash = "sha256:96ca72f1b886d148241348783498194c577fa30a8faac108586b14f1ba4473ac", size = 63004, upload-time = "2024-09-26T14:33:24.483Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/23/cd/066e86230ae37ed0be70aae89aabf03ca8d9f39c8aea0dec8029455b5540/opt_einsum-3.4.0-py3-none-any.whl", hash = "sha256:69bb92469f86a1565195ece4ac0323943e83477171b91d24c35afe028a90d7cd", size = 71932, upload-time = "2024-09-26T14:33:23.039Z" }, +] + +[[package]] +name = "packaging" +version = "26.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d7/f1/e7a6dd94a8d4a5626c03e4e99c87f241ba9e350cd9e6d75123f992427270/packaging-26.2.tar.gz", hash = "sha256:ff452ff5a3e828ce110190feff1178bb1f2ea2281fa2075aadb987c2fb221661", size = 228134, upload-time = "2026-04-24T20:15:23.917Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/df/b2/87e62e8c3e2f4b32e5fe99e0b86d576da1312593b39f47d8ceef365e95ed/packaging-26.2-py3-none-any.whl", hash = "sha256:5fc45236b9446107ff2415ce77c807cee2862cb6fac22b8a73826d0693b0980e", size = 100195, upload-time = "2026-04-24T20:15:22.081Z" }, +] + +[[package]] +name = "pluggy" +version = "1.6.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/f9/e2/3e91f31a7d2b083fe6ef3fa267035b518369d9511ffab804f839851d2779/pluggy-1.6.0.tar.gz", hash = "sha256:7dcc130b76258d33b90f61b658791dede3486c3e6bfb003ee5c9bfb396dd22f3", size = 69412, upload-time = "2025-05-15T12:30:07.975Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/54/20/4d324d65cc6d9205fabedc306948156824eb9f0ee1633355a8f7ec5c66bf/pluggy-1.6.0-py3-none-any.whl", hash = "sha256:e920276dd6813095e9377c0bc5566d94c932c33b27a3e3945d8389c374dd4746", size = 20538, upload-time = "2025-05-15T12:30:06.134Z" }, +] + +[[package]] +name = "pygments" +version = "2.20.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/c3/b2/bc9c9196916376152d655522fdcebac55e66de6603a76a02bca1b6414f6c/pygments-2.20.0.tar.gz", hash = "sha256:6757cd03768053ff99f3039c1a36d6c0aa0b263438fcab17520b30a303a82b5f", size = 4955991, upload-time = "2026-03-29T13:29:33.898Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/f4/7e/a72dd26f3b0f4f2bf1dd8923c85f7ceb43172af56d63c7383eb62b332364/pygments-2.20.0-py3-none-any.whl", hash = "sha256:81a9e26dd42fd28a23a2d169d86d7ac03b46e2f8b59ed4698fb4785f946d0176", size = 1231151, upload-time = "2026-03-29T13:29:30.038Z" }, +] + +[[package]] +name = "pytest" +version = "9.1.1" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "colorama", marker = "sys_platform == 'win32'" }, + { name = "iniconfig" }, + { name = "packaging" }, + { name = "pluggy" }, + { name = "pygments" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/e4/47/b9efed96c114afcfa3c9d3fe98a76a1d14c74a9e266d397cf6eb64be5e01/pytest-9.1.1.tar.gz", hash = "sha256:1088fbde8f2b49d95a549a195707afa7a76a3ce9bcadc26b6d71f0ffda5fe313", size = 1636369, upload-time = "2026-06-19T10:58:32.857Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/24/25/1de2678b631f5a49215c6c96fff41ba892b0a34df68d6d80292b1b48aa7f/pytest-9.1.1-py3-none-any.whl", hash = "sha256:37a86b45efb9a47a61a36449063e8e18d0cab3161329fc099eb21783169c4f0c", size = 386536, upload-time = "2026-06-19T10:58:31.347Z" }, +] + +[[package]] +name = "scipy" +version = "1.18.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/a7/25/c2700dfaf6442b4effaa91af24ebce5dc9d31bb4a69706313aae70d72cd0/scipy-1.18.0.tar.gz", hash = "sha256:67b2ad2ad54c72ca6d04975a9b2df8c3638c34ddd5b28738e94fc2b57929d378", size = 30774447, upload-time = "2026-06-19T15:01:43.456Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/05/52/9c0136c2de7ae0779b7b366447766cec6d9f0702c56bb8ffeb04c8fd3af4/scipy-1.18.0-cp313-cp313-macosx_10_15_x86_64.whl", hash = "sha256:09143f676d157d9f546d663504ef9c1becb819824f1afc018814176411942446", size = 31036107, upload-time = "2026-06-19T15:00:14.03Z" }, + { url = "https://files.pythonhosted.org/packages/02/73/0291a64843270f4efb86cdcf2ee0f2048631b65ec6b405398b2b4dbf11bf/scipy-1.18.0-cp313-cp313-macosx_12_0_arm64.whl", hash = "sha256:5efe260f69417b97ddae455bfb5a95e8359f7f66ad7fa9522a60feb66f169520", size = 28663303, upload-time = "2026-06-19T15:00:16.819Z" }, + { url = "https://files.pythonhosted.org/packages/d3/0f/10ffa0b697a572f4e0d48b92a88895d366422f019f723e7e14a84c050dac/scipy-1.18.0-cp313-cp313-macosx_14_0_arm64.whl", hash = "sha256:68363b7eaacd8b5dd426df56d782cc156468ac79a127a1b87ca597d6e2e82197", size = 20404960, upload-time = "2026-06-19T15:00:19.635Z" }, + { url = "https://files.pythonhosted.org/packages/7e/d2/e896cea21ba8edd6c81d4c55b1ffcc717e79698dcbebf9641b4cfb4c6622/scipy-1.18.0-cp313-cp313-macosx_14_0_x86_64.whl", hash = "sha256:c5557d8be5da8e41353fcd4d21491fdbab83b062fc579e94dc09a7c8ab4f669b", size = 23034074, upload-time = "2026-06-19T15:00:22.107Z" }, + { url = "https://files.pythonhosted.org/packages/ea/b2/e83ea34279a52c03374477c74006256ec78df65fc877baa4617d6de1d202/scipy-1.18.0-cp313-cp313-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:0d13bca67c096d89fb95ced0d8921807300fce0275643aef9533cc63a0773468", size = 33942038, upload-time = "2026-06-19T15:00:24.964Z" }, + { url = "https://files.pythonhosted.org/packages/f6/af/e8fe5fb136f51e2b01678b92cb4106d10d8cd68ec147ead2e7cb0ac75398/scipy-1.18.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:a46f9273dbd0eb1cefba61c9b8648b4dfe3cbc14a080176f9a73e44b8336dc7f", size = 35266390, upload-time = "2026-06-19T15:00:28.059Z" }, + { url = "https://files.pythonhosted.org/packages/3a/49/2c5cbb907b56695fc67517811d1db234dfd83381a84814ec220aded2794d/scipy-1.18.0-cp313-cp313-musllinux_1_2_aarch64.whl", hash = "sha256:5aba46108853ddfc77906b6557aac839d2b52e900c1d72a1180adaaab58d265f", size = 35551324, upload-time = "2026-06-19T15:00:31.014Z" }, + { url = "https://files.pythonhosted.org/packages/bb/73/eda39f7a2d306ff0ffc574afd13c0bbb6d10a603d9a413998ee269487a80/scipy-1.18.0-cp313-cp313-musllinux_1_2_x86_64.whl", hash = "sha256:b6f758e35f12757b5d95c00bc6de2438e229c2664b7a92e96f205959d9f2dfa4", size = 37404785, upload-time = "2026-06-19T15:00:34.072Z" }, + { url = "https://files.pythonhosted.org/packages/b7/d2/ae881ee28d014f38e0ccbfd974a06a919ba9af34f1f74bf42b5301891d63/scipy-1.18.0-cp313-cp313-win_amd64.whl", hash = "sha256:1afac4a847207c7ff8efd321734a50b06d0280b3b2a2c0fc2f413101747ad7c7", size = 36554943, upload-time = "2026-06-19T15:00:36.903Z" }, + { url = "https://files.pythonhosted.org/packages/70/3a/21154e2d54eb3639c6bf4dbae2e531c68356bfe95990daa30df33b30d556/scipy-1.18.0-cp313-cp313-win_arm64.whl", hash = "sha256:c5dbddf60e58c2312316d097271a8e73d40eaf2eabfa4d95ed7d3695bbf2ce7b", size = 24350911, upload-time = "2026-06-19T15:00:40.062Z" }, + { url = "https://files.pythonhosted.org/packages/78/b5/915a19b3de2f7430062b509653563db1633ddbb6f021b06731521115d4e2/scipy-1.18.0-cp314-cp314-macosx_10_15_x86_64.whl", hash = "sha256:4c256ee70c0d1a8a2ace807e199ccd4e3f57037433842abb3fb36bc17eaa9578", size = 31036253, upload-time = "2026-06-19T15:00:43.216Z" }, + { url = "https://files.pythonhosted.org/packages/d7/88/b72def7262e150d16be13fca37a96481138d624e700340bc3362a7588929/scipy-1.18.0-cp314-cp314-macosx_12_0_arm64.whl", hash = "sha256:2ef3abc54a4ffc53765374b0d5728532dfdd2585ed23f6b11c206a1f0b1b9af8", size = 28673758, upload-time = "2026-06-19T15:00:46.663Z" }, + { url = "https://files.pythonhosted.org/packages/91/02/2e636a61a525632c373cf6a9c24442a3ffb79e364d38e98b32042964ac32/scipy-1.18.0-cp314-cp314-macosx_14_0_arm64.whl", hash = "sha256:f2a6af57bd9e4a75d70e4117e78a1bbee84f79ae3fbb6d0111005d6ebcc4cb8d", size = 20415514, upload-time = "2026-06-19T15:00:49.399Z" }, + { url = "https://files.pythonhosted.org/packages/c9/b6/2135974442f6aba159d9d39d774a1c8cb19947016725d69fecc685df45bf/scipy-1.18.0-cp314-cp314-macosx_14_0_x86_64.whl", hash = "sha256:3f1ac564d3bf6c03d861d2cd87a1bea0da2887136f7fb1bf519c05a8971452d6", size = 23034398, upload-time = "2026-06-19T15:00:51.941Z" }, + { url = "https://files.pythonhosted.org/packages/f6/e6/ba89ec5abf6ee9257c0d1ec985573f3ae32742c24bc03e016388a40b1b15/scipy-1.18.0-cp314-cp314-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:40395a5fcd1abee49a5c7aaa98c29db393eedc835138560a588c47ec16156690", size = 33998032, upload-time = "2026-06-19T15:00:54.838Z" }, + { url = "https://files.pythonhosted.org/packages/7f/c4/bc41eb19b0fd0db868f4132920879019318d80cc522ad8f2bca4611af808/scipy-1.18.0-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:8ca01e8ae69f1b18e9a58d91afead31be3cef0dd905a10249dac559ee15460a0", size = 35283333, upload-time = "2026-06-19T15:00:58.152Z" }, + { url = "https://files.pythonhosted.org/packages/53/a4/cbdeef6eb3830a8462a9d4ada814de5fc984345cc9ecf17cbec51a036f1e/scipy-1.18.0-cp314-cp314-musllinux_1_2_aarch64.whl", hash = "sha256:7a7f3b01647384dbc3a711e8c6778e0aabbe93959249fef5c7393396bcac0867", size = 35610216, upload-time = "2026-06-19T15:01:01.155Z" }, + { url = "https://files.pythonhosted.org/packages/80/4d/b2b82502b65f661d1b789c1665dcdf315d5f12194e06fc0b37946294ebae/scipy-1.18.0-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:6aa94e78ec192a30063a5e72e561c28af769dc311190b24fe91774eff1969709", size = 37418960, upload-time = "2026-06-19T15:01:04.155Z" }, + { url = "https://files.pythonhosted.org/packages/93/3e/902d836831474b0ab5a37d16404f7bc5fafd9efba632890e271ba952635f/scipy-1.18.0-cp314-cp314-win_amd64.whl", hash = "sha256:2d8bbdc6c817f5b4006a54d799d4f5bab6f910193cbb9a1ff310833d4d270f61", size = 37288845, upload-time = "2026-06-19T15:01:07.822Z" }, + { url = "https://files.pythonhosted.org/packages/b6/43/8d73b337a3bdb14daa0314f0434210747c02d79d729ce1777574a817dcf6/scipy-1.18.0-cp314-cp314-win_arm64.whl", hash = "sha256:18e9575f1569b2c54174e6159d32942e03731177f63dce7975f0a0c88d102f5b", size = 24988971, upload-time = "2026-06-19T15:01:11.076Z" }, + { url = "https://files.pythonhosted.org/packages/b4/b4/f11918b0508a2787031a0499a03fbe3546f3bb5ca05d01038c45b278c09a/scipy-1.18.0-cp314-cp314t-macosx_10_15_x86_64.whl", hash = "sha256:f351e0dd702687d12a402b867a1b4146a256923e1c38317cbc472f6372b94707", size = 31399325, upload-time = "2026-06-19T15:01:13.723Z" }, + { url = "https://files.pythonhosted.org/packages/7b/d1/1f287b57c0ff0ee5185dff3946d92c8017d39b0e431f0ae79a3ff1859512/scipy-1.18.0-cp314-cp314t-macosx_12_0_arm64.whl", hash = "sha256:7c7a51b33ce387193c97f228320cf8e87361daa1bba750638677729598b3e677", size = 29092110, upload-time = "2026-06-19T15:01:16.908Z" }, + { url = "https://files.pythonhosted.org/packages/ff/1a/7b74eb6c392fdcb27d414c0e7558a6d0231eb3b6d73571f479bb81ea8794/scipy-1.18.0-cp314-cp314t-macosx_14_0_arm64.whl", hash = "sha256:84031d7b052a54fae2f8632e0ec802073d385476eb9a63079bce6e23ef9283d4", size = 20833811, upload-time = "2026-06-19T15:01:20.488Z" }, + { url = "https://files.pythonhosted.org/packages/7c/ad/f3941716320a7b9cb4d68734a903b45fe16eff5fb7da7e16f2e619304979/scipy-1.18.0-cp314-cp314t-macosx_14_0_x86_64.whl", hash = "sha256:56abf29a7c067dde59be8b9a22d606a4ea1b2f2a4b756d9d903c62818f5dacce", size = 23396644, upload-time = "2026-06-19T15:01:23.364Z" }, + { url = "https://files.pythonhosted.org/packages/22/22/1446b62ffe07f9719b7d9b1b6a4e05a772833ae8f441fe4c22c34c9b250f/scipy-1.18.0-cp314-cp314t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:1ad44305cfa24b1ba5803cbbebf033590ccbac1aa5d612d727b785325ab408b0", size = 34079318, upload-time = "2026-06-19T15:01:26.002Z" }, + { url = "https://files.pythonhosted.org/packages/56/3b/b87da667098bb470fa30c7011b0ba351ee976dd395c78798c66e941665a3/scipy-1.18.0-cp314-cp314t-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:945c1761b93f38d7f99ae81ae80c63e621471608c7eeead563f6df025585cd58", size = 35324320, upload-time = "2026-06-19T15:01:28.881Z" }, + { url = "https://files.pythonhosted.org/packages/f8/a1/c7932f91909759b0267f75fdea34e91309f96b895757534b76a90b6b4344/scipy-1.18.0-cp314-cp314t-musllinux_1_2_aarch64.whl", hash = "sha256:1a4441f15d620578772a49e5ab48c0ee1f7a0220e387110283062729136b2553", size = 35699541, upload-time = "2026-06-19T15:01:31.968Z" }, + { url = "https://files.pythonhosted.org/packages/f7/86/5185061a1fcc41d18c5dc2463969b3a3964b31d9ac67b2fb05d4c7ff7670/scipy-1.18.0-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:9aac6192fac56bf2ca534389d24623f07b39ff83317d58287285e7fbd622ff76", size = 37472480, upload-time = "2026-06-19T15:01:35.136Z" }, + { url = "https://files.pythonhosted.org/packages/31/8e/f04c68e39919a010d34f2ee1367fd705b0a25a02f609d755f0bfbc0a15fc/scipy-1.18.0-cp314-cp314t-win_amd64.whl", hash = "sha256:e40baea28ae7f5475c779741e2d90b1247c78531207b49c7030e698ff81cee3f", size = 37365390, upload-time = "2026-06-19T15:01:38.091Z" }, + { url = "https://files.pythonhosted.org/packages/d5/19/969dc072906c84dd0a3b05dcf57ea750936087d7873549e408b35cfc3f97/scipy-1.18.0-cp314-cp314t-win_arm64.whl", hash = "sha256:368e0a705903c466aa5f08eefb39e6b1b6b2d659e7352a31fd9e2438365be0f8", size = 25279661, upload-time = "2026-06-19T15:01:40.817Z" }, +]