Finite primitive and contracted Gaussian integral evaluation with certified screening

Inputs. A finite list of input basis contractions ${\varphi }_a={\sum }_{i=1}^{q_a}{d}_{ai}{g}_{ai}$, $1\le a\le n$, from Cartesian Gaussian primitives, contractions, and integral conventions, represented by finite dyadic intervals for real centers and positive exponents and finite dyadic intervals or complex rectangles for coefficients; a finite nuclear list $((C_{\mathrm{\ell }},w_{\mathrm{\ell }}){)}_{\mathrm{\ell }=1}^J$ whose centers and signed real weights are represented by finite dyadic intervals; a requested subset $\mathcal{I}\subseteq \{S,T,V,G\}$ of integral outputs (overlap, kinetic, nuclear attraction, and Coulomb); and Boolean energy and gradient requests, with a gradient request requiring an energy request. Require at least one requested output. Supply a positive integral tolerance ${\epsilon }_I$ for every requested integral and for integrals used internally by an energy request; supply a positive energy tolerance ${\epsilon }_E$ when energy is requested and, for a gradient request, an optional positive gradient tolerance defaulting to ${\epsilon }_E$. An optional positive screening tolerance is allowed only for an energy request. Also supply finite positive integer work and precision caps $W,P$ and a dyadic format of precision $2\le p\le P$ with exponent bounds and fixed elementary truncation index $K\ge 1$. A primitive request is the singleton-contraction case. An energy/gradient request additionally requires $n,r\ge 1$, real primitive-contraction coefficients $d_{ai}$ and real nuclear data, real matrices $C\in {\mathbb{R}}^{n\times r},D\in {\mathbb{R}}^{r\times r}$ given by dyadic boxes (exact matrices are singleton boxes), and uniform finite bounds $B\ge \sup \mathrm{\parallel }C{\mathrm{\parallel }}_F$, $M_D\ge \sup \mathrm{\parallel }D{\mathrm{\parallel }}_F$ from Finite tensor contractions and their variational energy and gradient. Internal primitive labels $(a,i)$ are distinct from the input basis labels $a$; $C$ is a subsequent transform of the input basis.

Outputs. A status success, inconclusive, domain, or range; outward real intervals or complex rectangles for the integrals in $\mathcal{I}$; an energy interval when energy is requested; a gradient center and Frobenius gradient-error bound when gradient is requested; and a separately accumulated discarded-energy bound when screening is enabled. Success encloses every parameter value in the input boxes. No status declares an omitted integral zero.

Finite calls and accounting. Use Outward-rounded interval, rectangle and ball operations for arithmetic. Every real elementary call means a typed DAG call to Certified range evaluation for elementary real and complex functions at the supplied format and $K$. Enclose $\pi$ by the Machin identity in Trigonometric identities, Machin's formula and reduced Taylor tails and the $K$-term arctangent series with geometric tails from Range-reduced exponential, logarithm and arctangent series have certified tails, using exact rational tail endpoints before final outward conversion. A shared counter charges each conversion, comparison, loop update, arithmetic operation, recurrence update and summand, including every primitive operation executed by a child. Check the counter before executing an operation and interrupt any child at exhaustion. A required Boys call receives remaining work $R$ as its finite work cap, precision cap $P$, the same format and $K$, absolute tolerance ${\epsilon }_I\mathrm{/}2$ and no relative request; its actual work is charged to the parent. A required child non-success on otherwise valid data or an intermediate overflow makes the current integral evaluation inconclusive; a rounded intermediate failing a domain certification on valid original data is also inconclusive. Return range only when a certified lower bound for a requested exact output component exceeds the largest finite magnitude of the format or a certified upper bound is below its negative. A mere overflow of a conservative intermediate or enclosure endpoint is not enough to claim range. Invalid original inputs return domain; cap exhaustion returns inconclusive. These choices are conservative and need not permit success under every cap or tolerance.

