#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
periodic-table-shape.py
Science Journaling Club - analysis for "Why the Periodic Table Is That Peculiar Shape"

Two independent calculations, both from first principles, both checked against
published measurement.

PART A. Build the shape of the periodic table out of quantum numbers alone.
        - enumerate every allowed orbital (n, l) with l < n
        - give each one 2(2l+1) seats, which is where the block widths come from
        - order them by the Madelung rule (increasing n+l, ties broken by lower n)
        - read off period lengths, block boundaries and the noble-gas atomic numbers
        - compare the predicted ground-state configuration of every element
          Z = 1..103 against the NIST-tabulated real one and list every violation

PART B. Effective nuclear charge from Slater's rules, and what it predicts.
        - Zeff across periods 2 and 3
        - first ionisation energy as a difference of Slater total energies,
          compared with NIST measured values
        - atomic radius as k * n*^2 / Zeff with one fitted scale constant,
          compared with Slater's (1964) empirical radii and Cordero's (2008)
          covalent radii
        - every discrepancy reported, not smoothed over

VALIDATION (printed side by side with the analytic / accepted answer):
        - hydrogen: Slater gives Zeff = 1 exactly, so the predicted ionisation
          energy must equal the analytic Rydberg, 13.6057 eV. NIST: 13.5984 eV.
        - He+ is a one-electron ion, so the analytic answer is 4 Ry = 54.4228 eV.
          NIST second ionisation energy of helium: 54.4178 eV.
        - the Madelung construction must reproduce the noble-gas atomic numbers
          2, 10, 18, 36, 54, 86, 118 with no input from chemistry at all.

