Multi-Electronic Atoms

1 · Kohn–Sham mean-field reduction

Kohn–Sham DFT replaces the interacting N-electron problem with N independent single-particle equations in a self-consistent effective potential V_eff. The equations must be solved iteratively (SCF loop) because V_eff depends on the density ρ, which depends on the orbitals φ_i. SPARC-atomSFE converges this loop using spectral finite elements with LDA_SVWN exchange-correlation.

2 · The effective potential — three contributions

V_nuc is exact and analytic. V_H is computed self-consistently on the radial grid (Poisson equation). V_xc is the only approximation: LDA_SVWN treats the exchange-correlation energy density as that of a uniform electron gas at the local density ρ(r) — accurate for core shells, slightly over-delocalizing for valence orbitals (the standard LDA self-interaction error, which is why the LDA valence radii are larger than Hartree-Fock Clementi-Raimondi values).

3 · Spherical symmetry — radial reduction and u(r)

The angular factor Y_ℓ^m separates exactly for a spherical V_eff(r). The substitution u = r·R eliminates the first-derivative term from the radial equation, turning it into a standard 1D Schrödinger form — the same shape as a particle in a potential, easy to solve numerically with finite elements or Numerov integration. The ℓ(ℓ+1)/(2r²) term is the centrifugal barrier: it pushes electrons with higher angular momentum away from the nucleus (p and d orbitals have no probability at r = 0).

4 · Why u(r)² is the probability we integrate

u(r)² already folds in the r² Jacobian of spherical coordinates — the same reason the hydrogen tables use [r·R(r)]² as their radial weight. The table row for a given (Z, n, ℓ) is therefore directly proportional to the probability of finding that subshell's electron at radius r; integrating it gives the inverse-CDF used for sampling. No extra r² factor is needed at draw time.

5 · The log-uniform radial grid

A log-uniform grid has a constant ratio between consecutive steps — Δr/r is constant — instead of a constant absolute step Δr. This is essential: the wavefunction oscillates on a scale of ~10⁻³ Bohr near a heavy nucleus (where u(r) rises sharply from zero), and decays exponentially over tens of Bohr in the valence region. A linear grid would either miss the core or waste 99% of its points on the tail. At 128 device points, the log-spacing keeps the mode radius accurate to better than 0.1% (measured against the 2001-point reference).

1 · Ground-state electron configuration

Madelung's rule gives the right filling order for most elements — it is a consequence of orbital energies being set more by n+ℓ than by n alone, due to electron-electron repulsion pushing inner shells inward. Known exceptions (Cr: [Ar]3d⁵4s¹, Cu: [Ar]3d¹⁰4s¹, and ~20 others for lanthanides and actinides) are substituted from a hard-coded table. Each subshell holds at most 2(2ℓ+1) electrons.

2 · Unsöld's theorem → full shells are isotropic

Unsöld's theorem (1927) is an exact result of spherical harmonics: when all 2ℓ+1 magnetic quantum numbers are equally occupied (as in a full subshell), the total angular density is a constant. The neon 2p⁶ shell is exactly spherical; the nitrogen 2p³ half-shell (one electron per m ∈ {−1, 0, +1}) is also spherical by the same argument — the three real spherical harmonics p_x, p_y, p_z together cover all three spatial directions equally. Carbon's 2p² is not spherical and must use Hund's rule instead.

3 · Hund's rule for partial subshells

Hund's first rule (maximum multiplicity) says electrons in a partially filled subshell prefer parallel spins, which forces them onto different orbitals. In practice this means each occupied m gets its own angular density profile — p_x, p_y, p_z for p², p³; d_z², d_xz, d_yz… for d³, etc. Each such orbital is sampled by the standard orbital sampler with the (ℓ, m)'s angular tables, not a uniform sphere. Without this, carbon's 2p² would incorrectly render as a sphere.

4 · Split the total point budget

N (kAtomNumPoints) is the total point budget for this element. Each drawing group's weight is its electron count — a doubly-occupied full 2p subshell (weight 6) gets three times as many points as a singly-occupied 2s (weight 2). The largest-remainder method prevents the sum from drifting off by one or two points due to floor truncation.

5 · Pick each group's radial source

