Tammes is my packing solver with the container replaced by the sphere surface: projection p/|p|, no wall term, Fibonacci-spiral starts. 0.9926 -> 0.9974.

-"""Baseline: Fibonacci-sphere start, then projected gradient descent on the soft-min energy
-sum (s/d_ij)^p with p stepped up 8 -> 64, keeping the best minimum distance seen. Lands a few
-percent below the records. Beat it."""
+"""Tammes: n points on the unit sphere, maximising the smallest pairwise distance.
 
+This is the same inflate-and-separate scheme I used for circle packing in a square, a
+polygon, a cube and a ball, with the container replaced by the sphere surface.  There the
+container entered as a projection; here the projection is simply p / |p|, and because every
+point lies on the boundary there is no wall term at all.  The objective is the minimum
+pairwise distance directly.
+
+Method.  Hold a target separation d, run Gauss-Seidel sweeps that push any pair closer than d
+apart along their chord and renormalise both, then ratchet d up whenever the achieved minimum
+improves and cut the step when it does not.  Starts are Fibonacci spirals at several offsets
+(near-optimal on a sphere and the natural analogue of the cropped lattice), plus jittered and
+random ones.  Then reinsertion: the pair realising the minimum is what caps D, so move one of
+the two to the emptiest point of the sphere and re-relax.
+"""
+
 import math
 import random
 import time
 
+_GOLDEN = math.pi * (3.0 - math.sqrt(5.0))
 
-def _fibonacci(n: int) -> list[list[float]]:
-    golden = math.pi * (3.0 - math.sqrt(5.0))
+
+def _fib(n, off):
     pts = []
     for i in range(n):
-        z = 1.0 - (2.0 * i + 1.0) / n
+        z = 1.0 - (2.0 * i + 2.0 * off) / n
+        if z > 1.0:
+            z = 1.0
+        elif z < -1.0:
+            z = -1.0
         r = math.sqrt(max(0.0, 1.0 - z * z))
-        pts.append([r * math.cos(golden * i), r * math.sin(golden * i), z])
+        a = _GOLDEN * i
+        pts.append([r * math.cos(a), r * math.sin(a), z])
     return pts
 
 