Pure standard library. No third-party imports. Run:  python periodic-table-shape.py
"""

from __future__ import annotations

import math

RY_EV = 13.605693122994      # Rydberg energy in eV (CODATA 2018)
A0_PM = 52.917721067         # Bohr radius in picometres

SUB = {0: "s", 1: "p", 2: "d", 3: "f", 4: "g", 5: "h", 6: "i", 7: "k"}

SYM = (
    "H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni "
    "Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe "
    "Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg "
    "Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg "
    "Bh Hs Mt Ds Rg Cn Nh Fl Mc Lv Ts Og"
).split()


# --------------------------------------------------------------------------
# PART A1.  Every orbital that quantum mechanics allows, and how wide it is.
# --------------------------------------------------------------------------

def orbitals(n_max: int = 8):
    """All (n, l) with 0 <= l < n, each carrying 2(2l+1) electron seats."""
    out = []
    for n in range(1, n_max + 1):
        for l in range(0, n):
            out.append({"n": n, "l": l, "seats": 2 * (2 * l + 1),
                        "name": "%d%s" % (n, SUB[l])})
    return out


def madelung_order(orbs):
    """Increasing n + l; ties broken by the smaller n. Nothing else."""
    return sorted(orbs, key=lambda o: (o["n"] + o["l"], o["n"]))


# --------------------------------------------------------------------------
# PART A2.  Fill the orbitals in that order and read the table off the result.
# --------------------------------------------------------------------------

def build_table(n_max: int = 8, z_stop: int = 120):
    """Walk the Madelung sequence, assigning electrons one at a time."""
    seq = madelung_order(orbitals(n_max))
    rows, z = [], 0
    for o in seq:
        if z >= z_stop:
            break
        start = z + 1
        z += o["seats"]
        rows.append({"name": o["name"], "n": o["n"], "l": o["l"],
                     "seats": o["seats"], "nl": o["n"] + o["l"],
                     "z_from": start, "z_to": min(z, z_stop)})
    return rows


def periods_from(rows):
    """A period starts at each ns orbital. That is the whole definition."""
    periods, cur = [], None
    for r in rows:
        if r["l"] == 0:
            if cur:
                periods.append(cur)
            cur = {"n": r["n"], "orbitals": [], "length": 0}
        cur["orbitals"].append(r["name"])
        cur["length"] += r["seats"]
    if cur:
        periods.append(cur)
    return periods


def predicted_config(z: int, n_max: int = 8):
    """Madelung-rule ground-state configuration for element Z, as a list."""
    seq = madelung_order(orbitals(n_max))
    left, cfg = z, []
    for o in seq:
        if left <= 0:
            break
        take = min(left, o["seats"])
        cfg.append((o["name"], take))
        left -= take
    return cfg


def cfg_str(cfg):
    return " ".join("%s%d" % (nm, k) for nm, k in cfg)


def valence_str(cfg, keep=3):
    """Last few subshells, which is where every disagreement lives."""
    return " ".join("%s%d" % (nm, k) for nm, k in cfg[-keep:])


# --------------------------------------------------------------------------
# PART A3.  The elements that refuse.
#   Measured ground-state configurations, NIST Atomic Spectra Database v5.12
#   (Kramida et al. 2024).  Only the tail of each configuration is stored;
#   the script checks it against the tail the Madelung rule predicts.
# --------------------------------------------------------------------------

NIST_ANOMALIES = {
    24:  ("Cr", "3d5 4s1"),
    29:  ("Cu", "3d10 4s1"),
    41:  ("Nb", "4d4 5s1"),
    42:  ("Mo", "4d5 5s1"),
    44:  ("Ru", "4d7 5s1"),
    45:  ("Rh", "4d8 5s1"),
    46:  ("Pd", "4d10"),
    47:  ("Ag", "4d10 5s1"),
    57:  ("La", "5d1 6s2"),
    58:  ("Ce", "4f1 5d1 6s2"),
    64:  ("Gd", "4f7 5d1 6s2"),
    78:  ("Pt", "5d9 6s1"),
    79:  ("Au", "5d10 6s1"),
    89:  ("Ac", "6d1 7s2"),
    90:  ("Th", "6d2 7s2"),
    91:  ("Pa", "5f2 6d1 7s2"),
    92:  ("U",  "5f3 6d1 7s2"),
    93:  ("Np", "5f4 6d1 7s2"),
    96:  ("Cm", "5f7 6d1 7s2"),
    103: ("Lr", "7s2 7p1"),
}

# What the Madelung rule says the same twenty elements should look like,
# written out by hand so the script can confirm the disagreement is genuine
# rather than an artefact of how the strings are formatted.
MADELUNG_TAIL = {
    24:  "3d4 4s2",   29: "3d9 4s2",   41: "4d3 5s2",   42: "4d4 5s2",
    44:  "4d6 5s2",   45: "4d7 5s2",   46: "4d8 5s2",   47: "4d9 5s2",
    57:  "4f1 6s2",   58: "4f2 6s2",   64: "4f8 6s2",   78: "5d8 6s2",
    79:  "5d9 6s2",   89: "5f1 7s2",   90: "5f2 7s2",   91: "5f3 7s2",
    92:  "5f4 7s2",   93: "5f5 7s2",   96: "5f8 7s2",  103: "6d1 7s2",
}


def parse_tail(tail: str):
    """'3d5 4s1' -> {'3d': 5, '4s': 1}."""
    out = {}
    for tok in tail.split():
        out[tok[:2]] = int(tok[2:])
    return out


def check_prediction(z: int, tail: str) -> bool:
    """Confirm the hand-written Madelung tail really is what the generator
    produces, so the anomaly table cannot quietly drift away from the code."""
    gen = dict(predicted_config(z))
    want = parse_tail(tail)
    for name, k in want.items():
        if gen.get(name, 0) != k:
            return False
    seats = {o["name"]: o["seats"] for o in orbitals()}
    for name, k in gen.items():
        if name not in want and k != seats[name]:
            return False
    return True


def block_of(z: int, n_max: int = 8) -> str:
    """Which block the Madelung rule puts element Z in."""
    cfg = predicted_config(z, n_max)
    return cfg[-1][0][1]


# --------------------------------------------------------------------------
# PART B1.  Slater's rules (Slater 1930, Phys. Rev. 36, 57).
# --------------------------------------------------------------------------

SLATER_NSTAR = {1: 1.0, 2: 2.0, 3: 3.0, 4: 3.7, 5: 4.0, 6: 4.2, 7: 4.3}


def slater_groups(cfg):
    """Slater's bracketing: (1s)(2s,2p)(3s,3p)(3d)(4s,4p)(4d)(4f)(5s,5p)..."""
    groups = {}
    for name, k in cfg:
        n, sub = int(name[0]), name[1]
        key = (n, "sp") if sub in "sp" else (n, sub)
        groups[key] = groups.get(key, 0) + k
    return groups


def _sort_key(key):
    n, tag = key
    return (n, 0 if tag == "sp" else 1 + "dfg".index(tag))


