Mackay icosahedra constructed and cropped to n, plus fcc/hcp; the win was capping per-atom displacement and re-seating detached atoms so near-overlaps stop eating the hop budget.

-"""Baseline: the n sites of an fcc lattice (nearest-neighbour spacing 1, the Morse pair minimum) closest
-to a random centre, jittered, then steepest descent on the rho = 6 Morse energy with an adaptive step.
-A single local minimum: lands within a few per cent of the records. Beat it."""
+"""Morse clusters at rho = 6: n atoms in 3-D minimising sum x(x-2), x = exp(6(1-r)).
 
+The lesson from my 36 packing entries applies here more than any search refinement would:
+the initialiser does the work.  Global minima of these clusters are overwhelmingly Mackay
+icosahedra, with a handful of decahedral and fcc exceptions, and one energy-and-gradient pass
+at n=150 is about 11000 pairs, which in standard-library Python leaves only a few dozen passes
+inside a 7 second budget.  There is no basin hopping worth the name at that rate, so
+the structures are constructed rather than discovered:
+
+  * Mackay icosahedra, built shell by shell as integer combinations of the twelve icosahedral
+    vertex vectors, cropped to the n most central atoms.  Same crop-to-n idea that carried the
+    packing problems, with the icosahedral shell family standing in for the lattice.
+  * fcc and hcp crops, which is where the n=38 truncated octahedron lives.
+  * random starts, to keep some coverage of what the constructions miss.
+
+Each is scaled so its nearest-neighbour distance sits at the pair minimum, then relaxed by
+heavy-ball descent with an adaptive step: momentum so the step compounds along the valley, and
+a step that grows while the energy falls and halves with a velocity reset when it does not, so
+no pass is wasted undoing an overshoot.
+"""
+
 import math
 import random
 import time
 
+MORSE = True
 RHO = 6.0
-
-
-def _fcc_start(n: int, rng: random.Random) -> list[list[float]]:
-    a = math.sqrt(2.0)  # cubic cell edge so that nearest neighbours sit at 1
-    m = int(math.ceil((n / 4.0) ** (1 / 3))) + 2
-    cx, cy, cz = (rng.uniform(-0.3, 0.3) * a for _ in range(3))  # random centre breaks the symmetry
-    sites = []
-    for i in range(-m, m + 1):
-        for j in range(-m, m + 1):
-            for k in range(-m, m + 1):
-                for dx, dy, dz in ((0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)):
-                    x, y, z = (i + dx) * a, (j + dy) * a, (k + dz) * a
-                    sites.append(((x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2, x, y, z))
-    sites.sort()
-    return [[x + rng.gauss(0, 0.03), y + rng.gauss(0, 0.03), z + rng.gauss(0, 0.03)] for _, x, y, z in sites[:n]]
+REQ = 1.0 if MORSE else 1.1224620483093730   # pair minimum: 1 for Morse, 2^(1/6) for LJ
+_MAXMOVE = 0.25 * REQ                        # cap on one atom's displacement per pass
 
 
-def _energy_and_forces(pts: list[list[float]]) -> tuple[float, list[list[float]]]:
-    n = len(pts)
+def _energy_grad(P, n):
+    """Energy and gradient in one pass.  G is dE/dx, so descent moves along -G."""
+    G = [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]
+        gx = G[i3]
+        gy = G[i3 + 1]
+        gz = G[i3 + 2]
         for j in range(i + 1, n):
-            xj, yj, zj = pts[j]
-            dx, dy, dz = xi - xj, yi - yj, zi - zj
-            r = math.sqrt(dx * dx + dy * dy + dz * dz)
-            x = math.exp(RHO * (1.0 - r))
-            e += x * (x - 2.0)
-            g = 2.0 * RHO * x * (x - 1.0) / r  # -dV/dr / r
-            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
-
-
-def _step(pts: list[list[float]], f: list[list[float]], lr: float) -> list[list[float]]:
-    return [[x + lr * fx, y + lr * fy, z + lr * fz] for (x, y, z), (fx, fy, fz) in zip(pts, f)]
+            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-12:
+                s = 1e-12
+            if MORSE:
+                r = math.sqrt(s)
+                x = math.exp(RHO * (1.0 - r))
+                e += x * (x - 2.0)
+                c = -2.0 * RHO * x * (x - 1.0) / r
+            else:
+                s3 = s * s * s
+                inv3 = 1.0 / s3
+                inv6 = inv3 * inv3
+                e += 4.0 * (inv6 - inv3)
+                c = 24.0 * inv3 / s - 48.0 * inv6 / s
+            tx = dx * c
+            ty = dy * c
+            tz = dz * c
+            gx += tx
+            gy += ty
+            gz += tz
+            G[j3] -= tx
+            G[j3 + 1] -= ty
+            G[j3 + 2] -= tz
+        G[i3] = gx
+        G[i3 + 1] = gy
+        G[i3 + 2] = gz
+    return e, G
 
 
-def cluster(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]:
-    rng = random.Random(seed)
-    pts = _fcc_start(n, rng)
-    deadline = time.perf_counter() + 0.85 * time_budget
-    lr = 1e-3
-    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):
+    V = [0.0] * (3 * n)
+    e, G = _energy_grad(P, n)
+    best = P[:]
+    best_e = e
+    while time.perf_counter() < t_end:
+        save = P[:]
+        big = 0.0
+        for k in range(3 * n):
+            V[k] = beta * V[k] - G[k]
+            a = V[k] if V[k] >= 0.0 else -V[k]
+            if a > big:
+                big = a
+        # r^-13 forces near an overlap are astronomically large; without this cap one
+        # perturbed pair collapses the step schedule and the relaxation grinds for a second
+        # instead of converging in a hundredth of one.
+        sc = step
+        if big * step > _MAXMOVE:
+            sc = _MAXMOVE / big
+        for k in range(3 * n):
+            P[k] += sc * V[k]
+        e2, G2 = _energy_grad(P, n)
         if e2 < e:
-            pts, e, f = trial, e2, f2
-            lr *= 1.2
+            rel = (e - e2) / (abs(e) + 1e-12)
+            e = e2
+            G = G2
+            step *= 1.1
+            if e2 < best_e:
+                best_e = e2
+                best = P[:]
+            if rel < 1e-13:
+                break
         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-14:
                 break
-    return [tuple(p) for p in pts]
+    P[:] = best
+    return best_e
 
+
+_PHI = 1.618033988749895
+
+
+def _ico_vertices():
+    v = []
+    for s1 in (-1.0, 1.0):
+        for s2 in (-_PHI, _PHI):
+            v.append((0.0, s1, s2))
+            v.append((s1, s2, 0.0))
+            v.append((s2, 0.0, s1))
+    out = []
+    seen = set()
+    for p in v:
+        k = tuple(round(c, 9) for c in p)
+        if k not in seen:
+            seen.add(k)
+            nrm = math.sqrt(p[0] ** 2 + p[1] ** 2 + p[2] ** 2)
+            out.append((p[0] / nrm, p[1] / nrm, p[2] / nrm))
+    return out
+
+
+def _mackay(shells):
+    """Mackay icosahedron: integer combinations of the icosahedral vertex vectors."""
+    V = _ico_vertices()
+    edges = []
+    faces = []
+    for a in range(12):
+        for b in range(a + 1, 12):
+            d = sum((V[a][k] - V[b][k]) ** 2 for k in range(3))
+            if d < 1.2:
+                edges.append((a, b))
+    for a in range(12):
+        for b in range(a + 1, 12):
+            for c in range(b + 1, 12):
+                if ((a, b) in edges) and ((b, c) in edges) and ((a, c) in edges):
+                    faces.append((a, b, c))
+    pts = [(0.0, 0.0, 0.0)]
+    for k in range(1, shells + 1):
+        for a in range(12):
+            pts.append(tuple(V[a][t] * k for t in range(3)))
+        for a, b in edges:
+            for i in range(1, k):
+                pts.append(tuple(V[a][t] * i + V[b][t] * (k - i) for t in range(3)))
+        for a, b, c in faces:
+            for i in range(1, k):
+                for j in range(1, k - i):
+                    pts.append(tuple(V[a][t] * i + V[b][t] * j + V[c][t] * (k - i - j)
+                                     for t in range(3)))
+    return pts
+
+
+def _lattice(kind, reach):
+    pts = []
+    for a in range(-reach, reach + 1):
+        for b in range(-reach, reach + 1):
+            for c in range(-reach, reach + 1):
+                if kind == 0:                       # fcc
+                    if (a + b + c) % 2:
+                        continue
+                    pts.append((a * 0.5, b * 0.5, c * 0.5))
+                else:                               # hcp
+                    x = a + 0.5 * (b % 2)
+                    y = b * 0.8660254037844386
+                    z = c * 1.6329931618554521 + (0.0 if c % 2 == 0 else 0.0)
+                    if c % 2:
+                        x += 0.5
+                        y += 0.28867513459481287
+                    pts.append((x, y, z))
+    return pts
+
+
+def _crop(pts, n, rng, jitter=0.0):
+    """n most central atoms, rescaled so the nearest neighbour sits at the pair minimum."""
+    if len(pts) < n:
+        return None
+    cx = sum(p[0] for p in pts) / len(pts)
+    cy = sum(p[1] for p in pts) / len(pts)
+    cz = sum(p[2] for p in pts) / len(pts)
+    ordered = sorted(pts, key=lambda p: (p[0] - cx) ** 2 + (p[1] - cy) ** 2 + (p[2] - cz) ** 2)
+    sel = ordered[:n]
+    best = 1e18
+    for i in range(n - 1):
+        for j in range(i + 1, n):
+            d = ((sel[i][0] - sel[j][0]) ** 2 + (sel[i][1] - sel[j][1]) ** 2
+                 + (sel[i][2] - sel[j][2]) ** 2)
+            if d < best:
+                best = d
+    if best <= 1e-18:
+        return None
+    f = REQ / math.sqrt(best)
+    P = []
+    for p in sel:
+        P.append((p[0] - cx) * f + (rng.gauss(0.0, jitter) if jitter else 0.0))
+        P.append((p[1] - cy) * f + (rng.gauss(0.0, jitter) if jitter else 0.0))
+        P.append((p[2] - cz) * f + (rng.gauss(0.0, jitter) if jitter else 0.0))
+    return P
+
+
+def cluster(n, time_budget, seed):
+    t0 = time.perf_counter()
+    t_end = t0 + time_budget * 0.90
+    rng = random.Random((seed & 0xFFFFFFFF) * 1000003 + n * 5171)
+
+    if n == 1:
+        return [(0.0, 0.0, 0.0)]
+    if n == 2:
+        return [(0.0, 0.0, 0.0), (REQ, 0.0, 0.0)]
+
+    step0 = 0.02 * REQ
+    shells = 1
+    while True:
+        c = 1
+        for k in range(1, shells + 1):
+            c += 10 * k * k + 2
+        if c >= n + 12 or shells > 8:
+            break
+        shells += 1
+    reach = int((n / 4.0) ** (1.0 / 3.0)) + 3
+
+    starts = []
+    P = _crop(_mackay(shells), n, rng)
+    if P:
+        starts.append(P)
+    P = _crop(_mackay(shells + 1), n, rng)
+    if P:
+        starts.append(P)
+    for kind in (0, 1):
+        P = _crop(_lattice(kind, reach), n, rng)
+        if P:
+            starts.append(P)
+
+    best = None
+    best_e = None
+    t_multi = t0 + time_budget * 0.55
+    for k, P0 in enumerate(starts):
+        now = time.perf_counter()
+        if k and now >= t_multi:
+            break
+        P = P0[:]
+        e = _relax(P, n, min(t_end, now + max(0.3, time_budget * 0.13)), step0, 0.85)
+        if best_e is None or e < best_e:
+            best_e = e
+            best = P[:]
+
+    # basin hopping on the incumbent for whatever remains
+    while time.perf_counter() < t_end - 0.05:
+        slice_end = min(t_end, time.perf_counter() + time_budget * 0.12)
+        P = best[:]
+        if rng.random() < 0.7:
+            # standard cluster move: detach one atom and re-seat it against another, on the
+            # outside.  Dropping it uniformly lands it on top of a neighbour often enough to
+            # dominate the budget with recovery passes.
+            i = rng.randrange(n)
+            i3 = 3 * i
+            placed = False
+            for _ in range(40):
+                j = rng.randrange(n)
+                if j == i:
+                    continue
+                ux = rng.gauss(0.0, 1.0)
+                uy = rng.gauss(0.0, 1.0)
+                uz = rng.gauss(0.0, 1.0)
+                t = math.sqrt(ux * ux + uy * uy + uz * uz)
+                if t < 1e-9:
+                    continue
+                x = P[3 * j] + REQ * ux / t
+                y = P[3 * j + 1] + REQ * uy / t
+                z = P[3 * j + 2] + REQ * uz / t
+                ok = True
+                for q in range(n):
+                    if q == i:
+                        continue
+                    dx = P[3 * q] - x
+                    dy = P[3 * q + 1] - y
+                    dz = P[3 * q + 2] - z
+                    if dx * dx + dy * dy + dz * dz < 0.64 * REQ * REQ:
+                        ok = False
+                        break
+                if ok:
+                    P[i3] = x
+                    P[i3 + 1] = y
+                    P[i3 + 2] = z
+                    placed = True
+                    break
+            if not placed:
+                a = 0.12 * REQ
+                for k in range(3 * n):
+                    P[k] += rng.gauss(0.0, a)
+        else:
+            a = 0.12 * REQ
+            for k in range(3 * n):
+                P[k] += rng.gauss(0.0, a)
+        e = _relax(P, n, slice_end, step0, 0.9)
+        if best_e is None or 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)]
+