1. Validate the request specification, positive exponent lower endpoints, finite input boxes, indices in ${\mathbb{N}}_0^3$, supported format and positive tolerances/caps. Invalid original data return domain. If energy/gradient is requested, verify realness and dimensions and certify the norm bounds by outward Frobenius upper bounds on the complete boxes. If this sufficient norm test cannot certify the supplied bounds, return inconclusive; do not claim that an unknown represented matrix violates its bound. Form finite primitive-pair and quartet lists for the requested integrals and any energy/gradient dependencies, lexicographically. Include all shifted overlap pairs required by the kinetic raising/lowering formula; a negative index contributes zero.

2. For each ordered primitive pair with exponents $a,b$, indices $\alpha ,\beta$ and centers $A,B$, compute $p=a+b$, $P=(aA+bB)\mathrm{/}p$, and $K_{ab}=\exp (-ab\mathrm{\parallel }A-B{\mathrm{\parallel }}^2\mathrm{/}p)$ by outward arithmetic and the real elementary call. Compute the Hermite coefficients by Finite Hermite expansion and coefficient recurrences for a Cartesian Gaussian pair, raising the first primitive in x,y,z order and then the second, storing each intermediate degree table until its successors have been formed. The final table has ${\prod }_{j=1}^3({\alpha }_j+{\beta }_j+1)$ entries; every intermediate table update is also charged. Compute overlap by the finite shifted polynomial expansion and the moments in Cartesian Gaussian moments and exponent-parameter derivatives, and kinetic by Cartesian Gaussian raising, lowering, and kinetic-energy recurrences. Certified square roots and integer powers evaluate all half-integer factors in those moment formulas.

3. To evaluate an attraction primitive pair, put $L=\mathrm{\mid }\alpha \mathrm{\mid }+\mathrm{\mid }\beta \mathrm{\mid }$ and use exponent $\eta =p$ and center difference $P-C_{\mathrm{n}\mathrm{u}\mathrm{c}}$. To evaluate a primitive Coulomb quartet, including any diagonal self quartet, put $L=\mathrm{\mid }\alpha \mathrm{\mid }+\mathrm{\mid }\beta \mathrm{\mid }+\mathrm{\mid }\gamma \mathrm{\mid }+\mathrm{\mid }\delta \mathrm{\mid }$, $\eta =pq\mathrm{/}(p+q)$ and center difference $P-Q$, where the second pair has exponent $q$. Form the nonnegative Boys argument interval $\mathcal{T}=\eta \mathrm{\parallel }P-Q{\mathrm{\parallel }}^2$, with the evident attraction substitution, using interval squaring and the positive exponent bounds. Request Boys orders $0,\dots ,L$ from Finite-budget certified evaluation of Boys function tables. If rounding prevents certification of a required positive denominator, return inconclusive.

4. Initialize $R_{000}^j=(-2\eta {)}^j{F}_j(\mathcal{T})$ for $0\le j\le L$. For total derivative degree $d=1,\dots ,L$, enumerate $(t,u,v)$ lexicographically with $t+u+v=d$, choose its first nonzero coordinate in x,y,z order, and apply that coordinate's recurrence in Hermite recurrence formulas for Gaussian nuclear-attraction and two-electron Coulomb integrals for $0\le j\le L-d$. Negative-index terms are zero. For example, when $t>0$, use $R_{tuv}^j=(P_1-Q_1)R_{t-1,u,v}^{j+1}+(t-1)R_{t-2,u,v}^{j+1}$. Each predecessor is at a smaller derivative degree and within its stored Boys-order range. Contract the finite Hermite sums in that theorem to get each primitive integral.

5. Contract primitive integrals to input-basis integrals by the finite coefficient sums. For overlap, kinetic and attraction the coefficient is $\overline{d_{ai}}{d}_{bj}$. For the defined pair-density Coulomb convention, the primitive quartet coefficient is $d_{ai}\overline{d_{bj}}\overline{d_{ck}}{d}_{dl}$: expand $\overline{f_{ab}(x)}{f}_{cd}(y)$ in Cartesian Gaussian primitives, contractions, and integral conventions. In the real energy branch these are the four real contraction coefficients. Cache only identical ordered parameter boxes, angular indices, format and truncation context; an $R$ cache record also stores its cutoff $L$. Reuse a table only if its stored cutoff is at least the needed cutoff; otherwise compute the larger table under the counter.

