{"id":"morse-clusters-rho14","name":"Morse clusters at rho = 14, minimum energy","family":"cluster-optimisation","description":"Place N atoms in 3-D to minimise the Morse energy sum e^(14(1-r))(e^(14(1-r)) - 2), the short-ranged case, for thirteen N between 20 and 80. Scored against the Cambridge Cluster Database putative global minima at rho = 14; none is proven optimal.","metric":"record_ratio","direction":"maximize","tolerance":0.05,"eval_timeout_seconds":300,"agent_timeout_seconds":1800,"mutable":["cluster.py"],"runtime":"python>=3.11, standard library only (math, random, itertools, functools, collections, heapq, time)","decomposable":true,"status":"active","captain":null,"parent_problem":null,"program_md":"# Morse clusters at rho = 14, minimum energy\n\n## Goal\n\n`cluster.py` exposes `cluster(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]`:\n`n` atoms `(x, y, z)` anywhere in 3-D space. Minimise the Morse energy\n\n    E = sum over i < j of  x_ij (x_ij - 2),   x_ij = exp(rho (1 - r_ij)),   rho = 14,   r_ij = |x_i - x_j|\n\nin reduced units (pair well depth 1, equilibrium pair separation 1). The range parameter `rho`\nsets how short-ranged the potential is; `rho = 14` is the hardest column of the Cambridge Cluster\nDatabase's Morse table (`morse-clusters-rho6` is the Lennard-Jones-like sibling). With such a narrow\nwell, strain is expensive and nearest-neighbour count is nearly everything: the global minima are\nclose-packed (fcc and hcp) or decahedral fragments rather than icosahedra, the energy landscape is\nfar rougher, and the number of local minima explodes. Thirteen values of `n` between 20 and 80 are\nthe benchmark. The database authors describe locating every Morse global minimum up to `n = 80`\nas \"a significant achievement and one which no unbiased global optimization algorithm has yet\nmanaged\", and the short-ranged clusters are the ones that defeat most methods.\n\n## Metric\n\n    metric = mean over n in NS of  E(n) / record(n)\n\nBoth energies are negative, so 1.0 means matching every record and above 1.0 means beating at\nleast one (a lower, more negative energy gives a larger ratio). A cluster with `E >= 0` scores 0\nfor that `n`. The eval recomputes the energy itself from the returned coordinates; a wrong count, a\nnon-finite coordinate, two coincident atoms (`r < 1e-9`) or a non-finite energy is a failed run.\nNothing your solver reports is used. `ZT_EVAL_SEED` only changes the `seed` handed to `cluster`,\nso your method must be robust to its starting point.\n\n## Records\n\nBest-known energies at `rho = 14` from the Cambridge Cluster Database \"Morse Clusters: Table of\nGlobal Minima\" (J. P. K. Doye, D. J. Wales and R. S. Berry, J. Chem. Phys. 103, 4234 (1995);\nJ. P. K. Doye and D. J. Wales, J. Chem. Soc. Faraday Trans. 93, 4233 (1997), with later additions\ncredited on the page; https://www-wales.ch.cam.ac.uk/~wales/CCD/jon/structures/Morse/tables.html,\nfetched 2026-09-07), exactly as printed there (6 decimals; the entry in bold in the `rho = 14`\ncolumn is the global minimum for that `n`). Every entry was recomputed from the database's own\npoints file: the files store the geometry at the `rho` given on their first line (6 or 10 for these\nentries, 14 for 30D'), so each was relaxed locally at `rho = 14` (ramping `rho` up in steps of 1,\nthen L-BFGS to a gradient below 1e-5), and every relaxed energy agrees with the printed value to\nwithin 5e-7. None is proven optimal. An energy below a record by more than 1e-5 is listed in\n`records_beaten`.\n\n| n | record E | database label | point group | nearest-neighbour pairs | rho of points file |\n|---|---|---|---|---|---|\n| 20 | -64.791953 | 20B | C2v | 65 | 6 |\n| 23 | -77.302495 | 23C | Cs | 78 | 10 |\n| 30 | -106.835790 | 30D' | Cs | 107 | 14 |\n| 38 | -144.321054 | 38D | Oh | 144 | 6 |\n| 45 | -174.511633 | 45E | C2v | 175 | 6 |\n| 50 | -198.455633 | 50D | D3h | 198 | 6 |\n| 55 | -220.646208 | 55C | C2v | 221 | 6 |\n| 59 | -240.572493 | 59E | Td | 240 | 6 |\n| 64 | -264.587042 | 64C | C2v | 265 | 6 |\n| 68 | -282.683003 | 68H | C3v | 282 | 6 |\n| 72 | -302.413229 | 72E | C1 | 303 | 6 |\n| 75 | -318.407330 | 75C | D5h | 319 | 6 |\n| 80 | -340.811371 | 80G | Cs | 340 | 6 |\n\nStructural notes from the table: `n = 38` is the fcc truncated octahedron (38D, Oh), `n = 50` a\nD3h close-packed fragment, `n = 59` a Td fragment, `n = 75` the Marks decahedron (75C, D5h),\n`n = 80` a Cs close-packed structure; the low-symmetry entries (`n = 72`, C1) are close-packed with\nstacking faults. 30D' is a later addition to the page, not credited there by name (the intro\ncredits the other primed entries to Branz, to Leary and Chan, and to Locatelli and Schoen). The\n`rho = 6` file of 32E does not relax into the printed `rho = 14` minimum, so `n = 32` is excluded.\n\n## Constraints\n\n- Standard library only. No numpy, no scipy. The eval rejects other imports.\n- Respect `time_budget` (seconds, per call). The eval fails a call that runs past 1.25x + 3 s.\n- Deterministic given `seed`: use `random.Random(seed)`.\n\n## Iterating\n\n- `ZT_EVAL_NS=38,75` runs a subset of n; `ZT_EVAL_PER_N_SECONDS=2` shortens the per-n budget. The\n  default is all thirteen n at 7 s each (about 95 s). Run `python eval.py` in your workspace.\n- The per-n detail in the eval output shows which n are furthest from their record.\n\n## Ideas that are known to matter (check the journal before repeating one)\n\n- One descent from a lattice cut lands a few per cent above the record. Basin hopping (perturb,\n  locally minimise, Metropolis accept) with a properly converged local minimiser is what found\n  these entries; the database authors used exactly that, and at `rho = 14` they needed seeding\n  from structural families, not just random kicks.\n- Structure follows the range: at `rho = 14` strain is prohibitive, so the winners are fcc and hcp\n  fragments (truncated octahedra, twinned close-packed pieces with stacking faults) and Marks\n  decahedra. Enumerate close-packed fragments by nearest-neighbour count (`nn` in the table), then\n  relax; icosahedral starts are usually wrong here.\n- Surface moves beat random kicks: remove the atom with the lowest number of contacts and re-place\n  it at the site with the most, then re-minimise. With a well this narrow, moves that keep the\n  contact count are the only ones that survive.\n- The ratio only rounds to 1.0 once you are within about 1e-8 relative, so finish each minimum\n  with a well-converged local optimiser (conjugate gradient or L-BFGS). The Hessian is very stiff\n  at `rho = 14`: steepest descent needs a tiny step, so a quasi-Newton polish pays off much more\n  than at `rho = 6`.\n- The pure-Python pair loop is the bottleneck: `exp` is the expensive call, so cache `x` per pair\n  and compute energy and forces in one pass. Pairs beyond `r = 2` contribute less than 1e-6 each\n  at `rho = 14` and can be skipped during minimisation, but keep the full sum for the final energy.\n\nWrite one honest line in `NOTES.md`: the idea, and which `n` it helped.\n\nSimpler is better: all else equal prefer the shorter solver, and treat removing code for an\nequal score as a win. Log every experiment, including discards, in your results.tsv.\n","eval_py":"\"\"\"Eval for morse-clusters-rho14. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe Morse energy at rho = 14 is recomputed here from the coordinates the solver returns; nothing\nthe solver reports is trusted.\n\nEnv:\n  ZT_EVAL_SEED             seed handed to cluster()\n  ZT_EVAL_PER_N_SECONDS    time budget per n (default 7)\n  ZT_EVAL_NS               comma-separated n values (default \"20,23,30,38,45,50,55,59,64,68,72,75,80\")\n\"\"\"\n\nfrom __future__ import annotations\n\nimport ast\nimport json\nimport math\nimport os\nimport random\nimport sys\nimport time\nfrom pathlib import Path\n\nSEED = os.environ.get(\"ZT_EVAL_SEED\", \"dev-seed\")\nBUDGET = float(os.environ.get(\"ZT_EVAL_PER_N_SECONDS\", \"7\"))\nNS = [int(x) for x in os.environ.get(\"ZT_EVAL_NS\", \"20,23,30,38,45,50,55,59,64,68,72,75,80\").split(\",\")]\nSTDLIB_ALLOW = {\"math\", \"random\", \"itertools\", \"functools\", \"collections\", \"heapq\", \"time\", \"sys\", \"typing\", \"operator\"}\nEPS = 1e-9\nFORBIDDEN_NAMES = {\"__import__\", \"importlib\", \"builtins\", \"__builtins__\", \"open\", \"exec\", \"eval\", \"compile\",\n                   \"globals\", \"__loader__\", \"__spec__\", \"breakpoint\", \"input\", \"memoryview\", \"vars\"}\n\nRHO = 14.0\n\n# Lowest known Morse energy E = sum_{i<j} x_ij (x_ij - 2), x_ij = exp(rho (1 - r_ij)), rho = 14, in units of\n# the pair well depth and equilibrium pair separation, exactly as printed (6 decimals, the bold rho = 14\n# entry per N) in the Cambridge Cluster Database \"Morse Clusters: Table of Global Minima\" (J. P. K. Doye,\n# D. J. Wales and R. S. Berry, J. Chem. Phys. 103, 4234 (1995); J. P. K. Doye and D. J. Wales, J. Chem.\n# Soc. Faraday Trans. 93, 4233 (1997), with later additions credited on the page;\n# https://www-wales.ch.cam.ac.uk/~wales/CCD/jon/structures/Morse/tables.html, fetched 2026-09-07).\n# Every entry was recomputed here from the database points file after a local relaxation at rho = 14\n# (the files hold the rho = 6 or rho = 10 geometry of the same basin, except 30D' which is stored at\n# rho = 14) and agrees to within 5e-7. All are putative. An energy more than BEAT_TOL below a record is\n# flagged in records_beaten; the table is rounded to 1e-6, so BEAT_TOL sits well above that rounding.\nRECORDS = {\n    20: -64.791953,  # 20B  C2v\n    23: -77.302495,  # 23C  Cs\n    30: -106.835790, # 30D' Cs\n    38: -144.321054, # 38D  Oh\n    45: -174.511633, # 45E  C2v\n    50: -198.455633, # 50D  D3h\n    55: -220.646208, # 55C  C2v\n    59: -240.572493, # 59E  Td\n    64: -264.587042, # 64C  C2v\n    68: -282.683003, # 68H  C3v\n    72: -302.413229, # 72E  C1\n    75: -318.407330, # 75C  D5h\n    80: -340.811371, # 80G  Cs\n}\nBEAT_TOL = 1e-5\n\n\ndef fail(msg: str, kind: str = \"error\") -> None:\n    print(json.dumps({\"metric\": 0.0, \"error\": msg, \"kind\": kind}))\n    sys.exit(1)\n\n\ndef check_imports(path: Path) -> None:\n    try:\n        tree = ast.parse(path.read_text(encoding=\"utf-8\"))\n    except SyntaxError as e:\n        fail(f\"syntax error in pack.py: {e}\", \"compile_error\")\n    for node in ast.walk(tree):\n        names = []\n        if isinstance(node, ast.Import):\n            names = [a.name.split(\".\")[0] for a in node.names]\n        elif isinstance(node, ast.ImportFrom) and node.module:\n            names = [node.module.split(\".\")[0]]\n        for nm in names:\n            if nm not in STDLIB_ALLOW:\n                fail(f\"import of '{nm}' is not allowed (stdlib subset only: {sorted(STDLIB_ALLOW)})\", \"compile_error\")\n        # dynamic imports and raw file/process access are not part of the problem either\n        ident = node.id if isinstance(node, ast.Name) else node.attr if isinstance(node, ast.Attribute) else None\n        if ident in FORBIDDEN_NAMES:\n            fail(f\"use of '{ident}' is not allowed in a solver\", \"compile_error\")\n        if isinstance(node, ast.ImportFrom) and node.level:\n            fail(\"relative imports are not allowed in a solver\", \"compile_error\")\n\n\ndef atoms(points: list, n: int) -> list[tuple[float, float, float]]:\n    \"\"\"n finite 3-D positions. Anything else is a wrong answer.\"\"\"\n    if not isinstance(points, (list, tuple)) or len(points) != n:\n        fail(f\"cluster({n}) must return {n} atoms\", \"wrong_answer\")\n    out = []\n    for p in points:\n        try:\n            x, y, z = float(p[0]), float(p[1]), float(p[2])\n        except Exception:\n            fail(f\"cluster({n}) returned a non-point {p!r}\", \"wrong_answer\")\n        if not all(map(math.isfinite, (x, y, z))):\n            fail(f\"cluster({n}) returned a non-finite coordinate\", \"wrong_answer\")\n        out.append((x, y, z))\n    return out\n\n\ndef energy(pts: list[tuple[float, float, float]], n: int) -> float:\n    e = 0.0\n    for i in range(n):\n        xi, yi, zi = pts[i]\n        for j in range(i + 1, n):\n            xj, yj, zj = pts[j]\n            r = math.sqrt((xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2)\n            if r < EPS:\n                fail(f\"cluster({n}): atoms {i} and {j} coincide\", \"wrong_answer\")\n            x = math.exp(RHO * (1.0 - r))\n            e += x * (x - 2.0)\n    if not math.isfinite(e):\n        fail(f\"cluster({n}): energy is not finite\", \"wrong_answer\")\n    return e\n\n\ndef main() -> None:\n    here = Path(__file__).parent\n    check_imports(here / \"cluster.py\")\n    sys.path.insert(0, str(here))\n    try:\n        import cluster as cand  # noqa: E402\n    except SystemExit:\n        raise\n    except Exception as e:\n        fail(f\"import cluster.py failed: {e!r}\", \"compile_error\")\n    if not hasattr(cand, \"cluster\"):\n        fail(\"cluster.py must define cluster(n, time_budget, seed)\", \"compile_error\")\n\n    seed_int = random.Random(f\"morse14|{SEED}\").getrandbits(32)\n    per_n, beaten = {}, []\n    for n in NS:\n        if n not in RECORDS:\n            fail(f\"no record for n={n}\", \"error\")\n        t0 = time.perf_counter()\n        try:\n            points = cand.cluster(n, BUDGET, seed_int)\n        except SystemExit:\n            raise\n        except Exception as e:\n            fail(f\"cluster({n}) raised {e!r}\", \"runtime_error\")\n        elapsed = time.perf_counter() - t0\n        if elapsed > 1.25 * BUDGET + 3:\n            fail(f\"cluster({n}) took {elapsed:.1f}s against a {BUDGET:.0f}s budget\", \"timeout\")\n        e = energy(atoms(points, n), n)\n        # Both energies are negative: ratio = E / record is 1.0 at the record, larger when E is lower\n        # (more negative), and 0 for a non-binding cluster (E >= 0), which keeps it in [0, inf).\n        ratio = e / RECORDS[n] if e < 0 else 0.0\n        per_n[n] = {\"energy\": round(e, 10), \"record\": RECORDS[n], \"ratio\": round(ratio, 9),\n                    \"seconds\": round(elapsed, 2)}\n        if e < RECORDS[n] - BEAT_TOL:\n            beaten.append(n)\n    metric = sum(v[\"ratio\"] for v in per_n.values()) / len(per_n)\n    print(json.dumps({\"metric\": round(metric, 9), \"per_n\": per_n, \"records_beaten\": beaten}))\n\n\nif __name__ == \"__main__\":\n    main()\n","baseline":{"cluster.py":"\"\"\"Baseline: the n sites of an fcc lattice (nearest-neighbour spacing 1, the Morse pair minimum) closest\nto a random centre, jittered, then steepest descent on the rho = 14 Morse energy with an adaptive step.\nA single local minimum: lands within a few per cent of the records (the rho = 14 well is narrow, so the\njitter is kept small). Beat it.\"\"\"\n\nimport math\nimport random\nimport time\n\nRHO = 14.0\n\n\ndef _fcc_start(n: int, rng: random.Random) -> list[list[float]]:\n    a = math.sqrt(2.0)  # cubic cell edge so that nearest neighbours sit at 1\n    m = int(math.ceil((n / 4.0) ** (1 / 3))) + 2\n    cx, cy, cz = (rng.uniform(-0.3, 0.3) * a for _ in range(3))  # random centre breaks the symmetry\n    sites = []\n    for i in range(-m, m + 1):\n        for j in range(-m, m + 1):\n            for k in range(-m, m + 1):\n                for dx, dy, dz in ((0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)):\n                    x, y, z = (i + dx) * a, (j + dy) * a, (k + dz) * a\n                    sites.append(((x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2, x, y, z))\n    sites.sort()\n    return [[x + rng.gauss(0, 0.02), y + rng.gauss(0, 0.02), z + rng.gauss(0, 0.02)] for _, x, y, z in sites[:n]]\n\n\ndef _energy_and_forces(pts: list[list[float]]) -> tuple[float, list[list[float]]]:\n    n = len(pts)\n    e = 0.0\n    f = [[0.0, 0.0, 0.0] for _ in range(n)]\n    for i in range(n):\n        xi, yi, zi = pts[i]\n        fi = f[i]\n        for j in range(i + 1, n):\n            xj, yj, zj = pts[j]\n            dx, dy, dz = xi - xj, yi - yj, zi - zj\n            r = math.sqrt(dx * dx + dy * dy + dz * dz)\n            x = math.exp(RHO * (1.0 - r))\n            e += x * (x - 2.0)\n            g = 2.0 * RHO * x * (x - 1.0) / r  # -dV/dr / r\n            fi[0] += dx * g; fi[1] += dy * g; fi[2] += dz * g\n            fj = f[j]\n            fj[0] -= dx * g; fj[1] -= dy * g; fj[2] -= dz * g\n    return e, f\n\n\ndef _step(pts: list[list[float]], f: list[list[float]], lr: float) -> list[list[float]]:\n    return [[x + lr * fx, y + lr * fy, z + lr * fz] for (x, y, z), (fx, fy, fz) in zip(pts, f)]\n\n\ndef cluster(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]:\n    rng = random.Random(seed)\n    pts = _fcc_start(n, rng)\n    deadline = time.perf_counter() + 0.85 * time_budget\n    lr = 1e-3\n    e, f = _energy_and_forces(pts)\n    while time.perf_counter() < deadline:\n        trial = _step(pts, f, lr)\n        e2, f2 = _energy_and_forces(trial)\n        if e2 < e:\n            pts, e, f = trial, e2, f2\n            lr *= 1.2\n        else:\n            lr *= 0.5\n            if lr < 1e-12:\n                break\n    return [tuple(p) for p in pts]\n"}}