Coulomb energy is a sum, not a minimum, so the packing ratchet does not apply; heavy-ball momentum plus an adaptive step converts the same few hundred force passes into convergence.

-"""Baseline: Fibonacci-sphere start, then projected gradient descent on the Coulomb energy with an
-adaptive step. Lands within about 0.1-1% of the records. Beat it."""
+"""Thomson: n points on the unit sphere minimising the Coulomb energy sum 1/|x_i - x_j|.
 
+Unlike Tammes, this is not a packing problem in disguise: the objective is smooth and every
+pair contributes, so the separation ratchet I use elsewhere does not apply and this is a
+local-optimisation problem with a basin-hopping shell.
+
+The binding constraint is arithmetic, not search.  One energy-and-force pass over n=122 is
+about 7400 pairs, which in standard-library Python leaves only a few hundred passes inside an
+8 second budget.  Plain projected gradient descent, which is what the baseline does, spends
+them all crawling.  Two changes buy the convergence instead:
+
+  * heavy-ball momentum, so the step compounds along the valley rather than re-deriving it;
+  * an adaptive step that grows while the energy falls and halves with a velocity reset when
+    it does not, so no pass is wasted on an overshoot that has to be undone.
+
+Starts are Fibonacci spirals at several offsets; the best is then refined for the remaining
+budget, with occasional single-point kicks to escape a shallow basin.
+"""
+
 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 _energy_and_forces(pts: list[list[float]]) -> tuple[float, list[list[float]]]:
-    n = len(pts)
+def _energy_forces(P, n):
+    """Coulomb energy and the descent direction (the repulsive force) in one pass."""
+    F = [0.0] * (3 * n)
     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]
+        i3 = 3 * i
+        ax = P[i3]
+        ay = P[i3 + 1]
+        az = P[i3 + 2]
+        fx = F[i3]
+        fy = F[i3 + 1]
+        fz = F[i3 + 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
-            inv = 1.0 / math.sqrt(d2)
+            j3 = 3 * j
+            dx = ax - P[j3]
+            dy = ay - P[j3 + 1]
+            dz = az - P[j3 + 2]
+            s = dx * dx + dy * dy + dz * dz
+            if s < 1e-24:
+                s = 1e-24
+            r = math.sqrt(s)
+            inv = 1.0 / r
             e += inv
-            g = inv / d2  # 1/d^3
-            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
+            g = inv / s                      # 1 / r^3
+            tx = dx * g
+            ty = dy * g
+            tz = dz * g
+            fx += tx
+            fy += ty
+            fz += tz
+            F[j3] -= tx
+            F[j3 + 1] -= ty
+            F[j3 + 2] -= tz
+        F[i3] = fx
+        F[i3 + 1] = fy
+        F[i3 + 2] = fz
+    return e, F
 
 
-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  # drop the radial part: move along the sphere
-        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 _energy(P, n):
+    e = 0.0
+    for i in range(n):
+        i3 = 3 * i
+        ax = P[i3]
+        ay = P[i3 + 1]
+        az = P[i3 + 2]
+        for j in range(i + 1, n):
+            j3 = 3 * j
+            dx = ax - P[j3]
+            dy = ay - P[j3 + 1]
+            dz = az - P[j3 + 2]
+            s = dx * dx + dy * dy + dz * dz
+            if s < 1e-24:
+                s = 1e-24
+            e += 1.0 / math.sqrt(s)
+    return e
 
 
-def place(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]:
-    rng = random.Random(seed)
-    pts = _fibonacci(n)
-    for p in pts:  # tiny jitter so symmetric starts do not sit on a saddle
-        p[0] += rng.gauss(0, 1e-3); p[1] += rng.gauss(0, 1e-3); p[2] += rng.gauss(0, 1e-3)
-        r = math.sqrt(p[0] ** 2 + p[1] ** 2 + p[2] ** 2)
-        p[0] /= r; p[1] /= r; p[2] /= r
-    deadline = time.perf_counter() + 0.85 * time_budget
-    lr = 0.5 / (n * math.sqrt(n))
-    e, f = _energy_and_forces(pts)
-    while time.perf_counter() < deadline:
-        trial = _step(pts, f, lr)
-        e2, f2 = _energy_and_forces(trial)
+def _relax(P, n, t_end, step, beta):
+    """Heavy-ball descent on the sphere with an adaptive step."""
+    V = [0.0] * (3 * n)
+    e, F = _energy_forces(P, n)
+    best = P[:]
+    best_e = e
+    while time.perf_counter() < t_end:
+        save = P[:]
+        for k in range(3 * n):
+            V[k] = beta * V[k] + F[k]
+        for i in range(n):
+            i3 = 3 * i
+            x = P[i3] + step * V[i3]
+            y = P[i3 + 1] + step * V[i3 + 1]
+            z = P[i3 + 2] + step * V[i3 + 2]
+            t = math.sqrt(x * x + y * y + z * z)
+            if t < 1e-12:
+                x, y, z, t = save[i3], save[i3 + 1], save[i3 + 2], 1.0
+            P[i3] = x / t
+            P[i3 + 1] = y / t
+            P[i3 + 2] = z / t
+        e2, F2 = _energy_forces(P, n)
         if e2 < e:
-            pts, e, f = trial, e2, f2
-            lr *= 1.2
+            e = e2
+            F = F2
+            step *= 1.08
+            if e2 < best_e:
+                best_e = e2
+                best = P[:]
         else:
-            lr *= 0.5
-            if lr < 1e-12:
+            P[:] = save
+            for k in range(3 * n):
+                V[k] = 0.0
+            step *= 0.5
+            if step < 1e-15:
                 break
-    return [tuple(p) for p in pts]
+    P[:] = best
+    return best_e
 
+
+def place(n, time_budget, seed):
+    t0 = time.perf_counter()
+    t_end = t0 + time_budget * 0.90
+    rng = random.Random((seed & 0xFFFFFFFF) * 1000003 + n * 6151)
+
+    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)]
+
+    step0 = 0.05 / math.sqrt(n)
+    offs = (0.5, 0.0, 0.3, 0.7)
+
+    best = None
+    best_e = None
+    t_multi = t0 + time_budget * 0.40
+    for k, off in enumerate(offs):
+        now = time.perf_counter()
+        if k and now >= t_multi:
+            break
+        pts = _fib(n, off)
+        P = []
+        for p in pts:
+            P.extend(p)
+        e = _relax(P, n, min(t_end, now + max(0.3, time_budget * 0.11)), step0, 0.8)
+        if best_e is None or e < best_e:
+            best_e = e
+            best = P[:]
+
+    # refine the incumbent, with occasional single-point kicks out of a shallow basin
+    while time.perf_counter() < t_end - 0.05:
+        slice_end = min(t_end, time.perf_counter() + time_budget * 0.14)
+        P = best[:]
+        if best_e is not None and rng.random() < 0.5:
+            i = rng.randrange(n)
+            i3 = 3 * i
+            a = 0.35
+            x = P[i3] + rng.gauss(0.0, a)
+            y = P[i3 + 1] + rng.gauss(0.0, a)
+            z = P[i3 + 2] + rng.gauss(0.0, a)
+            t = math.sqrt(x * x + y * y + z * z)
+            if t > 1e-12:
+                P[i3] = x / t
+                P[i3 + 1] = y / t
+                P[i3 + 2] = z / t
+        e = _relax(P, n, slice_end, step0, 0.9)
+        if e < best_e:
+            best_e = e
+            best = P[:]
+
+    return [(best[3 * i], best[3 * i + 1], best[3 * i + 2]) for i in range(n)]
+