The table path uses the all-electron Kohn-Sham wavefunction; the fallback uses the hydrogen analytic formula with a rescaled nuclear charge. Both produce an R(r) — the rest of the pipeline is identical from this point.

6 · Build the radial inverse-CDF

Both branches integrate the same physical density u(r)². The oriented branch resamples on a linear grid first because the hydrogen angular sampler (which it shares) expects a linear-grid radial CDF — an asymmetry kept deliberately from the Python model. The isotropic branch stays on the log-grid and integrates directly for higher accuracy near r = 0.

7 · Sample this group's N_g points

The choice between isotropic and oriented was made back in step 3; both arrive here with a radial CDF and convert the same way to Cartesian. Every draw is three uniform randoms + three O(1) table lookups — no rejection loop, no per-point search, deterministic cost.

8 · Identify the valence subshell

The P90 radius of each subshell's sampled points is a robust proxy for its spatial extent. The outermost subshell by this metric is used as the valence: it drives screen scaling (step 9) and gets the brighter color (step below). Note that for transition metals, 3d lies interior to 4s (Fe: 3d mode ≈ 43 pm, 4s mode ≈ 240 pm) — the P90 correctly identifies 4s as valence even though 3d has higher ℓ.

9 · Color by shell, emphasize valence

Points are emitted in contiguous per-subshell runs (AtomSubshellRange), so the lerp is applied once per subshell, not per point. Without the dim/brighten split, the valence subshell — often a small fraction of the total point budget — would be nearly invisible against the much larger core shells.

10 · Scale calibration and screen fit

The κ_Z factor corrects the LDA over-delocalization: it is the ratio of the Clementi-Raimondi literature valence radius to the LDA table's valence mode radius, computed once offline by atom_size_calib_gen.py. The subsequent divide by r_ref^outer then normalises all elements to the same on-screen target radius — like the orbital pipeline's P90 auto-zoom — so lithium and uranium are both legible on the 240 px display, at the cost of hiding the real size trend between elements.

1 · Clementi–Raimondi table (Z ≤ 54, primary source)

Clementi-Raimondi Z_eff values are obtained by fitting the hydrogen-like wavefunction to the actual Hartree-Fock SCF orbital — they are therefore the best possible nuclear charge to use with the analytic hydrogen formula. Their advantage over Slater's rules is significant for the p-block (error drops from ~5% to <3% for period 2 p states). They cover all occupied subshells of elements 1 to 54. Beyond that, no published 3-parameter fit exists in this project's scope, so Slater's rules take over.

2 · Slater's shielding rules (Z > 54)

Slater's rules (1930) capture the dominant physics of shielding: inner-core electrons (two or more groups below) screen the nuclear charge nearly completely, while electrons in the same subshell only partially screen each other (0.35, not 1.00), because they occupy the same spatial region rather than surrounding the nucleus. The d/f rule (1.00 for one group below) is a Slater-specific choice that slightly overestimates shielding for d electrons compared to the Clementi-Raimondi fit, but is simple and close enough for the fallback regime.

3 · The n* correction for outer shells (n ≥ 4)

The hydrogen radial wavefunction's mode radius scales as n² / Z. But Slater's Z_eff was derived assuming the hydrogenic n is the real principal quantum number. For valence electrons in period 4 and beyond, the wavefunction must be orthogonal to the filled inner shells — this pushes the outer lobe inward. Slater introduced an "effective quantum number" n* < n to account for this: replacing n with n* in the mode-radius formula gives a smaller (more physical) radius. The code applies this as a Z_eff rescaling (equivalent to n → n*) to keep the radial engine unchanged.

4 · The variable substitution r → Z_eff · r

The substitution contracts all distances by a factor of Z_eff — the orbital shape (number of radial nodes, angular lobes) is unchanged, only the scale shrinks. In code this means passing z_eff as a multiplicative factor into init_orbital_sampler() or buildRadialSamplerRuntime(), which applies it when evaluating R(r) — the function orbitals.py itself is never modified and remains the single validated source of the hydrogenic math. One Z_eff per subshell (not per individual electron) is the standard Slater mean-field approximation.

Neither model includes relativistic effects (Dirac equation). Relevant for Z ≳ 55 — the contracted 6s orbital of gold responsible for its color and Au–Au bond strength cannot be reproduced without a relativistic treatment. Noted as a known limit; not planned for this project's scope.