def slater_zeff(z: int, cfg, key):
    """Screening constant and Zeff for one electron in the named group."""
    groups = slater_groups(cfg)
    n, tag = key
    s = 0.0
    for gk, cnt in groups.items():
        gn, gtag = gk
        if gk == key:
            same = cnt - 1                      # not counting the electron itself
            s += same * (0.30 if (n == 1 and tag == "sp") else 0.35)
        elif tag == "sp":
            if gn == n - 1:
                s += 0.85 * cnt
            elif gn < n - 1:
                s += 1.00 * cnt
            elif gn == n:                       # e.g. 3d screening a 3s electron
                s += 1.00 * cnt if gtag != "sp" and _sort_key(gk) < _sort_key(key) else 0.0
        else:                                   # a d or an f electron
            if _sort_key(gk) < _sort_key(key):
                s += 1.00 * cnt
    return z - s, s


def outer_key(cfg):
    name = cfg[-1][0]
    n, sub = int(name[0]), name[1]
    return (n, "sp") if sub in "sp" else (n, sub)


def slater_total_energy(z: int, cfg) -> float:
    """Sum of one-electron binding energies, in eV. Positive = bound."""
    groups = slater_groups(cfg)
    total = 0.0
    for key, cnt in groups.items():
        zeff, _ = slater_zeff(z, cfg, key)
        nstar = SLATER_NSTAR[key[0]]
        total += cnt * RY_EV * (zeff / nstar) ** 2
    return total


def cation_config(cfg):
    """Remove one electron from the highest-energy occupied subshell."""
    out = [[nm, k] for nm, k in cfg]
    out[-1][1] -= 1
    if out[-1][1] == 0:
        out.pop()
    return [(nm, k) for nm, k in out]


def slater_ionisation(z: int) -> float:
    cfg = predicted_config(z)
    return slater_total_energy(z, cfg) - slater_total_energy(z, cation_config(cfg))


# --------------------------------------------------------------------------
# PART B2.  Published measurements to test the predictions against.
# --------------------------------------------------------------------------

# First ionisation energies, eV. NIST Atomic Spectra Database v5.12.
IE_NIST = {
    1: 13.598434, 2: 24.587389, 3: 5.391715, 4: 9.322699, 5: 8.298019,
    6: 11.260288, 7: 14.534130, 8: 13.618055, 9: 17.422820, 10: 21.564541,
    11: 5.139077, 12: 7.646236, 13: 5.985769, 14: 8.151683, 15: 10.486686,
    16: 10.360001, 17: 12.967633, 18: 15.759611, 19: 4.340664, 20: 6.113158,
}
IE_HE_SECOND = 54.417765      # NIST, and analytically 4 Ry for a one-electron ion

# Empirical atomic radii, pm. Slater, J. C. (1964) J. Chem. Phys. 41, 3199.
R_SLATER64 = {
    3: 145, 4: 105, 5: 85, 6: 70, 7: 65, 8: 60, 9: 50,
    11: 180, 12: 150, 13: 125, 14: 110, 15: 100, 16: 100, 17: 100,
}

# Covalent radii, pm. Cordero et al. (2008) Dalton Trans. 2832.
R_CORDERO = {
    3: 128, 4: 96, 5: 84, 6: 76, 7: 71, 8: 66, 9: 57, 10: 58,
    11: 166, 12: 141, 13: 121, 14: 111, 15: 107, 16: 105, 17: 102, 18: 106,
}

# Zeff for the outermost electron from SCF wavefunctions.
# Clementi & Raimondi (1963) J. Chem. Phys. 38, 2686.
ZEFF_CLEMENTI = {
    3: 1.279, 4: 1.912, 5: 2.421, 6: 3.136, 7: 3.834, 8: 4.453,
    9: 5.100, 10: 5.758,
    11: 2.507, 12: 3.308, 13: 4.066, 14: 4.285, 15: 4.886, 16: 5.482,
    17: 6.116, 18: 6.764,
}


# --------------------------------------------------------------------------
# PART C.  Relativity, in the crudest honest form: a Bohr-model 1s electron.
# --------------------------------------------------------------------------

ALPHA = 7.2973525693e-3       # fine-structure constant, CODATA 2018


def relativistic_contraction(z: int):
    """v/c = Z*alpha for a 1s electron; radius scales as 1/gamma."""
    beta = z * ALPHA
    if beta >= 1:
        return beta, float("nan"), float("nan")
    gamma = 1.0 / math.sqrt(1.0 - beta * beta)
    return beta, gamma, 1.0 - 1.0 / gamma


