3.1. Numerics. Integrals O ia (eq. 6) were evaluated by numerical quadrature using a single-center Euler-MacLaurin-Lebedev (EML) grid, with N r = 300 radial points and N Ω = 302 angular points per atom, except in the calculations reported in Section 4 .1 where extremely diffuse basis functions were employed. There, a denser grid (N r = 400, N Ω = 434) was used to compute the O ia . Either grid is considerably denser than those that are typically used to integrate the exchange-correlation functional. For the latter purpose, the SG-1 quadrature grid was used except in Sections 4 .1 and 4 .4 where the SG-3 grid 92 was used instead. For most calculations, the integral drop tolerance was set at τ ints = 10 -12 E h and the shell-pair drop tolerance was set at τ shlpr = 10 -12 a.u., although in Sections 4 .1 and 4 .4 the tighter values τ ints = 10 -14 E h and τ shlpr = 10 -16 a.u. were used instead. (See Ref. 94 for a discussion of appropriate numerical thresholds for diffuse basis functions.) The self-consistent field (SCF) convergence criterion was set to τ SCF = 10 -8 E h in the norm of the orbital gradient an eq. 2 was considered to be converged when the maximum element of the residual vector fell below τ CIS = 10 -6 E h .