6. Optional Coulomb-energy screening applies to whole input-basis quartets, before their off-diagonal primitive quartets are evaluated. First compute every required contracted diagonal $G_{ab,ab}$ without screening by steps 2–5, including all coefficients $d$. Obtain finite nonnegative upper endpoints $U_{ab}$ from their real projections, justified by Positive Coulomb bilinear form, Cauchy–Schwarz, and screening bound. Form the outward real interval $W_{ab}={\sum }_{p,q}{C}_{ap}{D}_{pq}{C}_{bq}$ and its upper absolute bound $A_{ab}=\max (\mathrm{\mid }\inf W_{ab}\mathrm{\mid },\mathrm{\mid }\sup W_{ab}\mathrm{\mid })$. Compute nonnegative outward upper bounds $Q_{abcd}\ge \sqrt{U_{ab}{U}_{cd}}$ and $w_{abcd}\ge \frac{1}{2}{A}_{ab}{A}_{cd}{Q}_{abcd}$. Enumerate all ordered candidate basis quartets lexicographically, with no unrecorded symmetry compression. Skip a quartet only if the outward accumulated sum after adding $w_{abcd}$ is at most the screening tolerance; otherwise compute its primitive sums. Store a skipped real entry as $[-Q_{abcd},Q_{abcd}]$ and report the actual accumulated weighted sum. Thus entrywise radius uses $Q$, while the discarded-energy bound uses $w$. Even a skipped integral must pass its requested integral-tolerance test. No screening is required or performed when this optional real-energy branch is absent.

7. If energy or gradient is requested, form the one-body matrix on the input basis by $H_{ab}=T_{ab}+{\sum }_{\mathrm{\ell }=1}^J{w}_{\mathrm{\ell }}{V}_{C_{\mathrm{\ell }}}({\varphi }_a,{\varphi }_b)$, using the supplied signed nuclear weights, and form the Coulomb tensor $G$ on that basis. Use the midpoint tensors $\tilde{H},\tilde{G}$ and largest entrywise radii ${\epsilon }_2,{\epsilon }_4$, including every screened interval. Evaluate the finite polynomial $E(C;D,\tilde{H},\tilde{G})$ by outward arithmetic over the full $C,D$ boxes, and enlarge its interval by ${\delta }_E=M_{D}rn{B}^2{\epsilon }_2+\frac{1}{2}{M}_D^2{r}^2{n}^2{B}^4{\epsilon }_4$ from Explicit finite-dimensional contraction error bounds for energy and gradient. If gradient is requested, evaluate every formal partial derivative of this same finite polynomial by differentiating its two and four factors of $C$ and using outward finite sums. If their intervals have midpoints $c_{a\mu }$ and radii ${\rho }_{a\mu }$, return gradient center $c$ with outward upper Frobenius error bound $\sqrt{{\sum }_{a,\mu }{\rho }_{a\mu }^2}+{\delta }_{\mathrm{\nabla }}$, where ${\delta }_{\mathrm{\nabla }}=2nrB{M}_D{\epsilon }_2+2n^2{r}^2{B}^3{M}_D^2{\epsilon }_4$. This includes uncertainty in $C,D$ as well as integral error; nuclear and contraction-coefficient uncertainty was already included in the integral boxes. Count all these operations. If a final certified bound proves that a requested output range lies wholly outside the finite format, return range; a final endpoint overflow without this proof returns inconclusive. Return success only when every requested integral radius, energy radius and gradient-error bound meets its own declared tolerance; all remaining failures return inconclusive.

An $R$ table with cutoff $L$ stores ${\sum }_{d=0}^L(L-d+1)\left(\genfrac{}{}{0px}{}{d+2}{2}\right)$ entries. All input, pair, primitive-quartet and auxiliary tables are finite. Every update and child evaluation is also subject to the shared finite work counter, so valid inputs may return inconclusive when their tolerance or work requirement cannot be certified.