# --------------------------------------------------------------------------
# Output
# --------------------------------------------------------------------------

def rule(ch="-", n=74):
    print(ch * n)


def main():
    print("=" * 74)
    print("WHY THE PERIODIC TABLE IS THAT PECULIAR SHAPE")
    print("Science Journaling Club - analysis/periodic-table-shape.py")
    print("=" * 74)

    # ---------------- A1: block widths --------------------------------
    print()
    print("PART A1. BLOCK WIDTHS FROM THE ALLOWED QUANTUM NUMBERS")
    rule()
    print("For a given l, m_l runs from -l to +l, which is 2l+1 values.")
    print("Each of those takes two electrons, one per spin. Seats = 2(2l+1).")
    print()
    print("  l   subshell   m_l values   seats = 2(2l+1)   block width")
    for l in range(4):
        print("  %d      %s          %2d            %2d               %2d"
              % (l, SUB[l], 2 * l + 1, 2 * (2 * l + 1), 2 * (2 * l + 1)))
    print()
    print("Those four numbers, 2 / 6 / 10 / 14, are the widths of the s, p, d")
    print("and f blocks. Nothing else in the table is free to choose them.")

    # ---------------- A2: filling order -------------------------------
    rows = build_table()
    print()
    print("PART A2. THE MADELUNG FILLING ORDER (increasing n+l, then increasing n)")
    rule()
    print("  order  orbital  n   l   n+l  seats   Z range")
    for i, r in enumerate(rows, 1):
        print("   %2d     %-4s   %d   %d    %d     %2d     %3d - %3d"
              % (i, r["name"], r["n"], r["l"], r["nl"], r["seats"],
                 r["z_from"], r["z_to"]))

    # ---------------- A3: period lengths ------------------------------
    periods = periods_from(rows)
    lengths = [p["length"] for p in periods]
    print()
    print("PART A3. PERIOD LENGTHS, READ STRAIGHT OFF THE FILLING ORDER")
    rule()
    print("  period   orbitals filled           length   ends at Z")
    running = 0
    ends = []
    for i, p in enumerate(periods, 1):
        running += p["length"]
        ends.append(running)
        print("    %d      %-24s  %3d       %3d"
              % (i, " ".join(p["orbitals"]), p["length"], running))
    print()
    print("  predicted period lengths : %s" % lengths[:7])
    print("  real periods 1-7         : [2, 8, 8, 18, 18, 32, 32]")
    print("  agreement                : %s"
          % ("EXACT" if lengths[:7] == [2, 8, 8, 18, 18, 32, 32] else "MISMATCH"))

    # ---------------- VALIDATION 1 ------------------------------------
    print()
    print("VALIDATION 1. NOBLE-GAS ATOMIC NUMBERS, PREDICTED vs OBSERVED")
    rule()
    noble = [2, 10, 18, 36, 54, 86, 118]
    noble_sym = ["He", "Ne", "Ar", "Kr", "Xe", "Rn", "Og"]
    print("  period   predicted Z at period end   observed noble gas   diff")
    for i in range(7):
        print("    %d              %3d                    %-3s %3d           %+d"
              % (i + 1, ends[i], noble_sym[i], noble[i], ends[i] - noble[i]))
    print()
    print("  The construction used no chemistry. It used l < n, 2(2l+1) seats,")
    print("  and the ordering n+l. It lands on all seven noble gases exactly,")
    print("  and the seven periods sum to %d, the number of known elements."
          % sum(lengths[:7]))

    # ---------------- A4: the generated shape -------------------------
    print()
    print("PART A4. THE SHAPE ITSELF, GENERATED FROM THE RULES ALONE")
    rule()
    blocks = {}
    for z in range(1, 119):
        blocks.setdefault(block_of(z), []).append(z)
    for b in "spdf":
        zs = blocks.get(b, [])
        print("  %s-block : %3d elements, %2d columns, %d rows"
              % (b, len(zs), {"s": 2, "p": 6, "d": 10, "f": 14}[b],
                 len(zs) // {"s": 2, "p": 6, "d": 10, "f": 14}[b]))
    print("  total    : %d" % sum(len(v) for v in blocks.values()))
    print()
    print("  Laid out honestly, with the f-block in line, the table is")
    print("  2 + 14 + 10 + 6 = 32 columns wide. The familiar 18-column table")
    print("  is that same table with 14 columns cut out and parked underneath.")

    # ---------------- A5: the violations ------------------------------
    print()
    print("PART A5. WHERE THE RULE BREAKS (Z = 1..103, against NIST v5.12)")
    rule()
    print("   Z  sym   Madelung predicts      NIST measures         block")
    n_bad = 0
    for z in sorted(NIST_ANOMALIES):
        sym, actual = NIST_ANOMALIES[z]
        pred = MADELUNG_TAIL[z]
        assert check_prediction(z, pred), (
            "generator and hand-written prediction disagree at Z=%d (%s)" % (z, pred))
        assert parse_tail(pred) != parse_tail(actual), (
            "Z=%d is not actually an anomaly" % z)
        n_bad += 1
        print("  %3d  %-3s  %-21s  %-20s  %s"
              % (z, sym, pred, actual, block_of(z)))
    print()
    print("  violations found      : %d out of 103 elements (%.1f%%)"
          % (n_bad, 100.0 * n_bad / 103))
    print("  obedient elements     : %d out of 103 (%.1f%%)"
          % (103 - n_bad, 100.0 * (103 - n_bad) / 103))
    by_block = {}
    for z in NIST_ANOMALIES:
        by_block[block_of(z)] = by_block.get(block_of(z), 0) + 1
    print("  violations by block   : %s"
          % ", ".join("%s: %d" % (b, by_block.get(b, 0)) for b in "spdf"))
    census = {}
    for z in range(1, 104):
        census[block_of(z)] = census.get(block_of(z), 0) + 1
    print("  block census, Z=1..103: %s"
          % ", ".join("%s: %d" % (b, census.get(b, 0)) for b in "spdf"))
    print("  Every one sits in the d-block or the f-block. The s-block and the")
    print("  p-block, %d elements between them, never break the rule once."
          % (census.get("s", 0) + census.get("p", 0)))
    print()
    single = sum(1 for z in NIST_ANOMALIES
                 if NIST_ANOMALIES[z][1].rstrip().endswith("s1"))
    print("  of the %d, how many move exactly one electron out of an ns shell: %d"
          % (n_bad, single))
    print("  Pd (Z=46) moves two, emptying 5s completely: [Kr] 4d10.")

    # ---------------- B1: Zeff ----------------------------------------
    print()
    print("PART B1. EFFECTIVE NUCLEAR CHARGE FROM SLATER'S RULES")
    rule()
    print("  Zeff = Z - S. For an ns or np electron, S counts 0.35 from each")
    print("  other electron in the same shell, 0.85 from each in the shell below,")
    print("  and 1.00 from everything deeper.")
    print()
    print("   Z  sym   config (valence)    S      Zeff    Zeff(SCF)   Slater - SCF")
    zeff_rows = {}
    for z in list(range(3, 11)) + list(range(11, 19)):
        cfg = predicted_config(z)
        key = outer_key(cfg)
        zeff, s = slater_zeff(z, cfg, key)
        zeff_rows[z] = zeff
        cl = ZEFF_CLEMENTI[z]
        print("  %3d  %-3s  %-18s %6.2f  %6.3f    %6.3f      %+7.3f"
              % (z, SYM[z - 1], valence_str(cfg, keep=2), s, zeff, cl, zeff - cl))
        if z == 10:
            print("  " + "." * 68)
    p2 = [zeff_rows[z] for z in range(3, 11)]
    p3 = [zeff_rows[z] for z in range(11, 19)]
    print()
    print("  period 2: Zeff runs %.2f -> %.2f, a rise of %.2f over 7 steps"
          % (p2[0], p2[-1], p2[-1] - p2[0]))
    print("            Slater step = %.3f per element, exactly and always"
          % ((p2[-1] - p2[0]) / 7))
    cl2 = [ZEFF_CLEMENTI[z] for z in range(3, 11)]
    print("            SCF step    = %.3f per element on average"
          % ((cl2[-1] - cl2[0]) / 7))
    print("  period 3: Zeff runs %.2f -> %.2f, same %.2f step"
          % (p3[0], p3[-1], (p3[-1] - p3[0]) / 7))
    cl3 = [ZEFF_CLEMENTI[z] for z in range(11, 19)]
    print("            SCF step    = %.3f per element on average"
          % ((cl3[-1] - cl3[0]) / 7))
    mad2 = sum(abs(zeff_rows[z] - ZEFF_CLEMENTI[z]) for z in range(3, 11)) / 8
    mad3 = sum(abs(zeff_rows[z] - ZEFF_CLEMENTI[z]) for z in range(11, 19)) / 8
    print("  mean |Slater - SCF|, period 2: %.3f    period 3: %.3f" % (mad2, mad3))
    print("  Slater runs high through period 2 and low through period 3, and the")
    print("  period-3 error shrinks from %+.3f at Na to %+.3f at Ar."
          % (zeff_rows[11] - ZEFF_CLEMENTI[11], zeff_rows[18] - ZEFF_CLEMENTI[18]))

    # ---------------- VALIDATION 2 ------------------------------------
    print()
    print("VALIDATION 2. ONE-ELECTRON LIMIT, WHERE THE ANSWER IS ANALYTIC")
    rule()
    h_pred = slater_ionisation(1)
    print("  quantity                        model        analytic / NIST     diff")
    print("  H  first ionisation energy    %8.4f eV     %8.4f eV      %+7.4f eV"
          % (h_pred, IE_NIST[1], h_pred - IE_NIST[1]))
    print("     (analytic hydrogenic value  %8.4f eV = 1 Ry)" % RY_EV)
    he2_pred = RY_EV * (2.0 ** 2) / 1.0 ** 2
    print("  He+ ionisation energy         %8.4f eV     %8.4f eV      %+7.4f eV"
          % (he2_pred, IE_HE_SECOND, he2_pred - IE_HE_SECOND))
    print("     (analytic hydrogenic value  %8.4f eV = 4 Ry)" % (4 * RY_EV))
    print()
    print("  With one electron there is nothing to screen, so Slater's rules")
    print("  collapse to the exact hydrogenic formula. Agreement to %.2f%% and"
          % (100 * abs(h_pred - IE_NIST[1]) / IE_NIST[1]))
    print("  %.3f%% is the arithmetic behaving, not the chemistry."
          % (100 * abs(he2_pred - IE_HE_SECOND) / IE_HE_SECOND))

    # ---------------- B2: ionisation energies -------------------------
    print()
    print("PART B2. PREDICTED vs MEASURED FIRST IONISATION ENERGY")
    rule()
    print("  IE = E(atom) - E(cation), both summed as sum_i Ry*(Zeff_i/n*_i)^2.")
    print()
    print("   Z  sym   predicted   measured    diff      error")
    errs = []
    pred_ie = {}
    for z in range(1, 21):
        p = slater_ionisation(z)
        m = IE_NIST[z]
        pred_ie[z] = p
        e = 100.0 * (p - m) / m
        errs.append(abs(e))
        print("  %3d  %-3s  %8.3f    %8.3f   %+7.3f   %+7.1f%%"
              % (z, SYM[z - 1], p, m, p - m, e))
        if z in (2, 10, 18):
            print("  " + "." * 60)
    print()
    print("  mean absolute error, Z = 1..20 : %.1f%%" % (sum(errs) / len(errs)))
    e2 = [100.0 * (pred_ie[z] - IE_NIST[z]) / IE_NIST[z] for z in range(3, 11)]
    print("  mean absolute error, period 2  : %.1f%%"
          % (sum(abs(x) for x in e2) / len(e2)))
    e3 = [100.0 * (pred_ie[z] - IE_NIST[z]) / IE_NIST[z] for z in range(11, 19)]
    print("  mean absolute error, period 3  : %.1f%%"
          % (sum(abs(x) for x in e3) / len(e3)))
    print("  period 3 is where the model falls apart: every prediction is too")
    print("  high, because Slater's n* = 3.0 for the third shell is too small")
    print("  and the 0.85 screening from the n = 2 shell is far too generous.")
    print("  worst single element           : %s at %+.1f%%"
          % (SYM[max(range(1, 21), key=lambda z: abs(100.0 * (pred_ie[z] - IE_NIST[z]) / IE_NIST[z])) - 1],
             max((100.0 * (pred_ie[z] - IE_NIST[z]) / IE_NIST[z] for z in range(1, 21)), key=abs)))
    print()
    print("  Direction of travel across period 2:")
    up_pred = sum(1 for z in range(4, 11) if pred_ie[z] > pred_ie[z - 1])
    up_meas = sum(1 for z in range(4, 11) if IE_NIST[z] > IE_NIST[z - 1])
    print("    steps where the model rises   : %d of 7" % up_pred)
    print("    steps where measurement rises : %d of 7" % up_meas)
    print("    the two measured reversals    : Be -> B  (%.3f -> %.3f eV)"
          % (IE_NIST[4], IE_NIST[5]))
    print("                                    N  -> O  (%.3f -> %.3f eV)"
          % (IE_NIST[7], IE_NIST[8]))
    print("    the model predicts both as rises, so it misses both. Slater's")
    print("    rules know nothing about which orbital an electron sits in, and")
    print("    both reversals are orbital effects: the 2p electron of boron is")
    print("    higher in energy than beryllium's 2s, and oxygen's fourth 2p")
    print("    electron has to pair up in an orbital that already holds one.")

    # ---------------- B3: radii ---------------------------------------
    print()
    print("PART B3. PREDICTED vs MEASURED ATOMIC RADIUS")
    rule()
    print("  Hydrogenic scaling: r proportional to n*^2 / Zeff, with a single")
    print("  scale constant k fitted once by least squares over Li..F and Na..Cl")
    print("  against Slater's (1964) empirical radii. One free parameter total.")
    print()
    zs = sorted(R_SLATER64)
    xs, ys = [], []
    for z in zs:
        cfg = predicted_config(z)
        zeff, _ = slater_zeff(z, cfg, outer_key(cfg))
        nstar = SLATER_NSTAR[outer_key(cfg)[0]]
        xs.append(nstar ** 2 / zeff)
        ys.append(float(R_SLATER64[z]))
    k = sum(x * y for x, y in zip(xs, ys)) / sum(x * x for x in xs)
    print("  fitted k = %.2f pm   (Bohr radius a0 = %.2f pm, ratio %.3f)"
          % (k, A0_PM, k / A0_PM))
    print()
    print("   Z  sym   n*^2/Zeff   predicted   Slater'64   Cordero'08   resid vs S64")
    resid = []
    for z in zs:
        cfg = predicted_config(z)
        zeff, _ = slater_zeff(z, cfg, outer_key(cfg))
        nstar = SLATER_NSTAR[outer_key(cfg)[0]]
        x = nstar ** 2 / zeff
        p = k * x
        r = p - R_SLATER64[z]
        resid.append(abs(r))
        print("  %3d  %-3s   %8.4f   %8.1f    %6d      %6s       %+7.1f"
              % (z, SYM[z - 1], x, p, R_SLATER64[z],
                 R_CORDERO.get(z, "-"), r))
        if z == 9:
            print("  " + "." * 66)
    print()
    print("  mean absolute residual vs Slater 1964 : %.1f pm" % (sum(resid) / len(resid)))
    print()

    def ratio(z_lo, z_hi, table):
        return table[z_lo] / table[z_hi]

    def pred_r(z):
        cfg = predicted_config(z)
        zeff, _ = slater_zeff(z, cfg, outer_key(cfg))
        return k * SLATER_NSTAR[outer_key(cfg)[0]] ** 2 / zeff

    print("  Contraction across a period, measured as r(first) / r(last):")
    print("    period 2, model      Li/F  = %.2f" % (pred_r(3) / pred_r(9)))
    print("    period 2, Slater'64  Li/F  = %.2f" % ratio(3, 9, R_SLATER64))
    print("    period 2, Cordero'08 Li/F  = %.2f" % ratio(3, 9, R_CORDERO))
    print("    period 3, model      Na/Cl = %.2f" % (pred_r(11) / pred_r(17)))
    print("    period 3, Slater'64  Na/Cl = %.2f" % ratio(11, 17, R_SLATER64))
    print("    period 3, Cordero'08 Na/Cl = %.2f" % ratio(11, 17, R_CORDERO))
    print()
    print("  The model gets the direction right everywhere and the amount wrong")
    print("  in period 3, where it shrinks the atoms about %.1f times harder than"
          % ((pred_r(11) / pred_r(17)) / ratio(11, 17, R_CORDERO)))
    print("  the covalent radii do. Covalent radii are contact distances between")
    print("  bonded atoms, and near the right-hand edge they stop shrinking almost")
    print("  entirely. A hydrogenic 1/Zeff law has no way to know that.")
    print()
    print("  Down a group, where the model does better:")
    for a, b in ((3, 11), (4, 12), (5, 13), (6, 14), (9, 17)):
        print("    %-2s -> %-2s   model %5.1f -> %5.1f pm   Slater'64 %3d -> %3d pm"
              % (SYM[a - 1], SYM[b - 1], pred_r(a), pred_r(b),
                 R_SLATER64[a], R_SLATER64[b]))

    # ---------------- C: relativity -----------------------------------
    print()
    print("PART C. THE CRUDE RELATIVISTIC ESTIMATE")
    rule()
    print("  A Bohr-model 1s electron orbits at v/c = Z*alpha. Its mass goes up")
    print("  by gamma, and the orbit radius goes down by the same factor.")
    print()
    print("   Z  sym    v/c      gamma    1s contraction")
    for z, s in ((1, "H"), (26, "Fe"), (47, "Ag"), (54, "Xe"), (79, "Au"),
                 (80, "Hg"), (82, "Pb"), (92, "U"), (118, "Og")):
        beta, gamma, c = relativistic_contraction(z)
        print("  %3d  %-3s  %6.4f   %7.4f      %5.1f%%" % (z, s, beta, gamma, 100 * c))
    print()
    b_ag, g_ag, c_ag = relativistic_contraction(47)
    b_au, g_au, c_au = relativistic_contraction(79)
    print("  silver -> gold, one row apart in the same group:")
    print("    contraction goes from %.1f%% to %.1f%%, a factor of %.1f"
          % (100 * c_ag, 100 * c_au, c_au / c_ag))
    print("  This is the back-of-an-envelope version. Proper four-component")
    print("  Dirac-Fock calculations give a 6s contraction for gold of the same")
    print("  order, and they are what the chemistry arguments actually rest on")
    print("  (Desclaux 1973; Pyykko & Desclaux 1979).")
    print()
    print("  Photon energies, for the colour argument:")
    for ev, what in ((2.4, "gold, 5d -> 6s absorption edge"),
                     (3.9, "silver, 4d -> 5s absorption edge")):
        print("    %.1f eV  ->  %5.0f nm   %s" % (ev, 1239.841984 / ev, what))
    print("    visible band runs about 380-750 nm, so gold's edge falls inside it")
    print("    and silver's falls in the ultraviolet, outside it.")

    # ---------------- headline numbers --------------------------------
    print()
    print("HEADLINE NUMBERS USED IN THE ARTICLE AND THE INTERACTIVE MODEL")
    rule("=")
    print("  block widths                       : 2 / 6 / 10 / 14")
    print("  period lengths                     : 2, 8, 8, 18, 18, 32, 32")
    print("  period ends                        : 2, 10, 18, 36, 54, 86, 118")
    print("  full-width table                   : 32 columns")
    print("  columns cut out and parked below   : 14")
    print("  Madelung violations, Z = 1..103    : %d (%.1f%%)"
          % (n_bad, 100.0 * n_bad / 103))
    print("  violations in the s and p blocks   : 0")
    print("  Slater Zeff step per element       : 0.650 exactly")
    print("  SCF Zeff step, period 2            : %.3f" % ((cl2[-1] - cl2[0]) / 7))
    print("  IE mean abs error, period 2        : %.1f%%"
          % (sum(abs(x) for x in e2) / len(e2)))
    print("  IE mean abs error, period 3        : %.1f%%"
          % (sum(abs(x) for x in e3) / len(e3)))
    print("  IE reversals in period 2, measured : 2 (Be->B, N->O)")
    print("  IE reversals the model predicts    : 0")
    print("  radius law, one fitted constant k  : %.2f pm" % k)
    print("  radius mean abs residual           : %.1f pm" % (sum(resid) / len(resid)))
    print("  1s contraction, Ag / Au / Og       : %.1f%% / %.1f%% / %.1f%%"
          % (100 * relativistic_contraction(47)[2],
             100 * relativistic_contraction(79)[2],
             100 * relativistic_contraction(118)[2]))
    print("  gold absorption edge               : 2.4 eV = %.0f nm" % (1239.841984 / 2.4))
    print("  silver absorption edge             : 3.9 eV = %.0f nm" % (1239.841984 / 3.9))

    print()
    print("=" * 74)
    print("END. Every number printed above comes from this file.")
    print("=" * 74)


if __name__ == "__main__":
    main()
