This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org

Notebooks · lib/hf/boys

The Boys function, four ways, held against mpmath

Every Gaussian integral with a 1/r in it, nuclear attraction and electron repulsion alike, ends in the Boys function

F_n(x) = int_0^1 t^(2n) exp(-x t^2) dt

at orders n = 0 up to the total angular momentum of the integral, and at an argument x that grows with the exponents and the distance between centers. This module covers n <= 16 (an (gg|gg) quartet) and 0 <= x <= 200.

It uses exp, sqrt, cos and arithmetic, and no erf or gamma. Those two differ between the interpreter and C in the last bits, and a function called once per primitive quartet should not be where the two engines part.

Every member has the signature (nmax, x) -> F[0..nmax], because an integral wants all orders at once and they come cheaply from one another through

F_(n-1)(x) = (2x F_n(x) + e^-x) / (2n - 1)       (downward)
F_(n+1)(x) = ((2n + 1) F_n(x) - e^-x) / (2x)     (upward)

Downward adds two positive numbers, so it never loses digits. Upward subtracts, and loses digits when e^-x is close to (2n + 1) F_n, which is small x.

The members of boys_methods():

  • boys_series_down: the series F_n = e^-x sum_k (2x)^k / ((2n+1)(2n+3)...(2n+2k+1)) at the top order, then downward. Every term is positive. The terms grow until k is near x, so at x = 200 the sum needs about 330 of them.
  • boys_asym_up: at or above the switch x_n = boys_switch(nmax), F_0 = sqrt(pi/x)/2 (the erfc tail it drops is at most 2^-53 of F_0 there) and then upward, keeping the e^-x term. Below the switch it is the series.
  • boys_quadrature: 20-point Gauss-Legendre on [0, 1], with the panel count doubled until two passes agree. This is the only member that shares no code with the series, so it is the one that can catch a bug in it.
  • boys_taylor: a Taylor step from the nearest node of a grid with spacing 0.1. Production callers build the grid once with boys_grid_new inside their own integral def (so it is freed when that def returns) and call boys_grid_eval. The registry member builds a one-node grid on the spot, which is slow but runs the same evaluation code. The node values come from the series, so this member is not independent of boys_series_down.

boys_values is the production entry: the series below the switch and the asymptotic form above it, which is boys_asym_up.

Where the switch goes. The old code in lib/gto_md switched at 50 to the pure asymptotic form F_n = (2n-1)!!/2^(n+1) sqrt(pi/x^(2n+1)), which drops every term carrying e^-x. Relative to F_n those terms grow with n, so the error does too: 6.27 × 10−11 relative at n = 12 and 1.15 × 10−8 at n = 16, both at x = 50. A per-order switch would fix that, but keeping e^-x in the upward recursion fixes it better. Then the only thing dropped is the erfc tail of F_0, and upward recursion carries a relative error in F_0 into every order unchanged while 2x > 2n + 1, which holds for n <= 16 above 16.5. The largest relative error against mpmath, over orders 0 to 16 for sqrt(pi/x)/2 plus upward recursion, and at n = 16 for the pure closed form:

x                25        30        33.3      36        50        70
upward, e^-x     1.6e-12   9.8e-15   5.8e-16   8.5e-16   2.0e-16   2.2e-16
closed form      3.0e-2    2.8e-3    4.7e-4    1.0e-4    1.1e-8    3.9e-15

The closed form would need a switch near 70 at n = 16.

The switch per order. boys_switch(n) gives x_n for each order, and boys_asym_up uses the one for its top order, since it returns every order up to it. x_n is the smallest multiple of 0.1 at which two things hold for every order m <= n: the dropped erfc tail, carried up to F_m by the recursion, is at most 2^-53 of F_m, and 2x > 2m + 1, so the upward step does not amplify rounding. tools/boys_reference.py recomputes each x_n from mpmath and refuses to write the table if the list in boys_switch differs. The tail's share of F_m is erfc(sqrt(x)) times a factor that moves only in the fourth digit between m = 0 and m = 16 (at x = 30, 9.4857 × 10−15 against 9.5121 × 10−15), and the second condition only binds below 16.5. So the tail criterion gives x_n = 34.4 for every order up to 16. The switch is per order in its interface and in the test that fixes it, and constant in value over this range.

Worst relative error against the table, over every tabulated x and every top order up to 16: boys_values 1.7 × 10−15, boys_series_down 3.4 × 10−15, boys_quadrature 2.4 × 10−15, boys_taylor and a caller's grid 1.6e-15.

The mpmath table in lib/hfref/boys (written by tools/boys_reference.py) has points at x_n and 0.1 either side of it for every order, and on both sides of the old switch at 50. The tests below hold every member against it.

Used with import hf.boys.

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 91–328