-def _min_dist(pts: list[list[float]]) -> float:
-    best = 4.0
-    for i in range(len(pts)):
-        xi, yi, zi = pts[i]
-        for j in range(i + 1, len(pts)):
-            xj, yj, zj = pts[j]
-            d2 = (xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2
-            if d2 < best:
-                best = d2
+def _rand(n, rng):
+    pts = []
+    for _ in range(n):
+        while True:
+            x = rng.gauss(0.0, 1.0)
+            y = rng.gauss(0.0, 1.0)
+            z = rng.gauss(0.0, 1.0)
+            s = math.sqrt(x * x + y * y + z * z)
+            if s > 1e-9:
+                pts.append([x / s, y / s, z / s])
+                break
+    return pts
+
+
+def _min_dist(pts, n):
+    best = 8.0
+    for i in range(n - 1):
+        a = pts[i]
+        ax = a[0]
+        ay = a[1]
+        az = a[2]
+        for j in range(i + 1, n):
+            b = pts[j]
+            dx = b[0] - ax
+            dy = b[1] - ay
+            dz = b[2] - az
+            s = dx * dx + dy * dy + dz * dz
+            if s < best:
+                best = s
     return math.sqrt(best)
 
 
-def _energy_and_forces(pts: list[list[float]], p: float, s: float) -> tuple[float, list[list[float]]]:
-    """E = sum (s/d)^p; the force on i from j is p (s/d)^p / d^2 * (x_i - x_j)."""
-    n = len(pts)
-    e = 0.0
-    f = [[0.0, 0.0, 0.0] for _ in range(n)]
-    for i in range(n):
-        xi, yi, zi = pts[i]
-        fi = f[i]
+def _min_pair(pts, n):
+    best = 8.0
+    bi, bj = 0, 1
+    for i in range(n - 1):
+        a = pts[i]
+        ax = a[0]
+        ay = a[1]
+        az = a[2]
         for j in range(i + 1, n):
-            xj, yj, zj = pts[j]
-            dx, dy, dz = xi - xj, yi - yj, zi - zj
-            d2 = dx * dx + dy * dy + dz * dz
-            t = (s * s / d2) ** (p / 2)
-            e += t
-            g = p * t / d2
-            fi[0] += dx * g; fi[1] += dy * g; fi[2] += dz * g
-            fj = f[j]
-            fj[0] -= dx * g; fj[1] -= dy * g; fj[2] -= dz * g
-    return e, f
+            b = pts[j]
+            dx = b[0] - ax
+            dy = b[1] - ay
+            dz = b[2] - az
+            s = dx * dx + dy * dy + dz * dz
+            if s < best:
+                best = s
+                bi, bj = i, j
+    return bi, bj
 
 
-def _step(pts: list[list[float]], f: list[list[float]], lr: float) -> list[list[float]]:
-    out = []
-    for (x, y, z), (fx, fy, fz) in zip(pts, f):
-        rad = fx * x + fy * y + fz * z
-        nx, ny, nz = x + lr * (fx - rad * x), y + lr * (fy - rad * y), z + lr * (fz - rad * z)
-        r = math.sqrt(nx * nx + ny * ny + nz * nz)
-        out.append([nx / r, ny / r, nz / r])
-    return out
+def _separate(pts, n, d, sweeps, rng):
+    d2 = d * d
+    for _ in range(sweeps):
+        worst = 0.0
+        for i in range(n):
+            a = pts[i]
+            ax = a[0]
+            ay = a[1]
+            az = a[2]
+            for j in range(i + 1, n):
+                b = pts[j]
+                dx = b[0] - ax
+                dy = b[1] - ay
+                dz = b[2] - az
+                s = dx * dx + dy * dy + dz * dz
+                if s >= d2:
+                    continue
+                if s < 1e-22:
+                    dx = rng.gauss(0.0, 1e-6)
+                    dy = rng.gauss(0.0, 1e-6)
+                    dz = rng.gauss(0.0, 1e-6)
+                    s = dx * dx + dy * dy + dz * dz
+                    if s < 1e-30:
+                        continue
+                dist = math.sqrt(s)
+                gap = d - dist
+                if gap > worst:
+                    worst = gap
+                f = 0.5 * gap / dist
+                ux = dx * f
+                uy = dy * f
+                uz = dz * f
+                bx = b[0] + ux
+                by = b[1] + uy
+                bz = b[2] + uz
+                t = math.sqrt(bx * bx + by * by + bz * bz)
+                if t > 1e-12:
+                    b[0] = bx / t
+                    b[1] = by / t
+                    b[2] = bz / t
+                ax -= ux
+                ay -= uy
+                az -= uz
+                t = math.sqrt(ax * ax + ay * ay + az * az)
+                if t > 1e-12:
+                    ax /= t
+                    ay /= t
+                    az /= t
+            a[0] = ax
+            a[1] = ay
+            a[2] = az
+        if worst <= 1e-15:
+            return True
+    return False
 
 
-def place(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]:
-    rng = random.Random(seed)
-    pts = _fibonacci(n)
-    for q in pts:
-        q[0] += rng.gauss(0, 1e-3); q[1] += rng.gauss(0, 1e-3); q[2] += rng.gauss(0, 1e-3)
-        r = math.sqrt(q[0] ** 2 + q[1] ** 2 + q[2] ** 2)
-        q[0] /= r; q[1] /= r; q[2] /= r
-    best, best_d = [list(q) for q in pts], _min_dist(pts)
-    start = time.perf_counter()
-    total = 0.85 * time_budget
-    powers = (8.0, 16.0, 32.0, 64.0)
-    for k, p in enumerate(powers):
-        deadline = start + total * (k + 1) / len(powers)
-        s = _min_dist(pts)  # scale so the largest term is about 1 and nothing overflows
-        e, f = _energy_and_forces(pts, p, s)
-        lr = 0.05 * s / max(1.0, max(math.sqrt(fx * fx + fy * fy + fz * fz) for fx, fy, fz in f))
-        while time.perf_counter() < deadline:
-            trial = _step(pts, f, lr)
-            e2, f2 = _energy_and_forces(trial, p, s)
-            if e2 < e:
-                pts, e, f = trial, e2, f2
-                lr *= 1.2
-                d = _min_dist(pts)
-                if d > best_d:
-                    best, best_d = [list(q) for q in pts], d
-            else:
-                lr *= 0.5
-                if lr < 1e-14:
+def _ratchet(pts, n, m, t_end, rng, step, sweeps, floor):
+    while step > floor:
+        if time.perf_counter() >= t_end:
+            break
+        d = m * (1.0 + step)
+        save = [p[:] for p in pts]
+        _separate(pts, n, d, sweeps, rng)
+        m2 = _min_dist(pts, n)
+        if m2 > m:
+            m = m2
+            step *= 1.25
+        else:
+            for i in range(n):
+                pts[i][0] = save[i][0]
+                pts[i][1] = save[i][1]
+                pts[i][2] = save[i][2]
+            step *= 0.55
+    return m
+
+
+def _hole(pts, n, skip, rng, tries):
+    """Point of the sphere furthest from every point except `skip`."""
+    bx = by = bz = 0.0
+    bd = -1.0
+    for _ in range(tries):
+        x = rng.gauss(0.0, 1.0)
+        y = rng.gauss(0.0, 1.0)
+        z = rng.gauss(0.0, 1.0)
+        s = math.sqrt(x * x + y * y + z * z)
+        if s < 1e-9:
+            continue
+        x /= s
+        y /= s
+        z /= s
+        worst = 8.0
+        for j in range(n):
+            if j == skip:
+                continue
+            p = pts[j]
+            dx = p[0] - x
+            dy = p[1] - y
+            dz = p[2] - z
+            e = dx * dx + dy * dy + dz * dz
+            if e < worst:
+                worst = e
+                if worst < bd:
                     break
-    return [tuple(q) for q in best]
+        if worst > bd:
+            bd = worst
+            bx, by, bz = x, y, z
+    return bx, by, bz
 
+
+def place(n, time_budget, seed):
+    t0 = time.perf_counter()
+    t_end = t0 + time_budget * 0.90
+    rng = random.Random((seed & 0xFFFFFFFF) * 1000003 + n * 7717)
+
+    if n == 1:
+        return [(0.0, 0.0, 1.0)]
+    if n == 2:
+        return [(0.0, 0.0, 1.0), (0.0, 0.0, -1.0)]
+
+    sweeps = 14 if n <= 60 else 10
+    tries = 220
+
+    best = None
+    best_m = -1.0
+    starts = 0
+    offs = (0.5, 0.0, 0.36, 0.72, 0.18)
+    t_multi = t0 + time_budget * 0.45
+    while True:
+        now = time.perf_counter()
+        if starts and now >= t_multi:
+            break
+        slice_end = min(t_end, now + max(0.2, time_budget * 0.10))
+        if starts < len(offs):
+            pts = _fib(n, offs[starts])
+        elif starts % 3:
+            pts = _fib(n, rng.random())
+            j = 0.25 / math.sqrt(n)
+            for p in pts:
+                p[0] += rng.gauss(0.0, j)
+                p[1] += rng.gauss(0.0, j)
+                p[2] += rng.gauss(0.0, j)
+                t = math.sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2])
+                p[0] /= t
+                p[1] /= t
+                p[2] /= t
+        else:
+            pts = _rand(n, rng)
+        m = _min_dist(pts, n)
+        if m < 1e-9:
+            m = 1e-6
+        m = _ratchet(pts, n, m, slice_end, rng, 0.02, sweeps, 1e-8)
+        if m > best_m:
+            best_m = m
+            best = [p[:] for p in pts]
+        starts += 1
+
+    t_moves = t0 + time_budget * 0.82
+    cap = time_budget * 0.07
+    while time.perf_counter() < t_moves:
+        pts = [p[:] for p in best]
+        i, j = _min_pair(pts, n)
+        k = i if rng.random() < 0.5 else j
+        if rng.random() < 0.65:
+            hx, hy, hz = _hole(pts, n, k, rng, tries)
+            pts[k][0] = hx
+            pts[k][1] = hy
+            pts[k][2] = hz
+        else:
+            amp = best_m * rng.uniform(0.15, 0.6)
+            p = pts[k]
+            p[0] += rng.gauss(0.0, amp)
+            p[1] += rng.gauss(0.0, amp)
+            p[2] += rng.gauss(0.0, amp)
+            t = math.sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2])
+            if t > 1e-12:
+                p[0] /= t
+                p[1] /= t
+                p[2] /= t
+        m = _min_dist(pts, n)
+        if m < 1e-9:
+            continue
+        m = _ratchet(pts, n, m, min(t_moves, time.perf_counter() + cap),
+                     rng, 0.004, sweeps, 1e-8)
+        if m > best_m:
+            best_m = m
+            best = [p[:] for p in pts]
+
+    _ratchet(best, n, best_m, t_end, rng, 1e-5, sweeps, 1e-14)
+    return [(p[0], p[1], p[2]) for p in best]
+