Lennard-Jones Clusters: Finding Molecular Ground States

Put N neutral atoms in a box and let them settle into their lowest energy configuration. Sounds simple. It’s NP-hard.

The number of local minima grows exponentially with the number of atoms. For a 13-atom cluster, there are thousands of distinct local minima. For larger clusters, exhaustive search is impossible. This is the Lennard-Jones cluster problem, and it’s the first real test case for ParGA.

The Physics

The Lennard-Jones potential describes how two neutral atoms interact:

$$V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

The first term is repulsion at short range (atoms don’t like overlapping). The second term is attraction at longer range (van der Waals forces). The potential has a minimum at r = 2^(1/6)σ, which is the equilibrium bond distance.

For a cluster of N atoms, the total energy is the sum over all pairs:

def lennard_jones_energy(positions: np.ndarray, n_atoms: int) -> float:
    coords = positions.reshape(n_atoms, 3)
    energy = 0.0
    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            r2 = np.sum((coords[i] - coords[j]) ** 2)
            if r2 < 1e-10:
                return 1e10  # Prevent division by zero
            r6 = r2**3
            r12 = r6**2
            energy += 4.0 * (1.0 / r12 - 1.0 / r6)
    return energy

We use reduced units where ε = 1 and σ = 1, which is standard in the computational chemistry literature.

Why It’s Hard

A 2-atom cluster is trivial: the minimum is just two atoms at the equilibrium distance, energy -1.0.

A 3-atom cluster forms an equilateral triangle, energy -3.0.

A 4-atom cluster is a tetrahedron, energy -6.0.

But as N grows, the landscape becomes increasingly rugged. There are many arrangements that look similar but have different energies. Small perturbations can push you from one local minimum to another.

The famous 13-atom cluster forms an icosahedron with one atom at the center. Its energy is -44.327. But there are other arrangements that are nearly as good and many that are substantially worse. A gradient descent starting from a random configuration will, more often than not, find one of the bad ones.

Known global minima from the Cambridge Cluster Database:

These are reference values for comparing optimization algorithms. Getting close to them with a general-purpose optimizer is a reasonable achievement. Matching them exactly usually requires problem-specific techniques.

The Genome Representation

Each atom has three coordinates (x, y, z), so an N-atom cluster has 3N genes:

genome_length = 3 * n_atoms  # [x1, y1, z1, x2, y2, z2, ...]

The bounds are set based on expected cluster size:

bound = 2.0 * (n_atoms ** (1/3))

This scales roughly with the cluster radius, giving atoms room to spread out initially while keeping the search space manageable.

Running the Optimization

With ParGA:

from parga import minimize

def fitness(positions):
    return lennard_jones_energy_vectorized(positions, n_atoms)

result = minimize(
    fitness,
    genome_length=3 * n_atoms,
    bounds=(-bound, bound),
    population_size=100 + 20 * n_atoms,
    generations=200 + 50 * n_atoms,
    mutation_rate=0.05,  # Higher than default
    crossover_rate=0.8,
    seed=42,
)

best_energy = -result.best_fitness

A few things to note:

Higher mutation rate: The default 1% mutation is too low for this problem. The landscape is rugged enough that we need more exploration. 5% works better.

Scaling with problem size: Bigger clusters need larger populations and more generations. The formula above grows linearly, which is conservative but reliable.

Vectorized energy: For clusters larger than 5 atoms, a vectorized energy calculation is noticeably faster:

def lennard_jones_energy_vectorized(positions: np.ndarray, n_atoms: int) -> float:
    coords = positions.reshape(n_atoms, 3)

    # All pairwise distances at once
    diff = coords[:, np.newaxis, :] - coords[np.newaxis, :, :]
    r2 = np.sum(diff**2, axis=2)

    # Upper triangle only (i < j)
    upper_indices = np.triu_indices(n_atoms, k=1)
    r2_pairs = r2[upper_indices]
    r2_pairs = np.maximum(r2_pairs, 1e-10)

    r6 = r2_pairs**3
    r12 = r6**2
    return np.sum(4.0 * (1.0 / r12 - 1.0 / r6))

NumPy broadcasting handles all the pairwise computations in a single pass.

Results, and where they stop

The draft version of this post benchmarked clusters of three to seven atoms, reported errors in the noise, and concluded that ParGA solves this. Running the same code out to thirteen atoms tells a more useful story:

Up to eight atoms it is essentially solved. By ten it is missing by 5%. By thirteen it is missing by 31%, which is not a near-miss at all: it has found a completely different and much worse structure, after four minutes of work.

