From ddef158ad6fde21523e156d96427005527eda21e Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 17 Aug 2026 09:24:11 -0500 Subject: [PATCH 1/2] Add Sphinx background documentation for AO integrals Adds the docs/source/background tree covering overlap, kinetic, nuclear attraction, and electron repulsion integrals, including the density-fitting derivation. Wires up sphinxcontrib-bibtex with a references.bib, and excludes that file from the license header check. Co-Authored-By: Claude Opus 5 --- .licenserc.yaml | 1 + docs/README.md | 65 +++ docs/requirements.txt | 1 + .../background/eris/density_fitting.rst | 219 +++++++++ docs/source/background/eris/eris.rst | 38 ++ docs/source/background/eris/error_bounds.rst | 461 ++++++++++++++++++ docs/source/background/eris/index.rst | 24 + docs/source/background/eris/screening.rst | 316 ++++++++++++ docs/source/background/index.rst | 26 + docs/source/background/kinetic.rst | 28 ++ docs/source/background/nuclear_attraction.rst | 29 ++ docs/source/background/overlap.rst | 26 + docs/source/conf.py | 5 + docs/source/index.rst | 1 + docs/source/references.bib | 68 +++ 15 files changed, 1308 insertions(+) create mode 100644 docs/README.md create mode 100644 docs/source/background/eris/density_fitting.rst create mode 100644 docs/source/background/eris/eris.rst create mode 100644 docs/source/background/eris/error_bounds.rst create mode 100644 docs/source/background/eris/index.rst create mode 100644 docs/source/background/eris/screening.rst create mode 100644 docs/source/background/index.rst create mode 100644 docs/source/background/kinetic.rst create mode 100644 docs/source/background/nuclear_attraction.rst create mode 100644 docs/source/background/overlap.rst create mode 100644 docs/source/references.bib diff --git a/.licenserc.yaml b/.licenserc.yaml index 87c2fcca..11a9d8ae 100644 --- a/.licenserc.yaml +++ b/.licenserc.yaml @@ -24,5 +24,6 @@ header: - docs/requirements.txt - version.txt - build/ + - docs/source/references.bib comment: never diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 00000000..618d3178 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,65 @@ + + +Building the Integrals Documentation +===================================== + +This directory contains the source files for generating the Sphinx +documentation for `Integrals`. General instructions for building +documentation found throughout the NWChemEx project are available at: +https://github.com/NWChemEx/NWChemEx/blob/master/docs/README.md + +Obtaining the Documentation's Dependencies +------------------------------------------- + +The documentation's dependencies can be installed via Python's `pip` +command. Commands are assumed to be run from this directory (the same +directory as this README file). + +~~~.sh +# These first two steps are strongly recommended, but not required +python3 -m venv venv +. venv/bin/activate +pip3 install -r requirements.txt +~~~ + +Building the Documentation +---------------------------- + +With the dependencies installed, build the documentation with: + +~~~.sh +make html BUILDDIR=${BUILDDIR} +~~~ + +where `${BUILDDIR}` is the directory where you want the resulting HTML to be +placed (*e.g.* `build`). + +Viewing the Documentation Locally +------------------------------------ + +After building, the main index of the resulting HTML will be located at +`${BUILDDIR}/html/index.html` and can be viewed by pointing your web browser +of choice at that file, either by opening it directly or by using +`file:///path/to/${BUILDDIR}/html/index.html` as the URL. Alternatively, for +a locally served copy that resembles how the docs are hosted online, run a +simple HTTP server from the build directory: + +~~~.sh +python3 -m http.server --directory ${BUILDDIR}/html +~~~ + +and then navigate to `http://localhost:8000` in your browser. diff --git a/docs/requirements.txt b/docs/requirements.txt index 528ccda4..fc8972f1 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,4 @@ GitPython sphinx==v7.2.6 sphinx_rtd_theme==1.3.0 +sphinxcontrib-bibtex==2.5.0 diff --git a/docs/source/background/eris/density_fitting.rst b/docs/source/background/eris/density_fitting.rst new file mode 100644 index 00000000..d739b8ed --- /dev/null +++ b/docs/source/background/eris/density_fitting.rst @@ -0,0 +1,219 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +############### +Density Fitting +############### + +.. |eri4| replace:: :math:`\left(\mu\nu|\lambda\sigma\right)` +.. |eri3| replace:: :math:`\left(P|\kappa\lambda\right)` +.. |eri2| replace:: :math:`\left(P|Q\right)` +.. |Bcal| replace:: :math:`B_{P,\mu\nu}` + +Density fitting (DF), also known as the resolution-of-the-identity (RI) +approximation, reduces the cost of computing the four-center ERI |eri4| +(see :doc:`eris` for its definition and notation) by approximately +factorizing it into a product of three-index quantities. This page derives +the DF approximation using the Coulomb metric, which is the variant +implemented in the ``Integrals`` package (see +:ref:`df-implementation-notes` below). + +************************** +The Density-Fitting Ansatz +************************** + +DF/RI introduces an auxiliary basis set, +:math:`\lbrace\chi_P\left(\vec{r}\right)\rbrace_{P=1}^{N_{aux}}`, and +approximates each AO charge distribution (product density) +:math:`\rho_{\mu\nu}\left(\vec{r}\right) \equiv +\phi_\mu\left(\vec{r}\right)\phi_\nu\left(\vec{r}\right)` as a linear +expansion in that auxiliary basis :cite:`whitten1973,dunlap1979`: + +.. math:: + + \rho_{\mu\nu}\left(\vec{r}\right) \approx + \tilde{\rho}_{\mu\nu}\left(\vec{r}\right) = + \sum_P^{N_{aux}} c^P_{\mu\nu}\chi_P\left(\vec{r}\right). + +Typically :math:`N_{aux}` is chosen to be a modest multiple of :math:`N_b` +(rather than scaling as :math:`N_b^2`, the number of AO pairs), which is +what makes the approximation useful. + +************************************* +Fitting Coefficients: Coulomb Metric +************************************* + +The expansion coefficients :math:`c^P_{\mu\nu}` are determined +variationally, by minimizing the residual self-repulsion energy of the +error density :math:`\Delta\rho_{\mu\nu} = +\rho_{\mu\nu} - \tilde\rho_{\mu\nu}` under the Coulomb metric +:cite:`whitten1973,dunlap1979`: + +.. math:: + + E^{err}_{\mu\nu} = \int\int + \Delta\rho_{\mu\nu}\left(\vec{r}_1\right) + \frac{1}{r_{12}} + \Delta\rho_{\mu\nu}\left(\vec{r}_2\right) + d\vec{r}_1 d\vec{r}_2 \geq 0. + +Defining the two-center Coulomb-metric matrix over the auxiliary basis, + +.. math:: + + \left(P|Q\right) = \int\int + \chi_P\left(\vec{r}_1\right)\frac{1}{r_{12}} + \chi_Q\left(\vec{r}_2\right) d\vec{r}_1 d\vec{r}_2, + +and the three-center integral between the auxiliary basis and an AO pair, + +.. math:: + + \left(P|\mu\nu\right) = \int\int + \chi_P\left(\vec{r}_1\right)\frac{1}{r_{12}} + \rho_{\mu\nu}\left(\vec{r}_2\right) d\vec{r}_1 d\vec{r}_2, + +setting :math:`\partial E^{err}_{\mu\nu} / \partial c^P_{\mu\nu} = 0` gives +the well-known result: + +.. math:: + :label: fit-coeffs + + c^P_{\mu\nu} = \sum_Q^{N_{aux}} \left(P|Q\right)^{-1}\left(Q|\mu\nu\right). + +Substituting the fitted densities back into the definition of the +four-center ERI, Eq. :eq:`eri4`, yields the DF approximation +:cite:`vahtras1993,feyereisen1993`: + +.. math:: + :label: df-eri4 + + \left(\mu\nu|\lambda\sigma\right) \approx + \sum_P^{N_{aux}}\sum_Q^{N_{aux}} + \left(\mu\nu|P\right)\left(P|Q\right)^{-1}\left(Q|\lambda\sigma\right). + +Because the metric matrix :math:`\left(P|Q\right)` is symmetric +positive-definite, Eq. :eq:`df-eri4` can always be written as a contraction +of two three-index tensors, + +.. math:: + :label: df-three-index + + \left(\mu\nu|\lambda\sigma\right) \approx + \sum_P^{N_{aux}} B_{P,\mu\nu} B_{P,\lambda\sigma}, + +for *any* factorization :math:`\left(P|Q\right)^{-1} = \sum_R M_{PR}M_{QR}` +of the inverse metric, with :math:`B_{P,\mu\nu} = \sum_Q M_{PQ}\left(Q|\mu\nu\right)`. +The two common choices for :math:`M` are the symmetric inverse square root, +:math:`M = \left(P|Q\right)^{-1/2}`, and a Cholesky factor of the *inverse* +metric; the latter is used in this package and is discussed next. + +.. _df-implementation-notes: + +******************************************* +Implementation: Cholesky-Factored Metric +******************************************* + +Rather than diagonalizing :math:`\left(P|Q\right)` to form the symmetric +inverse square root, the ``Integrals`` package Cholesky-factorizes the +Coulomb metric, + +.. math:: + + \left(P|Q\right) = \sum_R L_{PR}L_{QR} = \left(LL^T\right)_{PQ}, + +with :math:`L` lower triangular, and forms :math:`B_{P,\mu\nu}` using +:math:`M = L^{-1}`: + +.. math:: + :label: cholesky-fit + + B_{P,\mu\nu} = \sum_Q \left(L^{-1}\right)_{PQ}\left(Q|\mu\nu\right). + +This is equivalent to Eq. :eq:`df-three-index` since +:math:`\sum_R\left(L^{-1}\right)_{PR}\left(L^{-1}\right)_{QR} = +\left[\left(LL^T\right)^{-1}\right]_{PQ} = \left(P|Q\right)^{-1}`, but avoids +an explicit eigendecomposition of the metric; the Cholesky-decomposition +route to the fitted three-index tensor is discussed and benchmarked against +the symmetric-inverse-square-root (RI) route by +:cite:t:`weigend2009`. Concretely: + +- The ``CoulombMetric`` module computes the two-center integral |eri2| (via + its "Two-center ERI" submodule, which defaults to the raw ``ERI2`` + integral driver), Cholesky-decomposes it, and returns :math:`M = L^{-1}`. +- The ``DFIntegral`` module computes the raw three-center integral |eri3| + (via its "Three-center ERI" submodule, defaulting to the raw ``ERI3`` + driver) and its "Coulomb Metric" submodule (defaulting to + ``CoulombMetric``), then contracts them to form |Bcal| following + Eq. :eq:`cholesky-fit`. + +********************** +Density-Fitted J and K +********************** + +In self-consistent-field methods the Coulomb (J) and exchange (K) matrices +are built by contracting |eri4| with an AO density matrix :math:`D_{\kappa\lambda}`: + +.. math:: + + J_{\mu\nu} = \sum_{\kappa\lambda}\left(\mu\nu|\kappa\lambda\right)D_{\kappa\lambda}, + \qquad + K_{\mu\nu} = \sum_{\kappa\lambda}\left(\mu\lambda|\kappa\nu\right)D_{\kappa\lambda}. + +Substituting the DF approximation, Eq. :eq:`df-three-index`, for |eri4| +avoids ever forming the four-index quantities explicitly. + +**Density-Fitted J** (``JDensityFitted``). Contract the density into the +fitted tensor to form an intermediate over only the auxiliary index, + +.. math:: + + j_P = \sum_{\kappa\lambda} D_{\kappa\lambda}B_{P,\kappa\lambda}, + +then contract back onto the AO pair: + +.. math:: + + J_{\mu\nu} = \sum_P^{N_{aux}} j_P B_{P,\mu\nu}. + +This costs :math:`O\left(N_{aux}N_b^2\right)` to form :math:`j_P` and +:math:`O\left(N_{aux}N_b^2\right)` to form :math:`J_{\mu\nu}`, i.e. no step +scales worse than cubically in the system size. + +**Density-Fitted K** (``KDensityFitted``). Exchange cannot be reduced to a +single index-:math:`P` intermediate because the density matrix and the two +fitted tensors do not share a common pair of AO indices. Instead, an +intermediate retaining the auxiliary index and one AO index is formed, + +.. math:: + + K^{aux}_{P,\mu,\lambda} = \sum_\kappa D_{\kappa\lambda} B_{P,\mu\kappa}, + +and then contracted with the second fitted tensor: + +.. math:: + + K_{\mu\nu} = \sum_P^{N_{aux}}\sum_\lambda K^{aux}_{P,\mu,\lambda}B_{P,\nu\lambda}. + +Because the intermediate :math:`K^{aux}_{P,\mu,\lambda}` retains three +indices (rather than the single auxiliary index :math:`j_P` used for J), +forming density-fitted K is more expensive than density-fitted J, scaling +as :math:`O\left(N_{aux}N_b^3\right)`. + +********** +References +********** + +.. bibliography:: + :style: unsrt diff --git a/docs/source/background/eris/eris.rst b/docs/source/background/eris/eris.rst new file mode 100644 index 00000000..a2ff066a --- /dev/null +++ b/docs/source/background/eris/eris.rst @@ -0,0 +1,38 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +################ +Four-Center ERIs +################ + +.. |eri4| replace:: :math:`\left(\mu\nu|\lambda\sigma\right)` + +Electron repulsion integrals (ERIs) over a set of :math:`N_b` atomic orbitals +(AOs), :math:`\lbrace\phi_\mu\left(\vec{r}\right)\rbrace`, are the +four-index quantities: + +.. math:: + :label: eri4 + + \left(\mu\nu|\lambda\sigma\right) = \int\int + \phi_\mu\left(\vec{r}_1\right)\phi_\nu\left(\vec{r}_1\right) + \frac{1}{r_{12}} + \phi_\lambda\left(\vec{r}_2\right)\phi_\sigma\left(\vec{r}_2\right) + d\vec{r}_1 d\vec{r}_2. + +There are formally :math:`O\left(N_b^4\right)` unique |eri4| and, even with +permutational and screening-based reductions, computing, storing, and +contracting them is the dominant cost of most Hartree-Fock and +post-Hartree-Fock methods. :doc:`density_fitting` discusses an approximation +that avoids forming |eri4| explicitly. diff --git a/docs/source/background/eris/error_bounds.rst b/docs/source/background/eris/error_bounds.rst new file mode 100644 index 00000000..369a57a9 --- /dev/null +++ b/docs/source/background/eris/error_bounds.rst @@ -0,0 +1,461 @@ +.. Copyright 2026 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +############################# +Bounding the Screening Error +############################# + +.. |eri4| replace:: :math:`\left(\mu\nu|\lambda\sigma\right)` +.. |delta| replace:: :math:`\delta_{\mu\nu\lambda\sigma}` +.. |Bq| replace:: :math:`B_q` + +:doc:`screening` describes how contributions to |eri4| are discarded. This +page asks what the discarding cost. The answer developed here is a rigorous +two-sided interval, +:math:`\left[\delta_{\min},\delta_{\max}\right]`, that is guaranteed to +contain the error in every screened integral. + +The construction has two ingredients. A magnitude bound, taken from the +Cauchy-Schwarz inequality, brackets each discarded contribution symmetrically +about zero. A *sign oracle* then collapses that symmetric interval to a +half-interval for those contributions whose sign can be established a priori, +tightening the result and -- more usefully -- making it asymmetric. The second +ingredient is optional: the bound is valid with no sign information at all, +and improves monotonically as more signs are determined. + +***************************** +The Exact Screening Error +***************************** + +Let :math:`t` be the screening threshold and :math:`S\left(t\right)` the set +of primitive quartets that the gates of :doc:`screening` discard. Write a +quartet as :math:`q = \left(ijkl\right)` and its exact contribution to Eq. +:eq:`scr-prim-expansion` as + +.. math:: + :label: err-tq + + t_q = c_ic_jc_kc_l\left[ij|kl\right]. + +The screened integral is the sum of Eq. :eq:`err-tq` over the quartets that +survive. Subtracting the exact result, which sums over all quartets, leaves + +.. math:: + :label: err-delta + + \delta_{\mu\nu\lambda\sigma} \equiv + \left(\mu\nu|\lambda\sigma\right)^{t} - + \left(\mu\nu|\lambda\sigma\right)^{\rm exact} = + -\sum_{q\in S\left(t\right)} t_q. + +The sign convention in Eq. :eq:`err-delta` is *computed minus exact*, and it +is used without exception below. A positive |delta| means the screened value +is too large. The leading minus sign is easy to lose: discarding a positive +contribution makes the computed integral too *small*. + +***************************************** +A Bracket from Per-Quartet Intervals +***************************************** + +Suppose that for each discarded quartet we can produce an interval +:math:`t_q \in \left[\ell_q, u_q\right]`. Because Eq. :eq:`err-delta` is a +plain sum, the intervals add, and negating flips the endpoints: + +.. math:: + :label: err-generic-bracket + + \delta_{\mu\nu\lambda\sigma} \in + \left[-\sum_{q\in S\left(t\right)}u_q,\; + -\sum_{q\in S\left(t\right)}\ell_q\right]. + +This is the entire structure of the model. Everything that follows is a +matter of supplying :math:`\ell_q` and :math:`u_q`. Nothing in Eq. +:eq:`err-generic-bracket` requires the intervals to be tight, symmetric, or +sign-resolved; it requires only that each one genuinely contain its +:math:`t_q`. + +********************** +The Magnitude Bound +********************** + +The default interval comes from Cauchy-Schwarz, Eq. :eq:`scr-cauchy-schwarz`. +Applying it to Eq. :eq:`err-tq`, + +.. math:: + :label: err-bq + + \left|t_q\right| \leq B_q \equiv + \left|c_ic_jc_kc_l\right| + \sqrt{\left[ij|ij\right]}\sqrt{\left[kl|kl\right]}, + +so that :math:`t_q \in \left[-B_q, +B_q\right]` unconditionally. + +The choice of Cauchy-Schwarz here is deliberate and is the single most +important design decision on this page. As :ref:`scr-rigor` records, the +coarse and fine metrics that drive the screening gates are derived from the +s-type closed form and retain no angular prefactor; for :math:`l > 0` the true +integral may exceed them. An error model built on the fine metric therefore +produces intervals that the true error can escape, which defeats the purpose. +Equation :eq:`scr-cauchy-schwarz` is a rigorous bound at every angular +momentum, so Eq. :eq:`err-bq` cannot be violated. + +In practice :math:`\sqrt{\left[ij|ij\right]}` is evaluated per shell block +rather than per AO. Writing :math:`M_{ij}` for the block of +:math:`\left(ij|ij\right)` values, any quantity at least as large as the +largest diagonal element of :math:`M_{ij}` may be substituted; the Frobenius +norm :math:`\left\|M_{ij}\right\|_F` computed by +``CauchySchwarzPrimitiveEstimator`` qualifies, giving the usable form + +.. math:: + :label: err-bq-blocked + + B_q = \left|c_ic_jc_kc_l\right| + \sqrt{\left\|M_{ij}\right\|_F}\sqrt{\left\|M_{kl}\right\|_F}. + +Note the square roots, and note that the coefficient factors appear once. +:ref:`scr-cs-implementation` describes what that module returns in its raw +form; the two differ, and the difference matters. + +**************************************** +Tightening with Sign Information +**************************************** + +The sign of a discarded contribution factorizes, + +.. math:: + :label: err-sign-factor + + \mathrm{sign}\left(t_q\right) = + \mathrm{sign}\left(c_ic_jc_kc_l\right)\cdot + \mathrm{sign}\left(\left[ij|kl\right]\right), + +and the first factor is free. ``PrimitiveNormalization`` emits one +renormalized coefficient per primitive, shared by all AO components of that +primitive, so the coefficient product is available wherever the contraction +itself is. It is also genuinely informative: standard contracted basis sets +carry coefficients of both signs, so this factor is not a constant. + +The second factor requires an argument. If the primitive charge distributions +:math:`\rho_{ij} = g_ig_j` and :math:`\rho_{kl} = g_kg_l` are each of one sign +throughout space, then :math:`\left[ij|kl\right]` is the Coulomb interaction +of two same-signed distributions and is therefore positive. This gives a +*sufficient* -- not necessary -- criterion, so the oracle must be permitted to +answer "unknown." + +Define :math:`\sigma\left(q\right) \in \left\{+1,-1,0\right\}`, with +:math:`0` denoting an indeterminate sign: + +Tier 0: all-s quartets +====================== + +If all four primitives have :math:`l = 0`, both distributions are strictly +positive -- Eq. :eq:`scr-ssss` shows this directly -- and + +.. math:: + :label: err-tier0 + + \sigma\left(q\right) = \mathrm{sign}\left(c_ic_jc_kc_l\right). + +This tier is available in any basis and requires no information beyond the +angular momenta and the coefficients. + +Tier 1: one-signed Cartesian products +===================================== + +For Cartesian primitives, :math:`\rho_{ij}` is a product of three independent +factors, one per axis. Along :math:`x` the factor is +:math:`\left(x-A_x\right)^{a_x}\left(x-B_x\right)^{b_x}`, which changes sign +at :math:`A_x` when :math:`a_x` is odd and at :math:`B_x` when :math:`b_x` is +odd. It is therefore one-signed if and only if + +* :math:`a_x` and :math:`b_x` are both even, **or** +* :math:`A_x = B_x` and :math:`a_x + b_x` is even, + +and likewise for :math:`y` and :math:`z`. When this holds for both the bra and +ket pairs, Eq. :eq:`err-tier0` again applies. + +The second clause carries most of the coverage in practice. A same-center +:math:`p_x \cdot p_x` pair gives :math:`\left(x-A_x\right)^2`, which is +one-signed, and the pairs that matter for screening are heavily concentrated +on or near a single center. The first clause alone would admit only functions +with all-even Cartesian exponents. + +This tier is unavailable for ``pure`` shells, which is a per-shell property in +this package. Real solid harmonics with :math:`l > 0` have angular nodes, so +their products are not one-signed and the criterion cannot fire. + +.. _err-cartesian-caveat: + +Cartesian does not guarantee a sign +=================================== + +It is worth being explicit that working in a Cartesian basis does *not* make +every sign determinable. Consider :math:`\left(p_xp_x|ss\right)` with the two +:math:`p` functions on distinct centers. The bra factor is +:math:`\left(x-A_x\right)\left(x-B_x\right)`, which is negative between the +centers and positive outside, so Tier 1 declines it and +:math:`\sigma\left(q\right) = 0`. + +Because one-signedness is sufficient rather than necessary, such a quartet may +still have a definite sign; the oracle simply cannot certify it cheaply. This +is a limit on coverage, not a failure of the bound. What a Cartesian basis +buys is far *more* coverage than a solid-harmonic one, where only +:math:`l = 0` qualifies at all. + +********************** +Assembling the Bracket +********************** + +Partition :math:`S\left(t\right)` by the oracle into :math:`S_+`, +:math:`S_-`, and :math:`S_0`. A determined sign collapses +:math:`\left[-B_q, +B_q\right]` to a half-interval anchored at zero: +:math:`t_q \in \left[0, B_q\right]` on :math:`S_+` and +:math:`t_q \in \left[-B_q, 0\right]` on :math:`S_-`. Substituting into Eq. +:eq:`err-generic-bracket`, + +.. math:: + :label: err-bracket + + \delta_{\max} = \sum_{q\in S_-}B_q + \sum_{q\in S_0}B_q, + \qquad + \delta_{\min} = -\sum_{q\in S_+}B_q - \sum_{q\in S_0}B_q. + +Four properties follow immediately and are worth stating together. + +**It cannot be violated.** The magnitude bound is rigorous at every angular +momentum and the oracle declines to guess. Containment is guaranteed by +construction, not by empirical calibration. + +**It is valid with no signs at all.** If :math:`S_0 = S\left(t\right)`, Eq. +:eq:`err-bracket` reduces to the symmetric envelope +:math:`\pm\sum_q B_q`. The sign oracle can therefore be introduced +incrementally, one tier at a time, with no risk of regression: adding sign +coverage can only shrink the interval. + +**It can become one-sided.** If every discarded quartet has a determined +positive sign, then :math:`\delta_{\max} = 0` and the bracket is +:math:`\left[-\sum_q B_q, 0\right]` -- the screened integral is *guaranteed* +to be an underestimate, with a rigorous bound on by how much. A symmetric +envelope can never produce a statement of this kind at any threshold. + +**It composes with interval arithmetic.** Equation :eq:`err-bracket` is an +asymmetric interval, so its midpoint is a natural best estimate and its +half-width a natural uncertainty, in the form the package's uncertainty +machinery already consumes. + +***************************** +How Much the Sign Buys +***************************** + +It is easy to overstate the value of the sign oracle, so it is worth +quantifying. The width of the bracket is + +.. math:: + :label: err-width + + W = \delta_{\max} - \delta_{\min} + = \sum_{q\in S_0}2B_q + \sum_{q\in S_+\cup S_-}B_q + = \left(1+f_0\right)\sum_{q\in S\left(t\right)}B_q, + +where + +.. math:: + :label: err-f0 + + f_0 = \frac{\sum_{q\in S_0}B_q}{\sum_{q\in S\left(t\right)}B_q} + +is the :math:`B`-weighted fraction of indeterminate quartets. Relative to the +symmetric envelope, which is the :math:`f_0 = 1` case, sign determination +narrows the bracket by a factor of :math:`2/\left(1+f_0\right)`. **The best +possible improvement in width is therefore a factor of two**, attained only +when every sign is known. + +A factor of two is real but bounded, and it is not the main reason to do this. +The qualitative payoff is asymmetry. A symmetric interval always contains +zero, so it can only ever support the statement "the error is no larger than +this." An asymmetric one may exclude zero, supporting the strictly stronger +statement "the computed value is too small, by no more than this." Direction +is information that no amount of tightening a symmetric envelope will produce. + +Two further consequences of Eq. :eq:`err-width` are worth noting. Tightening +:math:`B_q` itself -- for instance via the distance-including bounds discussed +below -- is unbounded in its potential benefit, and so is ultimately the more +valuable direction. And because :math:`f_0` rises with angular momentum, the +sign oracle contributes most in exactly the s- and p-dominated regimes where +the integrals are cheapest to begin with. + +*********************************************** +Cartesian Screening and the Pure-Basis Bracket +*********************************************** + +Tier 1 requires Cartesian primitives, but results are usually wanted in a +solid-harmonic basis. These are reconciled without needing the sign of the +transformed integral at all, because the transformation is linear and +*intervals propagate through linear maps exactly*. + +Let :math:`w_c` denote the product of the four transformation coefficients +relating a solid-harmonic quartet to the Cartesian quartet :math:`c`, and let +:math:`\left[\ell_c, u_c\right]` be the Cartesian bracket. Then + +.. math:: + :label: err-transform + + \ell_{\rm pure} = \sum_c \begin{cases} + w_c\ell_c & w_c > 0\\ w_cu_c & w_c < 0 + \end{cases}, + \qquad + u_{\rm pure} = \sum_c \begin{cases} + w_cu_c & w_c > 0\\ w_c\ell_c & w_c < 0 + \end{cases}. + +Each endpoint is obtained by pairing every coefficient with whichever end of +its interval extremizes the sum. The result is rigorous, and one-sidedness +survives wherever cancellation among the :math:`w_c` does not destroy it. + +The widening this introduces is worth naming. Equation :eq:`err-transform` +treats the Cartesian contributions as independent when they are in fact +correlated, so the resulting interval is wider than the true attainable range. +The effect is mild for :math:`d` shells, where only a few Cartesian components +contribute to each solid-harmonic function, and grows with :math:`l`. + +The transformation coefficients themselves need not be derived: libint exposes +them in sparse row-compressed form, together with a closed-form generator. The +only new machinery required is an application of those coefficients that +carries a pair of endpoints rather than a single value, since the routines +libint provides transform values. + +*************** +Diagnostics +*************** + +Three cheap quantities characterize how well the model is doing on a given +system, and are more informative than the interval width alone: + +* :math:`f_0` from Eq. :eq:`err-f0`, the :math:`B`-weighted fraction of + indeterminate quartets. This measures sign *coverage* and, by Eq. + :eq:`err-width`, fixes the width improvement exactly. +* The asymmetry ratio + :math:`\left|\delta_{\max}+\delta_{\min}\right|/ + \left(\delta_{\max}-\delta_{\min}\right)`, which is zero for a symmetric + envelope and one for a fully one-sided bracket. +* Whether the bracket excludes zero, per element. This is the binary form of + the previous item and the one that matters for reporting: it is the + condition under which a *direction* can be asserted. + +.. _err-tightness: + +***************************** +Practical Tightness Caveat +***************************** + +The dominant source of looseness in Eq. :eq:`err-bracket` is not the sign +oracle but the magnitude bound. Two effects compound. + +First, the shell-block Frobenius norm of :ref:`scr-cs-implementation` +overestimates the largest element of the block by up to the square root of the +block size. Replacing it with a per-element maximum over the block would +recover this without sacrificing rigor. + +Second, and more importantly, Cauchy-Schwarz is distance-blind. The right-hand +side of Eq. :eq:`err-bq` does not depend on +:math:`\left|\vec{P}-\vec{Q}\right|`, whereas the true integral decays as +:math:`1/\left|\vec{P}-\vec{Q}\right|` by Eq. :eq:`scr-ssss`. Since screening +discards well-separated quartets preferentially, the bound is loosest over +precisely the set it is summed across. + +*********************** +Refinements Deferred +*********************** + +Two directions were considered and are recorded here so the choice is +traceable. + +**Distance-including bounds.** The :math:`1/\left|\vec{P}-\vec{Q}\right|` +decay can be restored while preserving rigor, which addresses the second +effect in :ref:`err-tightness` directly and, per Eq. :eq:`err-width`, offers +unbounded improvement rather than the sign oracle's factor of two. Rigorous +and non-rigorous variants are both well documented in the literature +:cite:`gill1994,maurer2012,thompson2017`. + +**A Boys-damped signed estimate.** Restoring :math:`F_0\left(T\right)` to the +fine metric and keeping the coefficient signs yields, by Eq. +:eq:`scr-fine-identity`, an estimate that is *exact* for s-type quartets and +cheaper than Cauchy-Schwarz, requiring no hyper-diagonal integrals. It was +rejected as the starting point for two reasons: it is a prediction rather than +a bound, so it supports no guarantee; and it inherits the fine metric's +failure to bound :math:`l > 0` contributions, which is the very behavior that +motivated moving off the fine metric. It remains attractive as a *bias +estimate* to accompany the bracket, and could be added without disturbing the +bound. + +******************************************* +Relation to the Existing Modules +******************************************* + +``PrimitiveErrorModel`` implements the unsigned predecessors of this model. It +walks the same decontracted quartet loop as ``PrimitiveContractor`` and +applies the same gates, accumulating one of three per-quartet estimates for +each skipped quartet. In the language of this page, each is a choice of +:math:`B_q` with :math:`S_0 = S\left(t\right)` -- that is, a symmetric +envelope with no sign information: + +* ``"Tolerance"`` takes :math:`B_q = t`, which is not a bound at all but a + nominal per-quartet charge; +* ``"Coarse"`` takes :math:`B_q = \bar{K}_{ij}\bar{K}_{kl}`; +* ``"Fine"`` takes + :math:`B_q = \left|Q_{ij}Q_{kl}\right|/\sqrt{\gamma_{ij}+\gamma_{kl}}`, + subject to the :math:`l > 0` caveat of :ref:`scr-rigor`. + +Validation +========== + +``AnalyticError`` computes the true signed error directly, as the difference +between a screened result and a tight-threshold benchmark. It therefore +supplies exactly the quantity Eq. :eq:`err-delta` defines, element by element, +and the model can be checked against it without any statistical machinery. + +The acceptance criterion is containment, and it is absolute: every element of +the true error must lie within :math:`\left[\delta_{\min},\delta_{\max}\right]`, +with **zero** violations across the test set. Only once containment holds is +tightness meaningful, at which point the diagnostics above quantify it. A +model that is tight but occasionally violated is strictly less useful than one +that is loose but never violated, since only the latter supports a guarantee. +The ``primitive_error_models`` example provides a driver for these +comparisons. + +*************** +Limitations +*************** + +Sign coverage degrades as angular momentum rises. Tier 0 requires an all-s +quartet; Tier 1 requires one-signed Cartesian products and is unavailable for +``pure`` shells. For basis sets dominated by high-:math:`l` solid-harmonic +functions, :math:`f_0 \rightarrow 1` and Eq. :eq:`err-bracket` degenerates to +the symmetric envelope. The model does not fail in that limit -- it simply +stops adding value over what a magnitude bound alone provides. + +The bound also says nothing about how the error propagates into derived +quantities. Screening errors in |eri4| enter the Fock matrix and the total +energy through contractions with the density, and bounding those requires +tracking correlations between elements that Eq. :eq:`err-bracket` treats +independently. That analysis is outside the scope of this page. + +********** +References +********** + +.. bibliography:: + :style: unsrt + :filter: docname in docnames + :labelprefix: ERR diff --git a/docs/source/background/eris/index.rst b/docs/source/background/eris/index.rst new file mode 100644 index 00000000..8a768a49 --- /dev/null +++ b/docs/source/background/eris/index.rst @@ -0,0 +1,24 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +############################ +Electron Repulsion Integrals +############################ + +.. toctree:: + :maxdepth: 2 + :caption: Contents + + eris + density_fitting diff --git a/docs/source/background/eris/screening.rst b/docs/source/background/eris/screening.rst new file mode 100644 index 00000000..585b7a89 --- /dev/null +++ b/docs/source/background/eris/screening.rst @@ -0,0 +1,316 @@ +.. Copyright 2026 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +########################## +Primitive Screening Models +########################## + +.. |eri4| replace:: :math:`\left(\mu\nu|\lambda\sigma\right)` +.. |Kij| replace:: :math:`K_{ij}` +.. |Qij| replace:: :math:`Q_{ij}` + +Most of the |eri4| defined in :doc:`eris` are negligible. Screening is the +practice of predicting, before an integral is computed, that its contribution +falls below a threshold and may be skipped. This page collects the screening +models used by the ``Integrals`` package, states each one precisely, and -- +in :ref:`scr-rigor` -- records which of them are actually rigorous bounds. +That last question turns out to matter a great deal, because a model that is +adequate for *deciding what to skip* is not automatically adequate for +*quantifying what skipping cost*. :doc:`error_bounds` takes up the latter +problem. + +************************************** +Notation and the Primitive Expansion +************************************** + +Each contracted AO :math:`\phi_\mu` is a fixed linear combination of primitive +Gaussians. Writing the expansion out inside Eq. :eq:`eri4` turns a contracted +ERI into a sum over quartets of primitives: + +.. math:: + :label: scr-prim-expansion + + \left(\mu\nu|\lambda\sigma\right) = + \sum_{i\in\mu}\sum_{j\in\nu}\sum_{k\in\lambda}\sum_{l\in\sigma} + c_i c_j c_k c_l \left[ij|kl\right], + +where square brackets denote an integral over primitives and :math:`c_i` is +the contraction coefficient of primitive :math:`i`. Throughout, the +:math:`c_i` are the *renormalized* coefficients returned by the +``Normalize`` property type -- the same values ``PrimitiveContractor`` +consumes through its ``Primitive Normalization`` submodule. These absorb both +the per-primitive normalization factor and the unit-norm scaling of the +contracted shell, so Eq. :eq:`scr-prim-expansion` holds as written, with no +further normalization factors. + +Screening operates on Eq. :eq:`scr-prim-expansion` one quartet at a time, so +the quantities of interest are properties of *pairs*. For two primitives with +exponents :math:`\alpha_i`, :math:`\alpha_j` centered at :math:`\vec{A}`, +:math:`\vec{B}`, the Gaussian product theorem gives + +.. math:: + :label: scr-gpt + + \gamma_{ij} = \alpha_i + \alpha_j, \qquad + \rho_{ij} = \frac{\alpha_i\alpha_j}{\gamma_{ij}}, \qquad + \vec{P}_{ij} = \frac{\alpha_i\vec{A} + \alpha_j\vec{B}}{\gamma_{ij}}, + +and the product of the two radial parts is a single Gaussian of exponent +:math:`\gamma_{ij}` centered at :math:`\vec{P}_{ij}`, scaled by the overlap +factor + +.. math:: + :label: scr-kij + + K_{ij} = \exp\left(-\rho_{ij}\left|\vec{A}-\vec{B}\right|^2\right). + +|Kij| decays as the square of the distance between the two centers, and it is +the origin of essentially all screening: a pair of primitives on distant +centers has an exponentially small product density, so any integral involving +that pair is exponentially small. + +******************************* +The Primitive [ss|ss] Integral +******************************* + +For four s-type primitives the integral in Eq. :eq:`scr-prim-expansion` has a +closed form. With :math:`\vec{P}` and :math:`\vec{Q}` the product centers of +the bra and ket pairs, + +.. math:: + :label: scr-ssss + + \left[ij|kl\right] = + \frac{2\pi^{5/2}} + {\gamma_{ij}\gamma_{kl}\sqrt{\gamma_{ij}+\gamma_{kl}}} + K_{ij}K_{kl}\,F_0\left(T\right), \qquad + T = \frac{\gamma_{ij}\gamma_{kl}}{\gamma_{ij}+\gamma_{kl}} + \left|\vec{P}-\vec{Q}\right|^2, + +where :math:`F_0` is the zeroth-order Boys function :cite:`boys1950`, + +.. math:: + :label: scr-boys + + F_0\left(T\right) = \int_0^1 \exp\left(-Tu^2\right)du + = \frac{1}{2}\sqrt{\frac{\pi}{T}}\,\mathrm{erf}\left(\sqrt{T}\right). + +Two properties of :math:`F_0` are used repeatedly below. It is bounded, +:math:`0 < F_0\left(T\right) \leq F_0\left(0\right) = 1`, and for large +:math:`T` it decays as +:math:`F_0\left(T\right) \rightarrow \frac{1}{2}\sqrt{\pi/T}`, which is the +classical :math:`1/\left|\vec{P}-\vec{Q}\right|` Coulomb tail. Equation +:eq:`scr-ssss` also makes plain that the s-type primitive ERI is *strictly +positive*, a fact :doc:`error_bounds` relies on. + +********************** +Coarse Pair Estimates +********************** + +The cheapest useful estimate discards everything in Eq. :eq:`scr-ssss` except +the exponential factors. ``BlackBoxPrimitiveEstimator`` computes + +.. math:: + :label: scr-coarse + + \bar{K}_{ij} = \left|c_i\right|\left|c_j\right|K_{ij}, + +which is the quantity libint uses for coarse screening and which the +``Integrals`` package builds with the ``coarse_k_ij`` helper. It requires only +exponents, centers, and coefficients -- no integrals at all -- so it can be +tabulated once per basis-set pair and reused for every quartet. + +****************** +The Fine Metric +****************** + +The next refinement restores the exponent-dependent prefactor. libint's +``ShellPair`` data carries, and the ``fine_k_ij`` helper reproduces, + +.. math:: + :label: scr-fine-pair + + Q_{ij} = \bar{K}_{ij}\,\frac{\sqrt{2}\,\pi^{5/4}}{\gamma_{ij}}, + +from which the quartet-level fine screening quantity is formed as +:math:`\left|Q_{ij}Q_{kl}\right|/\sqrt{\gamma_{ij}+\gamma_{kl}}`. Multiplying +Eq. :eq:`scr-fine-pair` out for the bra and ket pairs gives + +.. math:: + :label: scr-fine-identity + + \frac{Q_{ij}Q_{kl}}{\sqrt{\gamma_{ij}+\gamma_{kl}}} = + \left|c_ic_jc_kc_l\right| + \frac{2\pi^{5/2}} + {\gamma_{ij}\gamma_{kl}\sqrt{\gamma_{ij}+\gamma_{kl}}} + K_{ij}K_{kl}, + +since :math:`\sqrt{2}\cdot\sqrt{2} = 2` and +:math:`\pi^{5/4}\cdot\pi^{5/4} = \pi^{5/2}`. Comparing with Eq. +:eq:`scr-ssss`, the fine metric is *exactly* the contracted contribution of an +s-type primitive quartet with :math:`F_0\left(T\right)` replaced by one. It is +the :math:`\left|\vec{P}-\vec{Q}\right|\rightarrow 0` limit of the true +integral. + +This identity explains both the strength and the weakness of the fine metric. +Because :math:`F_0 \leq 1`, it is an upper bound on the s-type contribution, +and a cheap one. But it is a *distance-blind* upper bound: it discards +:math:`F_0\left(T\right)` entirely, and :math:`F_0` is small precisely for the +well-separated quartets that screening is most interested in. The fine metric +is therefore loosest exactly where it is used most. + +******************** +Cauchy-Schwarz +******************** + +A qualitatively different estimate follows from the observation that the +electron-repulsion operator defines a positive-definite bilinear form on +charge distributions. The Cauchy-Schwarz inequality then applies directly +:cite:`haser1989`: + +.. math:: + :label: scr-cauchy-schwarz + + \left|\left[ij|kl\right]\right| \leq + \sqrt{\left[ij|ij\right]}\sqrt{\left[kl|kl\right]}. + +Unlike Eqs. :eq:`scr-coarse` and :eq:`scr-fine-pair`, this holds for *any* +angular momentum, Cartesian or solid-harmonic, with no restriction whatsoever +-- it is a property of the operator, not of the basis functions. The cost is +that it requires the hyper-diagonal integrals :math:`\left[ij|ij\right]`, +which must actually be computed, though only :math:`O(N^2)` of them. + +``CauchySchwarzPrimitiveEstimator`` computes this family of quantities. Its +implementation is described in :ref:`scr-cs-implementation` below, since the +precise quantity it returns matters when the estimate is reused as an error +bound. + +Cauchy-Schwarz shares the fine metric's distance blindness: the right-hand +side of Eq. :eq:`scr-cauchy-schwarz` does not depend on +:math:`\left|\vec{P}-\vec{Q}\right|` at all, so it cannot capture the +:math:`1/\left|\vec{P}-\vec{Q}\right|` decay of a well-separated quartet. +Recovering that decay while preserving rigor is the subject of the +distance-including estimates, which the package does not currently implement; +:doc:`error_bounds` surveys them and explains why they are the most valuable +direction for future work. + +.. _scr-cs-implementation: + +Implementation Notes +==================== + +``CauchySchwarzPrimitiveEstimator`` decontracts both basis sets, computes the +primitive quartet tensor for the hyper-diagonal, and reduces each +:math:`\left(ij|ij\right)` shell block to a single number with +``utils::rank2_shell_norm``. Two details are worth stating explicitly, because +both affect how the result may be used: + +#. The reduction is a **Frobenius norm over the shell block**, not a maximum + over its elements. Writing :math:`M_{ij}` for the matrix of + :math:`\left(ij|ij\right)` values within the block, the module forms + :math:`\left\|M_{ij}\right\|_F`. Since :math:`M_{ij}` is a Gram matrix, its + Frobenius norm is at least its largest diagonal element, so the result is a + valid bound on every element of the block, but a loose one -- loose by up to + the square root of the block size. + +#. The returned value is scaled by the contraction coefficients, + :math:`\left|c_ic_j\right|\left\|M_{ij}\right\|_F`, and is **not** + square-rooted. This is the natural form for the ``PairScreener`` use case, + where a per-pair magnitude is compared against a tolerance directly. It is + *not* the form that appears in Eq. :eq:`scr-cauchy-schwarz`: assembling a + Cauchy-Schwarz product from it requires taking the square root of the norm + and being careful not to apply the coefficient factors twice. + :doc:`error_bounds` writes the bound in terms of the unscaled norm to keep + this unambiguous. + +There is no double counting of normalization between the two factors. The +decontracted basis carries unit contraction coefficients but is converted to +libint with normalization embedded, so :math:`M_{ij}` is built from normalized +primitives, while :math:`c_i` supplies the contraction coefficient exactly +once. + +.. _scr-rigor: + +************************** +Rigor of These Estimates +************************** + +The three models above are not interchangeable, and the difference is not +merely one of tightness. + +Equations :eq:`scr-coarse` and :eq:`scr-fine-pair` are derived from the +s-type closed form, Eq. :eq:`scr-ssss`. Neither retains any part of the +angular prefactor that a primitive with :math:`l > 0` contributes. For an +s-type quartet the fine metric is a rigorous upper bound, by the identity +:eq:`scr-fine-identity` together with :math:`F_0 \leq 1`. For :math:`l > 0` +no such guarantee exists, and the true integral may exceed the estimate. This +is not a defect of the implementation: libint's ``ScreeningMethod::Original``, +which ``PrimitiveContractor`` reproduces, is a screening heuristic chosen for +its cost, and it is documented as such. Used at a tight threshold it discards +only genuinely negligible terms. Used as the basis of an *error model*, +however, it produces estimates that the true error can and does exceed. + +Equation :eq:`scr-cauchy-schwarz` has no such restriction. It is a rigorous +bound at every angular momentum. This is the reason :doc:`error_bounds` builds +on Cauchy-Schwarz rather than on the metric that the screening gates +themselves use. + +****************************** +The Gates as Implemented +****************************** + +``PrimitiveContractor`` evaluates Eq. :eq:`scr-prim-expansion` directly, and +applies four tests in sequence. Writing :math:`t` for the threshold, a quartet +is skipped if any of the following holds: + +#. :math:`\bar{K}_{ij} < t`, discarding the bra pair outright; +#. :math:`\bar{K}_{kl} < t`, discarding the ket pair outright; +#. :math:`\bar{K}_{ij}\bar{K}_{kl} \leq t`; +#. :math:`\left|Q_{ij}Q_{kl}\right|/\sqrt{\gamma_{ij}+\gamma_{kl}} < t`. + +The first three are the coarse tests of Eq. :eq:`scr-coarse` and the last is +the fine test of Eq. :eq:`scr-fine-pair`. The threshold is exposed as the +``Screening Threshold`` input and defaults to :math:`10^{-16}`; the sequence +matches libint's ``ScreeningMethod::Original``. + +A separate module, ``ScreenPrimitivePairs``, applies a single tolerance test +to a matrix of pair estimates supplied by any ``PrimitivePairEstimator`` +submodule, and returns the surviving pair indices. This is the composable +entry point: it is agnostic to which of the estimates above is used. + +**************************** +From Screening to Error +**************************** + +Screening and error estimation ask different questions of the same +quantities. Screening asks whether a contribution is small enough to skip, and +a good heuristic answer suffices -- if the estimate is occasionally too small, +the consequence is a slightly larger error than the threshold nominally +promises. Error estimation asks what the accumulated cost of every skipped +contribution actually was, and here a heuristic will not do: an estimate that +the true value can exceed provides no guarantee at all. + +Nothing requires the two to use the same model. The screening gates may +remain the cheap, distance-blind heuristics documented above while the error +model is built on the rigorous inequality of Eq. :eq:`scr-cauchy-schwarz`. +:doc:`error_bounds` develops that error model. + +********** +References +********** + +.. bibliography:: + :style: unsrt + :filter: docname in docnames + :labelprefix: SCR diff --git a/docs/source/background/index.rst b/docs/source/background/index.rst new file mode 100644 index 00000000..031a7370 --- /dev/null +++ b/docs/source/background/index.rst @@ -0,0 +1,26 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +########## +Background +########## + +.. toctree:: + :maxdepth: 2 + :caption: Contents + + overlap + kinetic + nuclear_attraction + eris/index diff --git a/docs/source/background/kinetic.rst b/docs/source/background/kinetic.rst new file mode 100644 index 00000000..8a93bad9 --- /dev/null +++ b/docs/source/background/kinetic.rst @@ -0,0 +1,28 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +######################## +Kinetic Energy Integrals +######################## + +Over a set of :math:`N_b` atomic orbitals (AOs), +:math:`\lbrace\phi_\mu\left(\vec{r}\right)\rbrace`, the electronic kinetic +energy integrals are: + +.. math:: + + T_{\mu\nu} = \int + \phi_\mu\left(\vec{r}\right) + \left(-\frac{1}{2}\nabla^2\right) + \phi_\nu\left(\vec{r}\right) d\vec{r}. diff --git a/docs/source/background/nuclear_attraction.rst b/docs/source/background/nuclear_attraction.rst new file mode 100644 index 00000000..c34d28b4 --- /dev/null +++ b/docs/source/background/nuclear_attraction.rst @@ -0,0 +1,29 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +##################################### +Electron-Nuclear Attraction Integrals +##################################### + +Over a set of :math:`N_b` atomic orbitals (AOs), +:math:`\lbrace\phi_\mu\left(\vec{r}\right)\rbrace`, and a set of nuclei with +charges :math:`Z_A` located at :math:`\vec{R}_A`, the electron-nuclear +attraction integrals are: + +.. math:: + + V^{en}_{\mu\nu} = -\sum_A Z_A \int + \phi_\mu\left(\vec{r}\right) + \frac{1}{\left|\vec{r}-\vec{R}_A\right|} + \phi_\nu\left(\vec{r}\right) d\vec{r}. diff --git a/docs/source/background/overlap.rst b/docs/source/background/overlap.rst new file mode 100644 index 00000000..4617378f --- /dev/null +++ b/docs/source/background/overlap.rst @@ -0,0 +1,26 @@ +.. Copyright 2025 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +################# +Overlap Integrals +################# + +Over a set of :math:`N_b` atomic orbitals (AOs), +:math:`\lbrace\phi_\mu\left(\vec{r}\right)\rbrace`, the overlap integrals +are: + +.. math:: + + S_{\mu\nu} = \int + \phi_\mu\left(\vec{r}\right)\phi_\nu\left(\vec{r}\right) d\vec{r}. diff --git a/docs/source/conf.py b/docs/source/conf.py index 28005c82..f5ada731 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -74,8 +74,13 @@ "sphinx.ext.autosummary", "sphinx_rtd_theme", "sphinx.ext.intersphinx", + "sphinxcontrib.bibtex", ] +# -- Options for sphinxcontrib-bibtex ----------------------------------------- + +bibtex_bibfiles = ["references.bib"] + # Add any paths that contain templates here, relative to this directory. # templates_path = ['_templates'] diff --git a/docs/source/index.rst b/docs/source/index.rst index 417c1238..f61dfbbf 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -19,5 +19,6 @@ Welcome to Integral's documentation! :maxdepth: 2 :caption: Contents: + background/index module_api/index C++ API diff --git a/docs/source/references.bib b/docs/source/references.bib new file mode 100644 index 00000000..d2dd585b --- /dev/null +++ b/docs/source/references.bib @@ -0,0 +1,68 @@ +@article{whitten1973, + author = {Whitten, J. L.}, + title = {Coulombic potential energy integrals and approximations}, + journal = {The Journal of Chemical Physics}, + year = {1973}, + volume = {58}, + number = {10}, + pages = {4496--4501}, + doi = {10.1063/1.1679012} +} + +@article{dunlap1979, + author = {Dunlap, B. I. and Connolly, J. W. D. and Sabin, J. R.}, + title = {On some approximations in applications of {X}$\alpha$ theory}, + journal = {The Journal of Chemical Physics}, + year = {1979}, + volume = {71}, + number = {8}, + pages = {3396--3402}, + doi = {10.1063/1.438728} +} + +@article{vahtras1993, + author = {Vahtras, Olav and Alml{\"o}f, Jan and Feyereisen, Michael W.}, + title = {Integral approximations for {LCAO-SCF} calculations}, + journal = {Chemical Physics Letters}, + year = {1993}, + volume = {213}, + number = {5-6}, + pages = {514--518}, + doi = {10.1016/0009-2614(93)89151-7} +} + +@article{feyereisen1993, + author = {Feyereisen, Michael and Fitzgerald, George and Komornicki, Andrew}, + title = {Use of approximate integrals in ab initio theory. {A}n + application in {MP2} energy calculations}, + journal = {Chemical Physics Letters}, + year = {1993}, + volume = {208}, + number = {5-6}, + pages = {359--363}, + doi = {10.1016/0009-2614(93)87156-W} +} + +@article{weigend2002, + author = {Weigend, Florian}, + title = {A fully direct {RI-HF} algorithm: {I}mplementation, optimised + auxiliary basis sets, demonstration of accuracy and efficiency}, + journal = {Physical Chemistry Chemical Physics}, + year = {2002}, + volume = {4}, + number = {18}, + pages = {4285--4291}, + doi = {10.1039/b204199p} +} + +@article{weigend2009, + author = {Weigend, Florian and Kattannek, Marco and Ahlrichs, Reinhart}, + title = {Approximated electron repulsion integrals: {C}holesky + decomposition versus resolution of the identity methods}, + journal = {The Journal of Chemical Physics}, + year = {2009}, + volume = {130}, + number = {16}, + pages = {164106}, + doi = {10.1063/1.3116103} +} From 2b8c3a5e3ce346410c2955ba0f22a455a020524b Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 17 Aug 2026 09:24:21 -0500 Subject: [PATCH 2/2] Add ERI primitive screening and error-bound documentation Adds two background pages under docs/source/background/eris/: - screening.rst documents the primitive screening models the package already implements (coarse K_ij, libint's fine metric, Cauchy-Schwarz) and records which are actually rigorous bounds. The coarse and fine metrics derive from the s-type closed form and retain no angular prefactor, so they do not bound l > 0 contributions; Cauchy-Schwarz does, at every angular momentum. - error_bounds.rst derives a rigorous two-sided bracket on the error from neglecting primitive contributions. A Cauchy-Schwarz magnitude bound brackets each discarded quartet symmetrically; a three-valued sign oracle collapses that to a half-interval where the sign can be established a priori. The bound is valid with no sign information and tightens monotonically as coverage grows. Scopes the existing bibliography directives per-document and gives each a label prefix so the three pages cite disjoint key sets without duplicate-label warnings under -W. Co-Authored-By: Claude Opus 5 --- .../background/eris/density_fitting.rst | 2 + docs/source/background/eris/index.rst | 2 + docs/source/references.bib | 60 +++++++++++++++++++ 3 files changed, 64 insertions(+) diff --git a/docs/source/background/eris/density_fitting.rst b/docs/source/background/eris/density_fitting.rst index d739b8ed..e912c7b3 100644 --- a/docs/source/background/eris/density_fitting.rst +++ b/docs/source/background/eris/density_fitting.rst @@ -217,3 +217,5 @@ References .. bibliography:: :style: unsrt + :filter: docname in docnames + :labelprefix: DF diff --git a/docs/source/background/eris/index.rst b/docs/source/background/eris/index.rst index 8a768a49..dabcb136 100644 --- a/docs/source/background/eris/index.rst +++ b/docs/source/background/eris/index.rst @@ -22,3 +22,5 @@ Electron Repulsion Integrals eris density_fitting + screening + error_bounds diff --git a/docs/source/references.bib b/docs/source/references.bib index d2dd585b..2f03614a 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -66,3 +66,63 @@ @article{weigend2009 pages = {164106}, doi = {10.1063/1.3116103} } + +@article{boys1950, + author = {Boys, S. F.}, + title = {Electronic wave functions. {I}. {A} general method of + calculation for the stationary states of any molecular system}, + journal = {Proceedings of the Royal Society of London. Series A}, + year = {1950}, + volume = {200}, + number = {1063}, + pages = {542--554}, + doi = {10.1098/rspa.1950.0036} +} + +@article{haser1989, + author = {H{\"a}ser, Marco and Ahlrichs, Reinhart}, + title = {Improvements on the direct {SCF} method}, + journal = {Journal of Computational Chemistry}, + year = {1989}, + volume = {10}, + number = {1}, + pages = {104--111}, + doi = {10.1002/jcc.540100111} +} + +@article{gill1994, + author = {Gill, Peter M. W. and Johnson, Benny G. and Pople, John A.}, + title = {A simple yet powerful upper bound for {C}oulomb integrals}, + journal = {Chemical Physics Letters}, + year = {1994}, + volume = {217}, + number = {1-2}, + pages = {65--68}, + doi = {10.1016/0009-2614(93)E1340-M} +} + +@article{maurer2012, + author = {Maurer, Simon A. and Lambrecht, Daniel S. and Flaig, Denis and + Ochsenfeld, Christian}, + title = {Distance-dependent {S}chwarz-based integral estimates for + two-electron integrals: {R}eliable tightness vs. rigorous upper + bounds}, + journal = {The Journal of Chemical Physics}, + year = {2012}, + volume = {136}, + number = {14}, + pages = {144107}, + doi = {10.1063/1.3693908} +} + +@article{thompson2017, + author = {Thompson, Travis H. and Ochsenfeld, Christian}, + title = {Distance-including rigorous upper bounds and tight estimates for + two-electron integrals over long- and short-range operators}, + journal = {The Journal of Chemical Physics}, + year = {2017}, + volume = {147}, + number = {14}, + pages = {144101}, + doi = {10.1063/1.4994190} +}