Radial functions: an experiment against Hartree–Fock

Question. The app draws and integrates valence orbitals with two different radial functions: Slater's single exponents in the orbital layer (π overlaps, Löwdin densities, the π maps) and effective-charge hydrogenic functions in the 3D viewer. Which one is closer to a real atomic orbital, and does the difference matter for what the app shows?

Forecast. Written down before any number was computed: FORECAST.md (10 Oct 2026, 10:45).

#The reference, and how it was checked

There are no published orbitals I could download here, so I computed them: Hartree–Fock for each free atom (atomic_hf.py), using Slater's average-of-configuration energy (1960) in a large even-tempered Gaussian basis (21–40 functions per angular momentum).

• Closed shells (Be, Ne, Mg, Ar, Kr, Xe) reproduce the numerical Hartree–Fock limit to 2 × 10⁻⁶ – 4 × 10⁻⁴ hartree. The Ne 2p ⟨r⟩ comes out at 0.965 bohr, the tabulated value.
• Open shells: from the same orbitals, the ground-term energy (³P, ⁴S or ²P) lands 0.01–1.6 mEh above the published Hartree–Fock limit, never below, as it must for orbitals not optimised for that term:

A first attempt with fractionally occupied open shells failed this check by 150–550 mEh: it lets an electron repel itself, which swells the open shell (carbon 2p ⟨r⟩ came out 23 % too large). Those numbers were discarded. That failure is why the check is in the script.

#Results

1 · The atomic functions. Overlap with Hartree–Fock (1 = identical), and in brackets the error in ⟨r⟩:

2 · What the orbital layer integrates. Two-centre overlaps of p orbitals with a carbon 2p at typical bond lengths, as a ratio to the Hartree–Fock value (both atoms described by the same model):

3 · What the viewer draws. sp^λ hybrids, isosurface at |ψ| = 0.05 a.u.: the fraction of the Hartree–Fock lobe volume that each model reproduces (intersection over union, positive / negative lobe):

#Forecast against outcome

#What it means for the app

1. Not affected: Hückel energies, charges, bond orders and the three pKa descriptors (Q_site, V_probe, ΔE_deloc). They never use a radial function.
2. Affected, orbital layer: the overlap matrix S and the Löwdin matrix M = S^-1/2 behind the π density, interference and orbital maps. Slater's 2p is 11–12 % too compact and the π overlaps are 26–38 % too small. How much that changes the maps themselves I have not measured yet.
3. Affected, 3D viewer: for Br and I the viewer uses Z_eff = 6.10, copied from chlorine, with n = 4 and 5. The result has essentially nothing in common with the real 4p and 5p (overlap 0.08 and 0.003). That is a bug. For P, S and Cl the hydrogenic 3s puts its nodes in the wrong place (overlap 0.53–0.58).
4. Affected, the viewer's isovalue rule: the viewer draws each surface at a fraction of max |ψ| on its grid, and the grid includes the nucleus. A realistic s function is largest at the nucleus: in a Hartree–Fock hybrid, |ψ| at the nucleus is 1.5–7.3 times the lobe maximum. Feed it realistic orbitals with the present rule and the nucleus would set the isovalue, shrinking the lobes. An absolute isovalue (|ψ| = 0.05 a.u., the usual choice in quantum chemistry) avoids that.

#Recommendation

Use one radial table for both the orbital layer and the viewer: the compact fits (2–6 Slater terms per orbital, overlap with Hartree–Fock ≥ 0.99992; parameters in results.json under compact_terms).

• The orbital layer keeps its quadrature: an overlap of two sums is the sum of the pairwise Slater overlaps it already computes.
• The viewer evaluates Σ c_k r^(N_k−1) e^(−ζ_k r) instead of one Laguerre polynomial: a few exponentials per grid point.
• At the same time, switch the viewer to an absolute isovalue.

The idea is old. Clementi and Raimondi (1963) refitted single exponents to Hartree–Fock atoms; Clementi and Roetti (1974) tabulated Hartree–Fock orbitals as short Slater sums. This experiment redoes it for the elements the app uses and measures what it changes on screen.

#Limits of this experiment

• The reference is the free atom. Orbitals contract somewhat in molecules; minimal basis sets optimised on molecules go the other way from these fits (Hehre, Stewart and Pople 1969 used ζ = 1.72 for carbon, more compact than Slater's 1.625). Free-atom Hartree–Fock is a clean, reproducible definition of an atomic orbital, not the last word on orbitals inside a molecule.
• Average-of-configuration orbitals, not term-specific ones (differences of order 1 % in ⟨r⟩).
• Non-relativistic: for iodine the 5p is slightly affected by relativity (spin–orbit splitting is ignored).

#What happened next

Version 16.1 adopted the recommendation. make_radial_table.py fits the valence ns and np orbitals of 16 elements (H, B, C, N, O, F, Na, Mg, Si, P, S, Cl, K, Ca, Br, I; every overlap with Hartree–Fock ≥ 0.9999) into qo_chem/data/radial_fits.json. The orbital layer's overlaps, the π maps and the 3D viewer's surfaces now read that one table, and the viewer draws every surface at a fixed |ψ| in atomic units. The hydrogenic functions, including the wrong bromine and iodine charges, are gone. The pKa descriptors did not move, as predicted: they never used a radial function.

#Reproduce

python tools/radial/atomic_hf.py          # validation against the Hartree–Fock limit
python tools/radial/radial_benchmark.py   # about 2½ minutes on one core: results.json, curves.npz
python tools/radial/plot_results.py       # the three figures
python tools/radial/make_radial_table.py  # the table the app uses (about 3 minutes)