Notebooks · lib/hf/iv_ints

Integrals with guaranteed bounds, two ways

lib/hf computes every integral by more than one algorithm and requires the results to agree. Agreement to 10⁻¹² says the algorithms are probably right; it does not say the common answer is the integral. This module computes intervals that are guaranteed to contain the exact integral, given the basis as data (exponents, raw coefficients and centers, each an exact double), using only the outward-rounded arithmetic of lib/hf/interval and the Boys enclosures of lib/hf/iv_boys.

It does this for s and p primitives, by two algorithms that share only the Boys enclosure:

McMurchie-Davidson (iv_md_*): Hermite expansion coefficients E(i, j, t) per axis and the Hermite Coulomb integrals R(t, u, v), as in lib/hf/hermite, lib/hf/md_onebody and lib/hf/md_eri.
Obara-Saika (iv_os_*): one-dimensional overlap and kinetic recursions, the vertical recursion on auxiliary integrals [a|b]^(m) for nuclear attraction and [e0|f0]^(m) for electron repulsion, and the horizontal transfer (a, b+1_i| = (a+1_i, b| + (A - B)_i (a, b|, as in lib/hf/os_onebody and lib/hf/hgp_eri.

Two enclosures of the same number must overlap. If they do not, one of the two algorithms is wrong, or the arithmetic under them is. Where they overlap, the exact value is in the common part.

What is enclosed. The unit-normalized Cartesian integral of lib/hf. The contracted coefficients are enclosed too (iv_cn): Nr(a, l) = (2a/pi)^(3/4) (4a)^(l/2) from square roots, and the unit-self-overlap rescale from the interval self-overlap. The point codes round those coefficients once; this module does not.

Width. For water STO-3G, measured in the tests below and in scf/iv_bridge: most enclosures of S, T, V and the sampled ERIs are about 2 × 10⁻¹⁴ wide relative to their larger end, and the widest is 5.0 × 10⁻¹², at the kinetic energy between O 1s and an H 1s (-0.00256), a sum whose terms cancel. An integral that is 0 by symmetry gets an enclosure around 0 at most 4 × 10⁻¹⁵ wide. Interval arithmetic treats a quantity used twice as two independent ones, and that costs little here because the recursions are short at l <= 1.

PySCF. The enclosures exclude some of PySCF's values: V between O pz and O 2s is 8.3 × 10⁻¹⁵ from the 40-digit value, which both enclosures hold and PySCF's does not. So PySCF 2.5.0 is off there by more than these widths.

Layout. An interval table is flat: entry k is [T[2k], T[2k+1]]. One-electron members return [n, n, 2]; ERI members return [m, 2] for a list of packed-8 indices. The sampled members use IV_ERI_NSAMPLE distinct indices drawn from IV_ERI_SEED.

Used with import hf.iv_ints.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 56–1119

from hf.iv_boys import iv_boys
from hf.interval import (
    iv_add, iv_contains, iv_div, iv_exp, iv_get, iv_mul, iv_neg, iv_outside, iv_overlaps, iv_pi,
    iv_powi, iv_pt, iv_put, iv_relw, iv_scale, iv_sqr, iv_sqrt, iv_sub, iv_width,
)
from hf.basis import (
    AO_ATOM, AO_LX, AO_LY, AO_LZ, AO_SH, AT_CHG, AT_X, AT_Y, AT_Z, bas_ao_table, bas_contract,
    bas_nao_cart, BAS_NDIM, bas_odfact, PR_COEF, PR_EXP, SH_ATOM, SH_L, SH_NPRIM, SH_P0,
)
from hf.eri_store import es_index8, es_len8, es_pair, es_rng
from hfref.water_sto3g import (
    watersto3g_atoms, watersto3g_eri8_cart, WATERSTO3G_ERI8_CART_LEN, WATERSTO3G_NAO_CART,
    WATERSTO3G_NPRIMTOT, watersto3g_prims, watersto3g_shells, watersto3g_v_unit,
)

IV_S = 0
IV_T = 1
IV_V = 2
IV_ERI_SEED = 20260923
IV_ERI_NSAMPLE = 100
# S, T and V: the integrals one AO pair gives, in the slots above
IV_ONE_LEN = 3
# where the E tables start in a row of iv_pair_data: after p, P, the
# coefficient, A - B and K_AB, 11 intervals of 2 doubles
IV_PD_E0 = 22

# ---- basis data as intervals --------------------------------------------------------

#| The contracted coefficients as intervals, flat over primitives: the raw
#| coefficient times `Nr(a, l)`, rescaled so each contracted shell has unit
#| self-overlap.
def iv_cn(shells: Arr, prims: Arr) {
    CN = zeros(2 * len(prims))
    for s to len(shells) {
        l = shells[s, SH_L]
        p0 = round(shells[s, SH_P0])
        np = round(shells[s, SH_NPRIM])
        ss = [0.0, 0.0]
        for i to np {
            for j to np {
                a = prims[p0 + i, PR_EXP]
                b = prims[p0 + j, PR_EXP]
                # overlap of two unit primitives: q^(l + 3/2), q = 2 sqrt(ab) / (a + b)
                q = iv_div(iv_scale(iv_sqrt(iv_mul(iv_pt(a), iv_pt(b))), 2.0), iv_add(iv_pt(a), iv_pt(b)))
                ov = iv_mul(iv_powi(q, l + 1), iv_sqrt(q))
                ss = iv_add(ss, iv_scale(iv_scale(ov, prims[p0 + i, PR_COEF]), prims[p0 + j, PR_COEF]))
            }
        }
        scale = iv_div([1.0, 1.0], iv_sqrt(ss))
        for i to np {
            a = prims[p0 + i, PR_EXP]
            y = iv_div(iv_pt(2.0 * a), iv_pi())
            nr = iv_mul(iv_mul(iv_sqrt(y), iv_sqrt(iv_sqrt(y))), iv_powi(iv_sqrt(iv_pt(4.0 * a)), l))
            iv_put(CN, p0 + i, iv_scale(iv_mul(scale, nr), prims[p0 + i, PR_COEF]))
        }
    }
    return CN
}

#| The component factor `Nc = 1 / sqrt((2lx-1)!! (2ly-1)!! (2lz-1)!!)`.

from hf.iv_boys import iv_boys
from hf.interval import (
    iv_add, iv_contains, iv_div, iv_exp, iv_get, iv_mul, iv_neg, iv_outside, iv_overlaps, iv_pi,
    iv_powi, iv_pt, iv_put, iv_relw, iv_scale, iv_sqr, iv_sqrt, iv_sub, iv_width,
)
from hf.basis import (
    AO_ATOM, AO_LX, AO_LY, AO_LZ, AO_SH, AT_CHG, AT_X, AT_Y, AT_Z, bas_ao_table, bas_contract,
    bas_nao_cart, BAS_NDIM, bas_odfact, PR_COEF, PR_EXP, SH_ATOM, SH_L, SH_NPRIM, SH_P0,
)
from hf.eri_store import es_index8, es_len8, es_pair, es_rng
from hfref.water_sto3g import (
    watersto3g_atoms, watersto3g_eri8_cart, WATERSTO3G_ERI8_CART_LEN, WATERSTO3G_NAO_CART,
    WATERSTO3G_NPRIMTOT, watersto3g_prims, watersto3g_shells, watersto3g_v_unit,
)

IV_S = 0
IV_T = 1
IV_V = 2
IV_ERI_SEED = 20260923
IV_ERI_NSAMPLE = 100
# S, T and V: the integrals one AO pair gives, in the slots above
IV_ONE_LEN = 3
# where the E tables start in a row of iv_pair_data: after p, P, the
# coefficient, A - B and K_AB, 11 intervals of 2 doubles
IV_PD_E0 = 22

# ---- basis data as intervals --------------------------------------------------------

#| The contracted coefficients as intervals, flat over primitives: the raw
#| coefficient times `Nr(a, l)`, rescaled so each contracted shell has unit
#| self-overlap.
def iv_cn(shells: Arr, prims: Arr) {
    CN = zeros(2 * len(prims))
    for s to len(shells) {
        l = shells[s, SH_L]
        p0 = round(shells[s, SH_P0])
        np = round(shells[s, SH_NPRIM])
        ss = [0.0, 0.0]
        for i to np {
            for j to np {
                a = prims[p0 + i, PR_EXP]
                b = prims[p0 + j, PR_EXP]
                # overlap of two unit primitives: q^(l + 3/2), q = 2 sqrt(ab) / (a + b)
                q = iv_div(iv_scale(iv_sqrt(iv_mul(iv_pt(a), iv_pt(b))), 2.0), iv_add(iv_pt(a), iv_pt(b)))
                ov = iv_mul(iv_powi(q, l + 1), iv_sqrt(q))
                ss = iv_add(ss, iv_scale(iv_scale(ov, prims[p0 + i, PR_COEF]), prims[p0 + j, PR_COEF]))
            }
        }
        scale = iv_div([1.0, 1.0], iv_sqrt(ss))
        for i to np {
            a = prims[p0 + i, PR_EXP]
            y = iv_div(iv_pt(2.0 * a), iv_pi())
            nr = iv_mul(iv_mul(iv_sqrt(y), iv_sqrt(iv_sqrt(y))), iv_powi(iv_sqrt(iv_pt(4.0 * a)), l))
            iv_put(CN, p0 + i, iv_scale(iv_mul(scale, nr), prims[p0 + i, PR_COEF]))
        }
    }
    return CN
}

#| The component factor `Nc = 1 / sqrt((2lx-1)!! (2ly-1)!! (2lz-1)!!)`.