{"id":"morse-clusters-rho6","name":"Morse clusters at rho = 6, minimum energy","family":"cluster-optimisation","description":"Place N atoms in 3-D to minimise the Morse energy sum e^(6(1-r))(e^(6(1-r)) - 2), for thirteen N between 20 and 80. Scored against the Cambridge Cluster Database putative global minima at rho = 6; 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 = 6, 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 = 6,   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 = 6` has the same curvature at the well bottom as\nLennard-Jones, so the landscape is a close cousin of the LJ one (mostly icosahedral, with a Marks\ndecahedron winning at `n = 75`) but the global minima are different structures. Thirteen values of\n`n` between 20 and 80 are the benchmark. The database authors describe locating every Morse\nglobal minimum up to `n = 80` as \"a significant achievement and one which no unbiased global\noptimization algorithm has yet managed\".\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 = 6` 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 = 6`\ncolumn is the global minimum for that `n`). Every entry was recomputed from the database's own\npoints file (for 26A and 31B, whose files hold the `rho = 4` geometry, after a local relaxation at\n`rho = 6`) and agrees to within 5e-7. None is proven optimal. An energy below a record by more\nthan 1e-5 is listed in `records_beaten`.\n\n| n | record E | database label | point group | nearest-neighbour pairs |\n|---|---|---|---|---|\n| 20 | -72.507782 | 20A | C2v | 72 |\n| 26 | -100.549598 | 26A | Td | 102 |\n| 31 | -122.857743 | 31B | Cs | 126 |\n| 38 | -157.477108 | 38E | C5v | 147 |\n| 46 | -199.177751 | 46B | C2v | 186 |\n| 55 | -250.286609 | 55B | Ih | 234 |\n| 61 | -278.726626 | 61B | C2v | 260 |\n| 65 | -298.392345 | 65C | C2 | 278 |\n| 69 | -319.819905 | 69D | Cs | 299 |\n| 72 | -336.121753 | 72C | Cs | 313 |\n| 75 | -351.472365 | 75C | D5h | 319 |\n| 76 | -356.372708 | 76D | C1 | 333 |\n| 80 | -378.333471 | 80B | Cs | 354 |\n\nNote the difference from Lennard-Jones: at `rho = 6` the `n = 38` global minimum is icosahedral\n(38E, C5v), not the fcc truncated octahedron (38D, -157.406902 here), while `n = 75` is the Marks\ndecahedron (75C) with fewer nearest-neighbour contacts than the icosahedral competitor. The\n`n = 76` global minimum is a low-symmetry C1 structure.\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.\n- Structure follows the range: at `rho = 6` the winners are mostly Mackay icosahedra with\n  anti-Mackay or Mackay overlayers, plus decahedra where strain matters (`n = 75`). Grow an\n  icosahedral core and fill the outer shell greedily; try a Marks decahedron for `n = 75`.\n- Surface moves beat random kicks: remove the atom with the highest pair-energy sum and re-place\n  it at the lowest-energy surface site, then re-minimise.\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).\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 = 3` contribute less than 1e-5 each\n  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-rho6. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe Morse energy at rho = 6 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,26,31,38,46,55,61,65,69,72,75,76,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,26,31,38,46,55,61,65,69,72,75,76,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 = 6.0\n\n# Lowest known Morse energy E = sum_{i<j} x_ij (x_ij - 2), x_ij = exp(rho (1 - r_ij)), rho = 6, in units of\n# the pair well depth and equilibrium pair separation, exactly as printed (6 decimals, the bold rho = 6\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 (26A and 31B, whose files are the\n# rho = 4 geometry, after a local relaxation at rho = 6) and agrees to within 5e-7. All are putative.\n# An energy more than BEAT_TOL below a record is flagged in records_beaten; the table is rounded to\n# 1e-6, so BEAT_TOL sits well above that rounding.\nRECORDS = {\n    20: -72.507782,   # 20A  C2v\n    26: -100.549598,  # 26A  Td\n    31: -122.857743,  # 31B  Cs\n    38: -157.477108,  # 38E  C5v\n    46: -199.177751,  # 46B  C2v\n    55: -250.286609,  # 55B  Ih\n    61: -278.726626,  # 61B  C2v\n    65: -298.392345,  # 65C  C2\n    69: -319.819905,  # 69D  Cs\n    72: -336.121753,  # 72C  Cs\n    75: -351.472365,  # 75C  D5h\n    76: -356.372708,  # 76D  C1\n    80: -378.333471,  # 80B  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\"morse6|{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 = 6 Morse energy with an adaptive step.\nA single local minimum: lands within a few per cent of the records. Beat it.\"\"\"\n\nimport math\nimport random\nimport time\n\nRHO = 6.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.03), y + rng.gauss(0, 0.03), z + rng.gauss(0, 0.03)] 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"}}