Morse clusters at rho = 6, minimum energy

# Morse clusters at rho = 6, minimum energy ## Goal `cluster.py` exposes `cluster(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]`: `n` atoms `(x, y, z)` anywhere in 3-D space. Minimise the Morse energy E = sum over i < j of x_ij (x_ij - 2), x_ij = exp(rho (1 - r_ij)), rho = 6, r_ij = |x_i - x_j| in reduced units (pair well depth 1, equilibrium pair separation 1). The range parameter `rho` sets how short-ranged the potential is; `rho = 6` has the same curvature at the well bottom as Lennard-Jones, so the landscape is a close cousin of the LJ one (mostly icosahedral, with a Marks decahedron winning at `n = 75`) but the global minima are different structures. Thirteen values of `n` between 20 and 80 are the benchmark. The database authors describe locating every Morse global minimum up to `n = 80` as "a significant achievement and one which no unbiased global optimization algorithm has yet managed". ## Metric metric = mean over n in NS of E(n) / record(n) Both energies are negative, so 1.0 means matching every record and above 1.0 means beating at least one (a lower, more negative energy gives a larger ratio). A cluster with `E >= 0` scores 0 for that `n`. The eval recomputes the energy itself from the returned coordinates; a wrong count, a non-finite coordinate, two coincident atoms (`r < 1e-9`) or a non-finite energy is a failed run. Nothing your solver reports is used. `ZT_EVAL_SEED` only changes the `seed` handed to `cluster`, so your method must be robust to its starting point. ## Records Best-known energies at `rho = 6` from the Cambridge Cluster Database "Morse Clusters: Table of Global Minima" (J. P. K. Doye, D. J. Wales and R. S. Berry, J. Chem. Phys. 103, 4234 (1995); J. P. K. Doye and D. J. Wales, J. Chem. Soc. Faraday Trans. 93, 4233 (1997), with later additions credited on the page; https://www-wales.ch.cam.ac.uk/~wales/CCD/jon/structures/Morse/tables.html, fetched 2026-09-07), exactly as printed there (6 decimals; the entry in bold in the `rho = 6` column is the global minimum for that `n`). Every entry was recomputed from the database's own points file (for 26A and 31B, whose files hold the `rho = 4` geometry, after a local relaxation at `rho = 6`) and agrees to within 5e-7. None is proven optimal. An energy below a record by more than 1e-5 is listed in `records_beaten`. | n | record E | database label | point group | nearest-neighbour pairs | |---|---|---|---|---| | 20 | -72.507782 | 20A | C2v | 72 | | 26 | -100.549598 | 26A | Td | 102 | | 31 | -122.857743 | 31B | Cs | 126 | | 38 | -157.477108 | 38E | C5v | 147 | | 46 | -199.177751 | 46B | C2v | 186 | | 55 | -250.286609 | 55B | Ih | 234 | | 61 | -278.726626 | 61B | C2v | 260 | | 65 | -298.392345 | 65C | C2 | 278 | | 69 | -319.819905 | 69D | Cs | 299 | | 72 | -336.121753 | 72C | Cs | 313 | | 75 | -351.472365 | 75C | D5h | 319 | | 76 | -356.372708 | 76D | C1 | 333 | | 80 | -378.333471 | 80B | Cs | 354 | Note the difference from Lennard-Jones: at `rho = 6` the `n = 38` global minimum is icosahedral (38E, C5v), not the fcc truncated octahedron (38D, -157.406902 here), while `n = 75` is the Marks decahedron (75C) with fewer nearest-neighbour contacts than the icosahedral competitor. The `n = 76` global minimum is a low-symmetry C1 structure. ## Constraints - Standard library only. No numpy, no scipy. The eval rejects other imports. - Respect `time_budget` (seconds, per call). The eval fails a call that runs past 1.25x + 3 s. - Deterministic given `seed`: use `random.Random(seed)`. ## Iterating - `ZT_EVAL_NS=38,75` runs a subset of n; `ZT_EVAL_PER_N_SECONDS=2` shortens the per-n budget. The default is all thirteen n at 7 s each (about 95 s). Run `python eval.py` in your workspace. - The per-n detail in the eval output shows which n are furthest from their record. ## Ideas that are known to matter (check the journal before repeating one) - One descent from a lattice cut lands a few per cent above the record. Basin hopping (perturb, locally minimise, Metropolis accept) with a properly converged local minimiser is what found these entries; the database authors used exactly that. - Structure follows the range: at `rho = 6` the winners are mostly Mackay icosahedra with anti-Mackay or Mackay overlayers, plus decahedra where strain matters (`n = 75`). Grow an icosahedral core and fill the outer shell greedily; try a Marks decahedron for `n = 75`. - Surface moves beat random kicks: remove the atom with the highest pair-energy sum and re-place it at the lowest-energy surface site, then re-minimise. - The ratio only rounds to 1.0 once you are within about 1e-8 relative, so finish each minimum with a well-converged local optimiser (conjugate gradient or L-BFGS). - The pure-Python pair loop is the bottleneck: `exp` is the expensive call, so cache `x` per pair and compute energy and forces in one pass. Pairs beyond `r = 3` contribute less than 1e-5 each and can be skipped during minimisation, but keep the full sum for the final energy. Write one honest line in `NOTES.md`: the idea, and which `n` it helped. Simpler is better: all else equal prefer the shorter solver, and treat removing code for an equal score as a win. Log every experiment, including discards, in your results.tsv.