Two things are happening at once, and both are worth seeing.

The search space wins

The number of distinct local minima on a Lennard-Jones surface grows roughly exponentially in the number of atoms. Estimates put it near 1,500 for N=13 and around 10^5 by N=20. A population of a few hundred exploring a space with thousands of basins is not sampling it, it is spot-checking it.

That is why the failure is abrupt rather than gradual. Below the point where the population can cover the basins, the GA finds the right one nearly every time. Above it, it finds a plausible-looking wrong one, reports a converged run, and gives you no signal that anything went wrong. The 31% at N=13 arrived without a warning.

Thirteen is worse than it looks

N=13 is not just the largest case here, it is a famously hard one. The global minimum is the icosahedron, a structurally distinct motif with a much deeper energy well than its neighbours, and it is not reachable by small perturbations of the structures that work for N=12. The search has to find a different shape, not a better arrangement of the same shape.

Landscapes like that are exactly where crossover should help, since combining two partial structures is how a GA is supposed to make a discontinuous jump. It did not manage it here with default settings and a few hundred individuals, which is a fair statement about the budget rather than about the method.

What the time column says

Time grows faster than the error does: 0.5 seconds at N=4, 16 at N=7, 231 at N=13. The energy evaluation itself is O(N²) pairwise terms, which is nothing, so almost all of that growth is the search needing more generations and larger populations to get anywhere.

That combination, cost climbing while quality falls, is the honest signature of a method running out of road. It is much more useful to know where that happens than to know that the easy cases work.

This is the landscape where islands earn their keep

On the standard benchmark functions, splitting a population into islands at a fixed budget makes results monotonically worse. Rastrigin at sixteen islands is fifteen times worse than one population spending the same evaluations.

Lennard-Jones inverts that completely. Same total budget, 400 individuals for 300 generations, median best energy over 8 seeds:

More islands is better, monotonically, and at thirteen atoms the eight-island run recovers three full units of energy over the single population for exactly the same spend.

The difference is what the local minima are made of. Rastrigin’s are cosine ripples on a bowl: shallow, uniform, and all pointing at the same global structure, so isolation just slows down the spread of a good answer. A Lennard-Jones minimum is a structural motif, a distinct arrangement of atoms with its own deep well, and the good arrangement is not reachable by nudging a bad one. That is precisely the situation islands were invented for. Keeping subpopulations apart lets several different structures develop far enough to be evaluated on their merits instead of being swallowed by whichever one got lucky first.

So the advice “use islands for multimodal problems” is not wrong, it is underspecified. Counting local minima does not tell you which case you are in. What matters is whether the basins represent genuinely different solutions or just texture on the way to the same one.

What Makes LJ Clusters Good for GAs

Several properties make this a good genetic algorithm problem:

No gradient required: You don’t need to differentiate through the energy function. The GA works entirely from energy evaluations.

Multimodal landscape: Many local minima means gradient descent fails frequently. The GA’s parallel exploration pays off.

Continuous representation: Real-valued genomes fit naturally. Crossover and mutation have intuitive geometric interpretations.

Known ground truth: Having reference energies lets you evaluate algorithm quality objectively.

Scalable difficulty: You can make the problem harder by adding atoms. The difficulty grows smoothly, not discretely.

Limitations

Genetic algorithms aren’t the state of the art for Lennard-Jones clusters. Basin-hopping and related techniques (which combine local optimization with stochastic jumps) outperform pure GAs on this problem.

The GA finds good configurations, not provably optimal ones. For small clusters where we know the exact answer, the GA gets very close. For larger clusters where the true minimum is unknown, you can’t be sure.

The scaling is also a concern. A 13-atom cluster is 39 genes. A 50-atom cluster is 150 genes. GAs slow down substantially as genome length increases, and the exponential growth of local minima makes the problem harder faster than the GA gets better.

What LJ Clusters Teach Us

This problem demonstrates several general principles:

Fitness function flexibility: The energy function is a simple formula, but it could equally well be a molecular dynamics simulation or a quantum chemistry calculation. The GA doesn’t care.

Exploration vs exploitation: High mutation rates help on rugged landscapes. The default parameters assume smoother problems.

Island models for reliability: When single-run variance is high, multiple populations help. The cost is roughly equivalent to running multiple independent trials, but migration shares information.

Problem-specific tuning matters: The settings that work for benchmark functions (Rastrigin, Sphere) aren’t necessarily right for real problems. LJ clusters benefit from more mutation and larger populations.

Next: protein folding on a lattice, where the search space is discrete and the physics is simplified but the problem is still NP-hard.