{"id":"lennard-jones-clusters","name":"Lennard-Jones clusters, minimum energy","family":"cluster-optimisation","description":"Place N atoms in 3-D to minimise the Lennard-Jones energy sum 4(r^-12 - r^-6), for fifteen N between 20 and 150 including the non-icosahedral cases N = 38, 75-77, 98, 102-104. Scored against the Cambridge Cluster Database putative global minima; 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":"# Lennard-Jones clusters, 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 Lennard-Jones energy\n\n E = sum over i < j of 4 (r_ij^-12 - r_ij^-6), r_ij = |x_i - x_j|\n\nin reduced units (pair well depth 1, sigma 1, so the pair minimum sits at r = 2^(1/6)). Fifteen\nvalues of `n` between 20 and 150 are the benchmark. Lennard-Jones clusters are the standard test\nbed for global optimisation of configurational problems (basin hopping was introduced on them,\nWales and Doye 1997). Most putative global minima are Mackay icosahedra; the exceptions, an fcc\ntruncated octahedron at `n = 38`, Marks decahedra at `n = 75, 76, 77, 102, 103, 104` and a\ntetrahedral structure at `n = 98`, are the classic hard cases because their basins are tiny\ncompared with the icosahedral funnel. All eight are in this set.\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 from the Cambridge Cluster Database \"Table of Lennard-Jones Cluster Global\nMinima\" (D. J. Wales and J. P. K. Doye, J. Phys. Chem. A 101, 5111 (1997) and the references it\ncredits; https://www-wales.ch.cam.ac.uk/~wales/CCD/jon/structures/LJ/tables.150.html, fetched\n2026-09-07), exactly as printed there (6 decimals). Every entry was recomputed from the database's\nown points file and agrees to within 1.2e-6. None is proven optimal: no LJ_N global minimum with\n`n > 4` has a proof. An energy below a record by more than 1e-5 is listed in `records_beaten`.\n\n| n | record E | point group | structure | first reported |\n|---|---|---|---|---|\n| 20 | -77.177043 | C2v | icosahedral | Hoare |\n| 26 | -108.315616 | Td | icosahedral | Hoare |\n| 31 | -133.586422 | Cs | icosahedral | Northby |\n| 38 | -173.928427 | Oh | fcc truncated octahedron | Gomez, Pillardy, Doye |\n| 55 | -279.248470 | Ih | Mackay icosahedron | Hoare |\n| 69 | -359.882566 | C5v | icosahedral | Wales, Barron, Leary |\n| 75 | -397.492331 | D5h | Marks decahedron | Doye |\n| 76 | -402.894866 | Cs | Marks decahedron | Doye |\n| 77 | -409.083517 | C2v | Marks decahedron | Doye |\n| 98 | -543.665361 | Td | tetrahedral | Leary |\n| 102 | -569.363652 | C2v | Marks decahedron | Doye |\n| 103 | -575.766131 | Cs | Marks decahedron | Doye |\n| 104 | -582.086642 | C2v | Marks decahedron | Doye |\n| 110 | -621.788224 | Cs | icosahedral | Northby |\n| 150 | -893.310258 | C3v | icosahedral | Northby |\n\nFor the non-icosahedral `n` the same table lists the lowest icosahedral minimum, which is what an\nunbiased search usually finds first: `n = 38` -173.252378, `75` -396.282249, `76` -402.384580,\n`77` -408.518265, `98` -543.642957, `102` -569.277721, `103` -575.658879, `104` -582.038429.\nThose gaps (0.02 to 1.2 in energy, ratio 0.99996 to 0.997) are the signal in this benchmark.\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 fifteen n at 7 s each (about 110 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- The landscape has exponentially many minima; one descent from a lattice cut or a random start\n lands a few per cent above the record. Basin hopping (perturb, locally minimise, Metropolis accept\n at kT about 0.8) is what found most of these entries; a local minimiser that converges properly\n (L-BFGS or conjugate gradient in a few hundred pair-loop evaluations) is the core of it.\n- Seed with the right shape: a Mackay icosahedron with the outer shell filled by a greedy\n anti-Mackay/Mackay growth reaches the icosahedral minima directly, and `n = 55` and `n = 13`\n are complete icosahedra. For the hard `n`, start from a truncated octahedron (38), a Marks\n decahedron (75-77, 102-104) or the Leary tetrahedron (98); the icosahedral funnel will not\n find them in a 7 s budget.\n- Surface moves beat random kicks: remove the highest-energy atom and re-place it at the lowest\n 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.\n- The pure-Python pair loop is the bottleneck: cache differences, avoid function-call overhead,\n compute energy and forces in one pass, and skip pairs beyond a cutoff (about 2.5 sigma) during\n minimisation while keeping 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 lennard-jones-clusters. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe Lennard-Jones energy is recomputed here from the coordinates the solver returns; nothing the\nsolver 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,55,69,75,76,77,98,102,103,104,110,150\")\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,55,69,75,76,77,98,102,103,104,110,150\").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\n# Lowest known energy E = sum_{i 4 is proven. N = 38 (fcc truncated octahedron), 75-77 and\n# 102-104 (Marks decahedra) and 98 (tetrahedral) are the famous non-icosahedral cases; the rest are\n# Mackay-icosahedral. An energy more than BEAT_TOL below a record is flagged in records_beaten; the\n# table is rounded to 1e-6, so BEAT_TOL sits well above that rounding.\nRECORDS = {\n 20: -77.177043, # C2v Hoare\n 26: -108.315616, # Td Hoare\n 31: -133.586422, # Cs Northby\n 38: -173.928427, # Oh Gomez / Pillardy / Doye\n 55: -279.248470, # Ih Hoare\n 69: -359.882566, # C5v Wales / Barron / Leary\n 75: -397.492331, # D5h Doye\n 76: -402.894866, # Cs Doye\n 77: -409.083517, # C2v Doye\n 98: -543.665361, # Td Leary\n 102: -569.363652, # C2v Doye\n 103: -575.766131, # Cs Doye\n 104: -582.086642, # C2v Doye\n 110: -621.788224, # Cs Northby\n 150: -893.310258, # C3v Northby\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 r2 = (xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2\n if r2 < EPS * EPS:\n fail(f\"cluster({n}): atoms {i} and {j} coincide (infinite energy)\", \"wrong_answer\")\n r6 = 1.0 / (r2 * r2 * r2)\n e += 4.0 * (r6 * r6 - r6)\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\"lj|{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 at the LJ pair minimum) closest to\na random centre, jittered, then steepest descent on the LJ energy with an adaptive step. A single local\nminimum: lands within a few per cent of the records. Beat it.\"\"\"\n\nimport math\nimport random\nimport time\n\n\ndef _fcc_start(n: int, rng: random.Random) -> list[list[float]]:\n a = 2 ** (1 / 6) * math.sqrt(2.0) # cubic cell edge so that nearest neighbours sit at 2^(1/6)\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 r2 = dx * dx + dy * dy + dz * dz\n r6 = 1.0 / (r2 * r2 * r2)\n e += 4.0 * (r6 * r6 - r6)\n g = 24.0 * r6 * (2.0 * r6 - 1.0) / r2 # -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"}}